<https://www.gnu.org/licenses/>. */
#include "v_math.h"
-#include "poly_advsimd_f64.h"
static const struct data
{
- float64x2_t poly[12];
- float64x2_t pi, pi_over_2;
+ double c1, c3, c5, c7, c9, c11;
+ float64x2_t c0, c2, c4, c6, c8, c10;
uint64x2_t abs_mask;
+ float64x2_t pi, pi_over_2;
} data = {
/* Polynomial approximation of (asin(sqrt(x)) - sqrt(x)) / (x * sqrt(x))
on [ 0x1p-106, 0x1p-2 ], relative error: 0x1.c3d8e169p-57. */
- .poly = { V2 (0x1.555555555554ep-3), V2 (0x1.3333333337233p-4),
- V2 (0x1.6db6db67f6d9fp-5), V2 (0x1.f1c71fbd29fbbp-6),
- V2 (0x1.6e8b264d467d6p-6), V2 (0x1.1c5997c357e9dp-6),
- V2 (0x1.c86a22cd9389dp-7), V2 (0x1.856073c22ebbep-7),
- V2 (0x1.fd1151acb6bedp-8), V2 (0x1.087182f799c1dp-6),
- V2 (-0x1.6602748120927p-7), V2 (0x1.cfa0dd1f9478p-6), },
- .pi = V2 (0x1.921fb54442d18p+1),
- .pi_over_2 = V2 (0x1.921fb54442d18p+0),
+ .c0 = V2 (0x1.555555555554ep-3), .c1 = 0x1.3333333337233p-4,
+ .c2 = V2 (0x1.6db6db67f6d9fp-5), .c3 = 0x1.f1c71fbd29fbbp-6,
+ .c4 = V2 (0x1.6e8b264d467d6p-6), .c5 = 0x1.1c5997c357e9dp-6,
+ .c6 = V2 (0x1.c86a22cd9389dp-7), .c7 = 0x1.856073c22ebbep-7,
+ .c8 = V2 (0x1.fd1151acb6bedp-8), .c9 = 0x1.087182f799c1dp-6,
+ .c10 = V2 (-0x1.6602748120927p-7), .c11 = 0x1.cfa0dd1f9478p-6,
+ .pi = V2 (0x1.921fb54442d18p+1), .pi_over_2 = V2 (0x1.921fb54442d18p+0),
.abs_mask = V2 (0x7fffffffffffffff),
};
acos(x) ~ pi/2 - (x + x^3 P(x^2)).
- The largest observed error in this region is 1.18 ulps,
+ The largest observed error in this region is 1.18 ulp:
_ZGVnN2v_acos (0x1.fbab0a7c460f6p-2) got 0x1.0d54d1985c068p+0
want 0x1.0d54d1985c069p+0.
acos(x) = y + y * z * P(z), with z = (1-x)/2 and y = sqrt(z).
- The largest observed error in this region is 1.52 ulps,
- _ZGVnN2v_acos (0x1.23d362722f591p-1) got 0x1.edbbedf8a7d6ep-1
- want 0x1.edbbedf8a7d6cp-1. */
+ The largest observed error in this region is 1.50 ulp:
+ _ZGVnN2v_acos (0x1.252a2cf3fb9acp-1) got 0x1.ec1a46aa82901p-1
+ want 0x1.ec1a46aa829p-1. */
float64x2_t VPCS_ATTR V_NAME_D1 (acos) (float64x2_t x)
{
const struct data *d = ptr_barrier (&data);
float64x2_t z = vbslq_f64 (a_le_half, ax, vsqrtq_f64 (z2));
/* Use a single polynomial approximation P for both intervals. */
+ float64x2_t z3 = vmulq_f64 (z2, z);
float64x2_t z4 = vmulq_f64 (z2, z2);
float64x2_t z8 = vmulq_f64 (z4, z4);
- float64x2_t z16 = vmulq_f64 (z8, z8);
- float64x2_t p = v_estrin_11_f64 (z2, z4, z8, z16, d->poly);
- /* Finalize polynomial: z + z * z2 * P(z2). */
- p = vfmaq_f64 (z, vmulq_f64 (z, z2), p);
+ /* Order-11 Estrin. */
+ float64x2_t c13 = vld1q_f64 (&d->c1);
+ float64x2_t c57 = vld1q_f64 (&d->c5);
+ float64x2_t c911 = vld1q_f64 (&d->c9);
+
+ float64x2_t p01 = vfmaq_laneq_f64 (d->c0, z2, c13, 0);
+ float64x2_t p23 = vfmaq_laneq_f64 (d->c2, z2, c13, 1);
+ float64x2_t p03 = vfmaq_f64 (p01, z4, p23);
+
+ float64x2_t p45 = vfmaq_laneq_f64 (d->c4, z2, c57, 0);
+ float64x2_t p67 = vfmaq_laneq_f64 (d->c6, z2, c57, 1);
+ float64x2_t p47 = vfmaq_f64 (p45, z4, p67);
+
+ float64x2_t p89 = vfmaq_laneq_f64 (d->c8, z2, c911, 0);
+ float64x2_t p1011 = vfmaq_laneq_f64 (d->c10, z2, c911, 1);
+ float64x2_t p811 = vfmaq_f64 (p89, z4, p1011);
+
+ float64x2_t p411 = vfmaq_f64 (p47, z8, p811);
+ float64x2_t p = vfmaq_f64 (p03, z8, p411);
+
+ /* Finalize polynomial: z + z3 * P(z2). */
+ p = vfmaq_f64 (z, z3, p);
/* acos(|x|) = pi/2 - sign(x) * Q(|x|), for |x| < 0.5
= 2 Q(|x|) , for 0.5 < x < 1.0
<https://www.gnu.org/licenses/>. */
#include "sv_math.h"
-#include "poly_sve_f64.h"
static const struct data
{
- float64_t poly[12];
- float64_t pi, pi_over_2;
+ float64_t c1, c3, c5, c7, c9, c11;
+ float64_t c0, c2, c4, c6, c8, c10;
+ float64_t pi_over_2;
} data = {
/* Polynomial approximation of (asin(sqrt(x)) - sqrt(x)) / (x * sqrt(x))
on [ 0x1p-106, 0x1p-2 ], relative error: 0x1.c3d8e169p-57. */
- .poly = { 0x1.555555555554ep-3, 0x1.3333333337233p-4, 0x1.6db6db67f6d9fp-5,
- 0x1.f1c71fbd29fbbp-6, 0x1.6e8b264d467d6p-6, 0x1.1c5997c357e9dp-6,
- 0x1.c86a22cd9389dp-7, 0x1.856073c22ebbep-7, 0x1.fd1151acb6bedp-8,
- 0x1.087182f799c1dp-6, -0x1.6602748120927p-7, 0x1.cfa0dd1f9478p-6, },
- .pi = 0x1.921fb54442d18p+1,
+ .c0 = 0x1.555555555554ep-3, .c1 = 0x1.3333333337233p-4,
+ .c2 = 0x1.6db6db67f6d9fp-5, .c3 = 0x1.f1c71fbd29fbbp-6,
+ .c4 = 0x1.6e8b264d467d6p-6, .c5 = 0x1.1c5997c357e9dp-6,
+ .c6 = 0x1.c86a22cd9389dp-7, .c7 = 0x1.856073c22ebbep-7,
+ .c8 = 0x1.fd1151acb6bedp-8, .c9 = 0x1.087182f799c1dp-6,
+ .c10 = -0x1.6602748120927p-7, .c11 = 0x1.cfa0dd1f9478p-6,
.pi_over_2 = 0x1.921fb54442d18p+0,
};
acos(x) ~ pi/2 - (x + x^3 P(x^2)).
- The largest observed error in this region is 1.18 ulps,
- _ZGVsMxv_acos (0x1.fbc5fe28ee9e3p-2) got 0x1.0d4d0f55667f6p+0
- want 0x1.0d4d0f55667f7p+0.
+ The largest observed error in this region is 1.18 ulp:
+ _ZGVsMxv_acos (0x1.fbb7c9079b429p-2) got 0x1.0d51266607582p+0
+ want 0x1.0d51266607583p+0.
For |x| in [0.5, 1.0], use same approximation with a change of variable
acos(x) = y + y * z * P(z), with z = (1-x)/2 and y = sqrt(z).
- The largest observed error in this region is 1.52 ulps,
- _ZGVsMxv_acos (0x1.24024271a500ap-1) got 0x1.ed82df4243f0dp-1
- want 0x1.ed82df4243f0bp-1. */
+ The largest observed error in this region is 1.50 ulp:
+ _ZGVsMxv_acos (0x1.252a2cf3fb9acp-1) got 0x1.ec1a46aa82901p-1
+ want 0x1.ec1a46aa829p-1. */
svfloat64_t SV_NAME_D1 (acos) (svfloat64_t x, const svbool_t pg)
{
const struct data *d = ptr_barrier (&data);
+ svbool_t ptrue = svptrue_b64 ();
svuint64_t sign = svand_x (pg, svreinterpret_u64 (x), 0x8000000000000000);
svfloat64_t ax = svabs_x (pg, x);
svfloat64_t z = svsqrt_m (ax, a_gt_half, z2);
/* Use a single polynomial approximation P for both intervals. */
- svfloat64_t z4 = svmul_x (pg, z2, z2);
- svfloat64_t z8 = svmul_x (pg, z4, z4);
- svfloat64_t z16 = svmul_x (pg, z8, z8);
- svfloat64_t p = sv_estrin_11_f64_x (pg, z2, z4, z8, z16, d->poly);
+ svfloat64_t z3 = svmul_x (ptrue, z2, z);
+ svfloat64_t z4 = svmul_x (ptrue, z2, z2);
+ svfloat64_t z8 = svmul_x (ptrue, z4, z4);
+
+ svfloat64_t c13 = svld1rq (ptrue, &d->c1);
+ svfloat64_t c57 = svld1rq (ptrue, &d->c5);
+ svfloat64_t c911 = svld1rq (ptrue, &d->c9);
+
+ svfloat64_t p01 = svmla_lane (sv_f64 (d->c0), z2, c13, 0);
+ svfloat64_t p23 = svmla_lane (sv_f64 (d->c2), z2, c13, 1);
+ svfloat64_t p03 = svmla_x (pg, p01, z4, p23);
+
+ svfloat64_t p45 = svmla_lane (sv_f64 (d->c4), z2, c57, 0);
+ svfloat64_t p67 = svmla_lane (sv_f64 (d->c6), z2, c57, 1);
+ svfloat64_t p47 = svmla_x (pg, p45, z4, p67);
+
+ svfloat64_t p89 = svmla_lane (sv_f64 (d->c8), z2, c911, 0);
+ svfloat64_t p1011 = svmla_lane (sv_f64 (d->c10), z2, c911, 1);
+ svfloat64_t p811 = svmla_x (pg, p89, z4, p1011);
+
+ svfloat64_t p411 = svmla_x (pg, p47, z8, p811);
+ svfloat64_t p = svmad_x (pg, p411, z8, p03);
/* Finalize polynomial: z + z * z2 * P(z2). */
- p = svmla_x (pg, z, svmul_x (pg, z, z2), p);
+ p = svmad_x (pg, p, z3, z);
/* acos(|x|) = pi/2 - sign(x) * Q(|x|), for |x| < 0.5
= 2 Q(|x|) , for 0.5 < x < 1.0
= pi - 2 Q(|x|) , for -1.0 < x < -0.5. */
- svfloat64_t y
- = svreinterpret_f64 (svorr_x (pg, svreinterpret_u64 (p), sign));
-
- svbool_t is_neg = svcmplt (pg, x, 0.0);
- svfloat64_t off = svdup_f64_z (is_neg, d->pi);
- svfloat64_t mul = svsel (a_gt_half, sv_f64 (2.0), sv_f64 (-1.0));
- svfloat64_t add = svsel (a_gt_half, off, sv_f64 (d->pi_over_2));
-
- return svmla_x (pg, add, mul, y);
+ svfloat64_t mul = svreinterpret_f64 (
+ svlsl_m (a_gt_half, svreinterpret_u64 (sv_f64 (1.0)), 10));
+ mul = svreinterpret_f64 (sveor_x (ptrue, svreinterpret_u64 (mul), sign));
+ svfloat64_t add = svreinterpret_f64 (
+ svorr_x (ptrue, sign, svreinterpret_u64 (sv_f64 (d->pi_over_2))));
+ add = svsub_m (a_gt_half, sv_f64 (d->pi_over_2), add);
+
+ return svmsb_x (pg, p, mul, add);
}
<https://www.gnu.org/licenses/>. */
#include "v_math.h"
-#include "poly_advsimd_f64.h"
static const struct data
{
- float64x2_t poly[12];
+ float64x2_t c0, c2, c4, c6, c8, c10;
float64x2_t pi_over_2;
uint64x2_t abs_mask;
+ double c1, c3, c5, c7, c9, c11;
} data = {
/* Polynomial approximation of (asin(sqrt(x)) - sqrt(x)) / (x * sqrt(x))
on [ 0x1p-106, 0x1p-2 ], relative error: 0x1.c3d8e169p-57. */
- .poly = { V2 (0x1.555555555554ep-3), V2 (0x1.3333333337233p-4),
- V2 (0x1.6db6db67f6d9fp-5), V2 (0x1.f1c71fbd29fbbp-6),
- V2 (0x1.6e8b264d467d6p-6), V2 (0x1.1c5997c357e9dp-6),
- V2 (0x1.c86a22cd9389dp-7), V2 (0x1.856073c22ebbep-7),
- V2 (0x1.fd1151acb6bedp-8), V2 (0x1.087182f799c1dp-6),
- V2 (-0x1.6602748120927p-7), V2 (0x1.cfa0dd1f9478p-6), },
- .pi_over_2 = V2 (0x1.921fb54442d18p+0),
- .abs_mask = V2 (0x7fffffffffffffff),
+ .c0 = V2 (0x1.555555555554ep-3), .c1 = 0x1.3333333337233p-4,
+ .c2 = V2 (0x1.6db6db67f6d9fp-5), .c3 = 0x1.f1c71fbd29fbbp-6,
+ .c4 = V2 (0x1.6e8b264d467d6p-6), .c5 = 0x1.1c5997c357e9dp-6,
+ .c6 = V2 (0x1.c86a22cd9389dp-7), .c7 = 0x1.856073c22ebbep-7,
+ .c8 = V2 (0x1.fd1151acb6bedp-8), .c9 = 0x1.087182f799c1dp-6,
+ .c10 = V2 (-0x1.6602748120927p-7), .c11 = 0x1.cfa0dd1f9478p-6,
+ .pi_over_2 = V2 (0x1.921fb54442d18p+0), .abs_mask = V2 (0x7fffffffffffffff),
};
#define AllMask v_u64 (0xffffffffffffffff)
asin(x) = pi/2 - (y + y * z * P(z)), with z = (1-x)/2 and y = sqrt(z).
The largest observed error in this region is 2.69 ulps,
- _ZGVnN2v_asin (0x1.044ac9819f573p-1) got 0x1.110d7e85fdd5p-1
- want 0x1.110d7e85fdd53p-1. */
+ _ZGVnN2v_asin (0x1.044e8cefee301p-1) got 0x1.1111dd54ddf96p-1
+ want 0x1.1111dd54ddf99p-1. */
float64x2_t VPCS_ATTR V_NAME_D1 (asin) (float64x2_t x)
{
const struct data *d = ptr_barrier (&data);
return special_case (x, x, AllMask);
#endif
- uint64x2_t a_lt_half = vcltq_f64 (ax, v_f64 (0.5));
+ uint64x2_t a_lt_half = vcaltq_f64 (x, v_f64 (0.5));
/* Evaluate polynomial Q(x) = y + y * z * P(z) with
z = x ^ 2 and y = |x| , if |x| < 0.5
float64x2_t z4 = vmulq_f64 (z2, z2);
float64x2_t z8 = vmulq_f64 (z4, z4);
float64x2_t z16 = vmulq_f64 (z8, z8);
- float64x2_t p = v_estrin_11_f64 (z2, z4, z8, z16, d->poly);
+
+ /* order-11 estrin. */
+ float64x2_t c13 = vld1q_f64 (&d->c1);
+ float64x2_t c57 = vld1q_f64 (&d->c5);
+ float64x2_t c911 = vld1q_f64 (&d->c9);
+
+ float64x2_t p01 = vfmaq_laneq_f64 (d->c0, z2, c13, 0);
+ float64x2_t p23 = vfmaq_laneq_f64 (d->c2, z2, c13, 1);
+ float64x2_t p03 = vfmaq_f64 (p01, z4, p23);
+
+ float64x2_t p45 = vfmaq_laneq_f64 (d->c4, z2, c57, 0);
+ float64x2_t p67 = vfmaq_laneq_f64 (d->c6, z2, c57, 1);
+ float64x2_t p47 = vfmaq_f64 (p45, z4, p67);
+
+ float64x2_t p89 = vfmaq_laneq_f64 (d->c8, z2, c911, 0);
+ float64x2_t p1011 = vfmaq_laneq_f64 (d->c10, z2, c911, 1);
+ float64x2_t p811 = vfmaq_f64 (p89, z4, p1011);
+
+ float64x2_t p07 = vfmaq_f64 (p03, z8, p47);
+ float64x2_t p = vfmaq_f64 (p07, z16, p811);
/* Finalize polynomial: z + z * z2 * P(z2). */
p = vfmaq_f64 (z, vmulq_f64 (z, z2), p);
<https://www.gnu.org/licenses/>. */
#include "sv_math.h"
-#include "poly_sve_f64.h"
static const struct data
{
- float64_t poly[12];
- float64_t pi_over_2f;
+ float64_t c1, c3, c5, c7, c9, c11;
+ float64_t c0, c2, c4, c6, c8, c10;
+ float64_t pi_over_2;
} data = {
/* Polynomial approximation of (asin(sqrt(x)) - sqrt(x)) / (x * sqrt(x))
on [ 0x1p-106, 0x1p-2 ], relative error: 0x1.c3d8e169p-57. */
- .poly = { 0x1.555555555554ep-3, 0x1.3333333337233p-4,
- 0x1.6db6db67f6d9fp-5, 0x1.f1c71fbd29fbbp-6,
- 0x1.6e8b264d467d6p-6, 0x1.1c5997c357e9dp-6,
- 0x1.c86a22cd9389dp-7, 0x1.856073c22ebbep-7,
- 0x1.fd1151acb6bedp-8, 0x1.087182f799c1dp-6,
- -0x1.6602748120927p-7, 0x1.cfa0dd1f9478p-6, },
- .pi_over_2f = 0x1.921fb54442d18p+0,
+ .c0 = 0x1.555555555554ep-3, .c1 = 0x1.3333333337233p-4,
+ .c2 = 0x1.6db6db67f6d9fp-5, .c3 = 0x1.f1c71fbd29fbbp-6,
+ .c4 = 0x1.6e8b264d467d6p-6, .c5 = 0x1.1c5997c357e9dp-6,
+ .c6 = 0x1.c86a22cd9389dp-7, .c7 = 0x1.856073c22ebbep-7,
+ .c8 = 0x1.fd1151acb6bedp-8, .c9 = 0x1.087182f799c1dp-6,
+ .c10 = -0x1.6602748120927p-7, .c11 = 0x1.cfa0dd1f9478p-6,
+ .pi_over_2 = 0x1.921fb54442d18p+0,
};
-#define P(i) sv_f64 (d->poly[i])
-
/* Double-precision SVE implementation of vector asin(x).
For |x| in [0, 0.5], use an order 11 polynomial P such that the final
approximation is an odd polynomial: asin(x) ~ x + x^3 P(x^2).
- The largest observed error in this region is 0.52 ulps,
- _ZGVsMxv_asin(0x1.d95ae04998b6cp-2) got 0x1.ec13757305f27p-2
- want 0x1.ec13757305f26p-2.
-
- For |x| in [0.5, 1.0], use same approximation with a change of variable
+ The largest observed error in this region is 0.98 ulp:
+ _ZGVsMxv_asin (0x1.d98f6a748ed8ap-2) got 0x1.ec4eb661a73d3p-2
+ want 0x1.ec4eb661a73d2p-2.
- asin(x) = pi/2 - (y + y * z * P(z)), with z = (1-x)/2 and y = sqrt(z).
+ For |x| in [0.5, 1.0], use same approximation with a change of variable:
+ asin(x) = pi/2 - (y + y * z * P(z)), with z = (1-x)/2 and y = sqrt(z).
- The largest observed error in this region is 2.69 ulps,
- _ZGVsMxv_asin(0x1.044ac9819f573p-1) got 0x1.110d7e85fdd5p-1
- want 0x1.110d7e85fdd53p-1. */
+ The largest observed error in this region is 2.66 ulp:
+ _ZGVsMxv_asin (0x1.04024f6e2a2fbp-1) got 0x1.10b9586f087a8p-1
+ want 0x1.10b9586f087abp-1. */
svfloat64_t SV_NAME_D1 (asin) (svfloat64_t x, const svbool_t pg)
{
const struct data *d = ptr_barrier (&data);
+ svbool_t ptrue = svptrue_b64 ();
svuint64_t sign = svand_x (pg, svreinterpret_u64 (x), 0x8000000000000000);
svfloat64_t ax = svabs_x (pg, x);
svfloat64_t z = svsqrt_m (ax, a_ge_half, z2);
/* Use a single polynomial approximation P for both intervals. */
+ svfloat64_t z3 = svmul_x (pg, z2, z);
svfloat64_t z4 = svmul_x (pg, z2, z2);
svfloat64_t z8 = svmul_x (pg, z4, z4);
- svfloat64_t z16 = svmul_x (pg, z8, z8);
- svfloat64_t p = sv_estrin_11_f64_x (pg, z2, z4, z8, z16, d->poly);
+
+ svfloat64_t c13 = svld1rq (ptrue, &d->c1);
+ svfloat64_t c57 = svld1rq (ptrue, &d->c5);
+ svfloat64_t c911 = svld1rq (ptrue, &d->c9);
+
+ /* Order-11 Estrin scheme. */
+ svfloat64_t p01 = svmla_lane (sv_f64 (d->c0), z2, c13, 0);
+ svfloat64_t p23 = svmla_lane (sv_f64 (d->c2), z2, c13, 1);
+ svfloat64_t p03 = svmla_x (pg, p01, z4, p23);
+
+ svfloat64_t p45 = svmla_lane (sv_f64 (d->c4), z2, c57, 0);
+ svfloat64_t p67 = svmla_lane (sv_f64 (d->c6), z2, c57, 1);
+ svfloat64_t p47 = svmla_x (pg, p45, z4, p67);
+
+ svfloat64_t p89 = svmla_lane (sv_f64 (d->c8), z2, c911, 0);
+ svfloat64_t p1011 = svmla_lane (sv_f64 (d->c10), z2, c911, 1);
+ svfloat64_t p811 = svmla_x (pg, p89, z4, p1011);
+
+ svfloat64_t p411 = svmla_x (pg, p47, z8, p811);
+ svfloat64_t p = svmla_x (pg, p03, z8, p411);
+
/* Finalize polynomial: z + z * z2 * P(z2). */
- p = svmla_x (pg, z, svmul_x (pg, z, z2), p);
+ p = svmla_x (pg, z, z3, p);
- /* asin(|x|) = Q(|x|) , for |x| < 0.5
- = pi/2 - 2 Q(|x|), for |x| >= 0.5. */
- svfloat64_t y = svmad_m (a_ge_half, p, sv_f64 (-2.0), d->pi_over_2f);
+ /* asin(|x|) = Q(|x|), for |x| < 0.5
+ = pi/2 - 2 Q(|x|), for |x| >= 0.5. */
+ svfloat64_t y = svmad_m (a_ge_half, p, sv_f64 (-2.0), d->pi_over_2);
- /* Copy sign. */
+ /* Reinsert the sign from the argument. */
return svreinterpret_f64 (svorr_x (pg, svreinterpret_u64 (y), sign));
}
<https://www.gnu.org/licenses/>. */
#include "v_math.h"
-#include "poly_advsimd_f32.h"
static const struct data
{
- float32x4_t poly[5];
+ float32x4_t c0, c2, c4;
+ float c1, c3;
float32x4_t pi_over_2f;
} data = {
/* Polynomial approximation of (asin(sqrt(x)) - sqrt(x)) / (x * sqrt(x)) on
[ 0x1p-24 0x1p-2 ] order = 4 rel error: 0x1.00a23bbp-29 . */
- .poly = { V4 (0x1.55555ep-3), V4 (0x1.33261ap-4), V4 (0x1.70d7dcp-5),
- V4 (0x1.b059dp-6), V4 (0x1.3af7d8p-5) },
- .pi_over_2f = V4 (0x1.921fb6p+0f),
+ .c0 = V4 (0x1.55555ep-3f), .c1 = 0x1.33261ap-4f,
+ .c2 = V4 (0x1.70d7dcp-5f), .c3 = 0x1.b059dp-6f,
+ .c4 = V4 (0x1.3af7d8p-5f), .pi_over_2f = V4 (0x1.921fb6p+0f),
};
#define AbsMask 0x7fffffff
-#define Half 0x3f000000
#define One 0x3f800000
#define Small 0x39800000 /* 2^-12. */
/* Single-precision implementation of vector asin(x).
- For |x| < Small, approximate asin(x) by x. Small = 2^-12 for correct
- rounding. If WANT_SIMD_EXCEPT = 0, Small = 0 and we proceed with the
- following approximation.
- For |x| in [Small, 0.5], use order 4 polynomial P such that the final
+ For |x| <0.5, use order 4 polynomial P such that the final
approximation is an odd polynomial: asin(x) ~ x + x^3 P(x^2).
The largest observed error in this region is 0.83 ulps,
#endif
float32x4_t ax = vreinterpretq_f32_u32 (ia);
- uint32x4_t a_lt_half = vcltq_u32 (ia, v_u32 (Half));
+ uint32x4_t a_lt_half = vcaltq_f32 (x, v_f32 (0.5f));
/* Evaluate polynomial Q(x) = y + y * z * P(z) with
z = x ^ 2 and y = |x| , if |x| < 0.5
z = (1 - |x|) / 2 and y = sqrt(z), if |x| >= 0.5. */
float32x4_t z2 = vbslq_f32 (a_lt_half, vmulq_f32 (x, x),
- vfmsq_n_f32 (v_f32 (0.5), ax, 0.5));
+ vfmsq_n_f32 (v_f32 (0.5f), ax, 0.5f));
float32x4_t z = vbslq_f32 (a_lt_half, ax, vsqrtq_f32 (z2));
/* Use a single polynomial approximation P for both intervals. */
- float32x4_t p = v_horner_4_f32 (z2, d->poly);
+
+ /* PW Horner 3 evaluation scheme. */
+ float32x4_t z4 = vmulq_f32 (z2, z2);
+ float32x4_t c13 = vld1q_f32 (&d->c1);
+ float32x4_t p01 = vfmaq_laneq_f32 (d->c0, z2, c13, 0);
+ float32x4_t p23 = vfmaq_laneq_f32 (d->c2, z2, c13, 1);
+ float32x4_t p = vfmaq_f32 (p23, d->c4, z4);
+ p = vfmaq_f32 (p01, p, z4);
/* Finalize polynomial: z + z * z2 * P(z2). */
p = vfmaq_f32 (z, vmulq_f32 (z, z2), p);
/* asin(|x|) = Q(|x|) , for |x| < 0.5
= pi/2 - 2 Q(|x|), for |x| >= 0.5. */
float32x4_t y
- = vbslq_f32 (a_lt_half, p, vfmsq_n_f32 (d->pi_over_2f, p, 2.0));
+ = vbslq_f32 (a_lt_half, p, vfmsq_n_f32 (d->pi_over_2f, p, 2.0f));
/* Copy sign. */
return vbslq_f32 (v_u32 (AbsMask), y, x);
#include "math_config.h"
#include "v_math.h"
-#include "poly_advsimd_f64.h"
static const struct data
{
+ double c1, c3, c5, c7, c9, c11, c13, c15, c17, c19;
float64x2_t c0, c2, c4, c6, c8, c10, c12, c14, c16, c18;
float64x2_t pi_over_2;
- double c1, c3, c5, c7, c9, c11, c13, c15, c17, c19;
- uint64x2_t zeroinfnan, minustwo;
+ uint64x2_t zeroinfnan;
} data = {
- /* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
- [2**-1022, 1.0]. */
- .c0 = V2 (-0x1.5555555555555p-2),
- .c1 = 0x1.99999999996c1p-3,
- .c2 = V2 (-0x1.2492492478f88p-3),
- .c3 = 0x1.c71c71bc3951cp-4,
- .c4 = V2 (-0x1.745d160a7e368p-4),
- .c5 = 0x1.3b139b6a88ba1p-4,
- .c6 = V2 (-0x1.11100ee084227p-4),
- .c7 = 0x1.e1d0f9696f63bp-5,
- .c8 = V2 (-0x1.aebfe7b418581p-5),
- .c9 = 0x1.842dbe9b0d916p-5,
- .c10 = V2 (-0x1.5d30140ae5e99p-5),
- .c11 = 0x1.338e31eb2fbbcp-5,
- .c12 = V2 (-0x1.00e6eece7de8p-5),
- .c13 = 0x1.860897b29e5efp-6,
- .c14 = V2 (-0x1.0051381722a59p-6),
- .c15 = 0x1.14e9dc19a4a4ep-7,
- .c16 = V2 (-0x1.d0062b42fe3bfp-9),
- .c17 = 0x1.17739e210171ap-10,
- .c18 = V2 (-0x1.ab24da7be7402p-13),
- .c19 = 0x1.358851160a528p-16,
+ /* Coefficients of polynomial P such that
+ atan(x)~x+x*P(x^2) on [2^-1022, 1.0]. */
+ .c0 = V2 (-0x1.555555555552ap-2),
+ .c1 = 0x1.9999999995aebp-3,
+ .c2 = V2 (-0x1.24924923923f6p-3),
+ .c3 = 0x1.c71c7184288a2p-4,
+ .c4 = V2 (-0x1.745d11fb3d32bp-4),
+ .c5 = 0x1.3b136a18051b9p-4,
+ .c6 = V2 (-0x1.110e6d985f496p-4),
+ .c7 = 0x1.e1bcf7f08801dp-5,
+ .c8 = V2 (-0x1.ae644e28058c3p-5),
+ .c9 = 0x1.82eeb1fed85c6p-5,
+ .c10 = V2 (-0x1.59d7f901566cbp-5),
+ .c11 = 0x1.2c982855ab069p-5,
+ .c12 = V2 (-0x1.eb49592998177p-6),
+ .c13 = 0x1.69d8b396e3d38p-6,
+ .c14 = V2 (-0x1.ca980345c4204p-7),
+ .c15 = 0x1.dc050eafde0b3p-8,
+ .c16 = V2 (-0x1.7ea70755b8eccp-9),
+ .c17 = 0x1.ba3da3de903e8p-11,
+ .c18 = V2 (-0x1.44a4b059b6f67p-13),
+ .c19 = 0x1.c4a45029e5a91p-17,
.pi_over_2 = V2 (0x1.921fb54442d18p+0),
.zeroinfnan = V2 (2 * 0x7ff0000000000000ul - 1),
- .minustwo = V2 (0xc000000000000000),
};
#define SignMask v_u64 (0x8000000000000000)
}
/* Fast implementation of vector atan2.
- Maximum observed error is 2.8 ulps:
- _ZGVnN2vv_atan2 (0x1.9651a429a859ap+5, 0x1.953075f4ee26p+5)
- got 0x1.92d628ab678ccp-1
- want 0x1.92d628ab678cfp-1. */
+ Maximum observed error is 1.97 ulps:
+ _ZGVnN2vv_atan2 (0x1.42337dba73768p+5, 0x1.422d748cd3e29p+5)
+ got 0x1.9224810264efcp-1 want 0x1.9224810264efep-1. */
float64x2_t VPCS_ATTR V_NAME_D2 (atan2) (float64x2_t y, float64x2_t x)
{
const struct data *d = ptr_barrier (&data);
uint64x2_t pred_xlt0 = vcltzq_f64 (x);
uint64x2_t pred_aygtax = vcagtq_f64 (y, x);
- /* Set up z for call to atan. */
- float64x2_t n = vbslq_f64 (pred_aygtax, vnegq_f64 (ax), ay);
- float64x2_t q = vbslq_f64 (pred_aygtax, ay, ax);
- float64x2_t z = vdivq_f64 (n, q);
-
- /* Work out the correct shift. */
- float64x2_t shift
- = vreinterpretq_f64_u64 (vandq_u64 (pred_xlt0, d->minustwo));
- shift = vbslq_f64 (pred_aygtax, vaddq_f64 (shift, v_f64 (1.0)), shift);
- shift = vmulq_f64 (shift, d->pi_over_2);
-
- /* Calculate the polynomial approximation.
- Use split Estrin scheme for P(z^2) with deg(P)=19. Use split instead of
- full scheme to avoid underflow in x^16.
- The order 19 polynomial P approximates
- (atan(sqrt(x))-sqrt(x))/x^(3/2). */
+ /* Set up z for evaluation of atan. */
+ float64x2_t num = vbslq_f64 (pred_aygtax, vnegq_f64 (ax), ay);
+ float64x2_t den = vbslq_f64 (pred_aygtax, ay, ax);
+ float64x2_t z = vdivq_f64 (num, den);
+
+ /* Work out the correct shift for atan2:
+ Multiplication by pi is done later.
+ -pi when x < 0 and ax < ay
+ -pi/2 when x < 0 and ax > ay
+ 0 when x >= 0 and ax < ay
+ pi/2 when x >= 0 and ax > ay. */
+ float64x2_t shift = vreinterpretq_f64_u64 (
+ vandq_u64 (pred_xlt0, vreinterpretq_u64_f64 (v_f64 (-2.0))));
+ float64x2_t shift2 = vreinterpretq_f64_u64 (
+ vandq_u64 (pred_aygtax, vreinterpretq_u64_f64 (v_f64 (1.0))));
+ shift = vaddq_f64 (shift, shift2);
+
+ /* Calculate the polynomial approximation. */
float64x2_t z2 = vmulq_f64 (z, z);
- float64x2_t x2 = vmulq_f64 (z2, z2);
- float64x2_t x4 = vmulq_f64 (x2, x2);
- float64x2_t x8 = vmulq_f64 (x4, x4);
+ float64x2_t z3 = vmulq_f64 (z2, z);
+ float64x2_t z4 = vmulq_f64 (z2, z2);
+ float64x2_t z8 = vmulq_f64 (z4, z4);
+ float64x2_t z16 = vmulq_f64 (z8, z8);
float64x2_t c13 = vld1q_f64 (&d->c1);
float64x2_t c57 = vld1q_f64 (&d->c5);
float64x2_t c1315 = vld1q_f64 (&d->c13);
float64x2_t c1719 = vld1q_f64 (&d->c17);
- /* estrin_7. */
+ /* Order-7 Estrin. */
float64x2_t p01 = vfmaq_laneq_f64 (d->c0, z2, c13, 0);
float64x2_t p23 = vfmaq_laneq_f64 (d->c2, z2, c13, 1);
- float64x2_t p03 = vfmaq_f64 (p01, x2, p23);
+ float64x2_t p03 = vfmaq_f64 (p01, z4, p23);
float64x2_t p45 = vfmaq_laneq_f64 (d->c4, z2, c57, 0);
float64x2_t p67 = vfmaq_laneq_f64 (d->c6, z2, c57, 1);
- float64x2_t p47 = vfmaq_f64 (p45, x2, p67);
+ float64x2_t p47 = vfmaq_f64 (p45, z4, p67);
- float64x2_t p07 = vfmaq_f64 (p03, x4, p47);
+ float64x2_t p07 = vfmaq_f64 (p03, z8, p47);
- /* estrin_11. */
+ /* Order-11 Estrin. */
float64x2_t p89 = vfmaq_laneq_f64 (d->c8, z2, c911, 0);
float64x2_t p1011 = vfmaq_laneq_f64 (d->c10, z2, c911, 1);
- float64x2_t p811 = vfmaq_f64 (p89, x2, p1011);
+ float64x2_t p811 = vfmaq_f64 (p89, z4, p1011);
float64x2_t p1213 = vfmaq_laneq_f64 (d->c12, z2, c1315, 0);
float64x2_t p1415 = vfmaq_laneq_f64 (d->c14, z2, c1315, 1);
- float64x2_t p1215 = vfmaq_f64 (p1213, x2, p1415);
+ float64x2_t p1215 = vfmaq_f64 (p1213, z4, p1415);
float64x2_t p1617 = vfmaq_laneq_f64 (d->c16, z2, c1719, 0);
float64x2_t p1819 = vfmaq_laneq_f64 (d->c18, z2, c1719, 1);
- float64x2_t p1619 = vfmaq_f64 (p1617, x2, p1819);
+ float64x2_t p1619 = vfmaq_f64 (p1617, z4, p1819);
- float64x2_t p815 = vfmaq_f64 (p811, x4, p1215);
- float64x2_t p819 = vfmaq_f64 (p815, x8, p1619);
+ float64x2_t p815 = vfmaq_f64 (p811, z8, p1215);
+ float64x2_t p819 = vfmaq_f64 (p815, z16, p1619);
- float64x2_t ret = vfmaq_f64 (p07, p819, x8);
+ float64x2_t poly = vfmaq_f64 (p07, p819, z16);
/* Finalize. y = shift + z + z^3 * P(z^2). */
- ret = vfmaq_f64 (z, ret, vmulq_f64 (z2, z));
- ret = vaddq_f64 (ret, shift);
+ float64x2_t ret = vfmaq_f64 (z, shift, d->pi_over_2);
+ ret = vfmaq_f64 (ret, z3, poly);
if (__glibc_unlikely (v_any_u64 (special_cases)))
return special_case (y, x, ret, sign_xy, special_cases);
/* Account for the sign of x and y. */
- ret = vreinterpretq_f64_u64 (
+ return vreinterpretq_f64_u64 (
veorq_u64 (vreinterpretq_u64_f64 (ret), sign_xy));
-
- return ret;
}
#include "math_config.h"
#include "sv_math.h"
-#include "poly_sve_f64.h"
static const struct data
{
- float64_t poly[20];
- float64_t pi_over_2;
+ float64_t c0, c2, c4, c6, c8, c10, c12, c14, c16, c18;
+ float64_t c1, c3, c5, c7, c9, c11, c13, c15, c17, c19;
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-1022, 1.0]. */
- .poly = { -0x1.5555555555555p-2, 0x1.99999999996c1p-3, -0x1.2492492478f88p-3,
- 0x1.c71c71bc3951cp-4, -0x1.745d160a7e368p-4, 0x1.3b139b6a88ba1p-4,
- -0x1.11100ee084227p-4, 0x1.e1d0f9696f63bp-5, -0x1.aebfe7b418581p-5,
- 0x1.842dbe9b0d916p-5, -0x1.5d30140ae5e99p-5, 0x1.338e31eb2fbbcp-5,
- -0x1.00e6eece7de8p-5, 0x1.860897b29e5efp-6, -0x1.0051381722a59p-6,
- 0x1.14e9dc19a4a4ep-7, -0x1.d0062b42fe3bfp-9, 0x1.17739e210171ap-10,
- -0x1.ab24da7be7402p-13, 0x1.358851160a528p-16, },
- .pi_over_2 = 0x1.921fb54442d18p+0,
+ .c0 = -0x1.555555555552ap-2, .c1 = 0x1.9999999995aebp-3,
+ .c2 = -0x1.24924923923f6p-3, .c3 = 0x1.c71c7184288a2p-4,
+ .c4 = -0x1.745d11fb3d32bp-4, .c5 = 0x1.3b136a18051b9p-4,
+ .c6 = -0x1.110e6d985f496p-4, .c7 = 0x1.e1bcf7f08801dp-5,
+ .c8 = -0x1.ae644e28058c3p-5, .c9 = 0x1.82eeb1fed85c6p-5,
+ .c10 = -0x1.59d7f901566cbp-5, .c11 = 0x1.2c982855ab069p-5,
+ .c12 = -0x1.eb49592998177p-6, .c13 = 0x1.69d8b396e3d38p-6,
+ .c14 = -0x1.ca980345c4204p-7, .c15 = 0x1.dc050eafde0b3p-8,
+ .c16 = -0x1.7ea70755b8eccp-9, .c17 = 0x1.ba3da3de903e8p-11,
+ .c18 = -0x1.44a4b059b6f67p-13, .c19 = 0x1.c4a45029e5a91p-17,
};
-
/* Special cases i.e. 0, infinity, nan (fall back to scalar calls). */
static svfloat64_t NOINLINE
special_case (svfloat64_t y, svfloat64_t x, svfloat64_t ret,
}
/* Fast implementation of SVE atan2. Errors are greatest when y and
- x are reasonably close together. The greatest observed error is 2.28 ULP:
- _ZGVsMxvv_atan2 (-0x1.5915b1498e82fp+732, 0x1.54d11ef838826p+732)
- got -0x1.954f42f1fa841p-1 want -0x1.954f42f1fa843p-1. */
-svfloat64_t SV_NAME_D2 (atan2) (svfloat64_t y, svfloat64_t x, const svbool_t pg)
+ x are reasonably close together. The greatest observed error is 1.94 ULP:
+ _ZGVsMxvv_atan2 (0x1.8a4bf7167228ap+5, 0x1.84971226bb57bp+5)
+ got 0x1.95db19dfef9ccp-1 want 0x1.95db19dfef9cep-1. */
+svfloat64_t SV_NAME_D2 (atan2) (svfloat64_t y, svfloat64_t x,
+ const svbool_t pg)
{
- const struct data *data_ptr = ptr_barrier (&data);
+ const struct data *d = ptr_barrier (&data);
svuint64_t ix = svreinterpret_u64 (x);
svuint64_t iy = svreinterpret_u64 (y);
+ svbool_t ptrue = svptrue_b64 ();
svbool_t cmp_x = zeroinfnan (ix, pg);
svbool_t cmp_y = zeroinfnan (iy, pg);
svbool_t pred_aygtax = svcmpgt (pg, ay, ax);
- /* Set up z for call to atan. */
- svfloat64_t n = svsel (pred_aygtax, svneg_x (pg, ax), ay);
- svfloat64_t d = svsel (pred_aygtax, ay, ax);
- svfloat64_t z = svdiv_x (pg, n, d);
-
- /* Work out the correct shift. */
+ /* Set up z for evaluation of atan. */
+ svfloat64_t num = svsel (pred_aygtax, svneg_x (pg, ax), ay);
+ svfloat64_t den = svsel (pred_aygtax, ay, ax);
+ svfloat64_t z = svdiv_x (pg, num, den);
+
+ /* Work out the correct shift for atan2:
+ Multiplication by pi is done later.
+ -pi when x < 0 and ax < ay
+ -pi/2 when x < 0 and ax > ay
+ 0 when x >= 0 and ax < ay
+ pi/2 when x >= 0 and ax > ay. */
svfloat64_t shift = svreinterpret_f64 (svlsr_x (pg, sign_x, 1));
+ svfloat64_t shift_mul = svreinterpret_f64 (
+ svorr_x (pg, sign_x, svreinterpret_u64 (sv_f64 (0x1.921fb54442d18p+0))));
shift = svsel (pred_aygtax, sv_f64 (1.0), shift);
- shift = svreinterpret_f64 (svorr_x (pg, sign_x, svreinterpret_u64 (shift)));
- shift = svmul_x (pg, shift, data_ptr->pi_over_2);
+ shift = svmla_x (pg, z, shift, shift_mul);
/* Use split Estrin scheme for P(z^2) with deg(P)=19. */
svfloat64_t z2 = svmul_x (pg, z, z);
- svfloat64_t x2 = svmul_x (pg, z2, z2);
- svfloat64_t x4 = svmul_x (pg, x2, x2);
- svfloat64_t x8 = svmul_x (pg, x4, x4);
+ svfloat64_t z3 = svmul_x (pg, z2, z);
+ svfloat64_t z4 = svmul_x (pg, z2, z2);
+ svfloat64_t z8 = svmul_x (pg, z4, z4);
+ svfloat64_t z16 = svmul_x (pg, z8, z8);
- svfloat64_t ret = svmla_x (
- pg, sv_estrin_7_f64_x (pg, z2, x2, x4, data_ptr->poly),
- sv_estrin_11_f64_x (pg, z2, x2, x4, x8, data_ptr->poly + 8), x8);
+ /* Order-7 Estrin. */
+ svfloat64_t c13 = svld1rq (ptrue, &d->c1);
+ svfloat64_t c57 = svld1rq (ptrue, &d->c5);
- /* y = shift + z + z^3 * P(z^2). */
- svfloat64_t z3 = svmul_x (pg, z2, z);
- ret = svmla_x (pg, z, z3, ret);
+ svfloat64_t p01 = svmla_lane (sv_f64 (d->c0), z2, c13, 0);
+ svfloat64_t p23 = svmla_lane (sv_f64 (d->c2), z2, c13, 1);
+ svfloat64_t p45 = svmla_lane (sv_f64 (d->c4), z2, c57, 0);
+ svfloat64_t p67 = svmla_lane (sv_f64 (d->c6), z2, c57, 1);
+
+ svfloat64_t p03 = svmla_x (pg, p01, z4, p23);
+ svfloat64_t p47 = svmla_x (pg, p45, z4, p67);
+ svfloat64_t p07 = svmla_x (pg, p03, z8, p47);
+
+ /* Order-11 Estrin. */
+ svfloat64_t c911 = svld1rq (ptrue, &d->c9);
+ svfloat64_t c1315 = svld1rq (ptrue, &d->c13);
+ svfloat64_t c1719 = svld1rq (ptrue, &d->c17);
- ret = svadd_m (pg, ret, shift);
+ svfloat64_t p89 = svmla_lane (sv_f64 (d->c8), z2, c911, 0);
+ svfloat64_t p1011 = svmla_lane (sv_f64 (d->c10), z2, c911, 1);
+ svfloat64_t p811 = svmla_x (pg, p89, z4, p1011);
+
+ svfloat64_t p1213 = svmla_lane (sv_f64 (d->c12), z2, c1315, 0);
+ svfloat64_t p1415 = svmla_lane (sv_f64 (d->c14), z2, c1315, 1);
+ svfloat64_t p1215 = svmla_x (pg, p1213, z4, p1415);
+
+ svfloat64_t p1617 = svmla_lane (sv_f64 (d->c16), z2, c1719, 0);
+ svfloat64_t p1819 = svmla_lane (sv_f64 (d->c18), z2, c1719, 1);
+ svfloat64_t p1619 = svmla_x (pg, p1617, z4, p1819);
+
+ svfloat64_t p815 = svmla_x (pg, p811, z8, p1215);
+ svfloat64_t p819 = svmla_x (pg, p815, z16, p1619);
+
+ svfloat64_t poly = svmla_x (pg, p07, z16, p819);
+
+ /* y = shift + z + z^3 * P(z^2). */
+ svfloat64_t ret = svmla_x (pg, shift, z3, poly);
/* Account for the sign of x and y. */
if (__glibc_unlikely (svptest_any (pg, cmp_xy)))
<https://www.gnu.org/licenses/>. */
#include "v_math.h"
-#include "poly_advsimd_f32.h"
static const struct data
{
- float32x4_t c0, pi_over_2, c4, c6, c2;
+ float32x4_t c0, c4, c6, c2;
float c1, c3, c5, c7;
uint32x4_t comp_const;
+ float32x4_t pi;
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-128, 1.0].
Generated using fpminimax between FLT_MIN and 1. */
- .c0 = V4 (-0x1.55555p-2f), .c1 = 0x1.99935ep-3f,
- .c2 = V4 (-0x1.24051ep-3f), .c3 = 0x1.bd7368p-4f,
- .c4 = V4 (-0x1.491f0ep-4f), .c5 = 0x1.93a2c0p-5f,
- .c6 = V4 (-0x1.4c3c60p-6f), .c7 = 0x1.01fd88p-8f,
- .pi_over_2 = V4 (0x1.921fb6p+0f), .comp_const = V4 (2 * 0x7f800000lu - 1),
+ .c0 = V4 (-0x1.5554dcp-2), .c1 = 0x1.9978ecp-3,
+ .c2 = V4 (-0x1.230a94p-3), .c3 = 0x1.b4debp-4,
+ .c4 = V4 (-0x1.3550dap-4), .c5 = 0x1.61eebp-5,
+ .c6 = V4 (-0x1.0c17d4p-6), .c7 = 0x1.7ea694p-9,
+ .pi = V4 (0x1.921fb6p+1f), .comp_const = V4 (2 * 0x7f800000lu - 1),
};
#define SignMask v_u32 (0x80000000)
zeroinfnan (uint32x4_t i, const struct data *d)
{
/* 2 * i - 1 >= 2 * 0x7f800000lu - 1. */
- return vcgeq_u32 (vsubq_u32 (vmulq_n_u32 (i, 2), v_u32 (1)), d->comp_const);
+ return vcgeq_u32 (vsubq_u32 (vshlq_n_u32 (i, 1), v_u32 (1)), d->comp_const);
}
/* Fast implementation of vector atan2f. Maximum observed error is
- 2.95 ULP in [0x1.9300d6p+6 0x1.93c0c6p+6] x [0x1.8c2dbp+6 0x1.8cea6p+6]:
- _ZGVnN4vv_atan2f (0x1.93836cp+6, 0x1.8cae1p+6) got 0x1.967f06p-1
- want 0x1.967f00p-1. */
+ 2.13 ULP in [0x1.9300d6p+6 0x1.93c0c6p+6] x [0x1.8c2dbp+6 0x1.8cea6p+6]:
+ _ZGVnN4vv_atan2f (0x1.14a9d4p-87, 0x1.0eb886p-87) got 0x1.97aea2p-1
+ want 0x1.97ae9ep-1. */
float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (atan2) (float32x4_t y, float32x4_t x)
{
const struct data *d = ptr_barrier (&data);
uint32x4_t pred_xlt0 = vcltzq_f32 (x);
uint32x4_t pred_aygtax = vcgtq_f32 (ay, ax);
- /* Set up z for call to atanf. */
- float32x4_t n = vbslq_f32 (pred_aygtax, vnegq_f32 (ax), ay);
- float32x4_t q = vbslq_f32 (pred_aygtax, ay, ax);
- float32x4_t z = vdivq_f32 (n, q);
-
- /* Work out the correct shift. */
+ /* Set up z for evaluation of atanf. */
+ float32x4_t num = vbslq_f32 (pred_aygtax, vnegq_f32 (ax), ay);
+ float32x4_t den = vbslq_f32 (pred_aygtax, ay, ax);
+ float32x4_t z = vdivq_f32 (num, den);
+
+ /* Work out the correct shift for atan2:
+ Multiplication by pi is done later.
+ -pi when x < 0 and ax < ay
+ -pi/2 when x < 0 and ax > ay
+ 0 when x >= 0 and ax < ay
+ pi/2 when x >= 0 and ax > ay. */
float32x4_t shift = vreinterpretq_f32_u32 (
- vandq_u32 (pred_xlt0, vreinterpretq_u32_f32 (v_f32 (-2.0f))));
- shift = vbslq_f32 (pred_aygtax, vaddq_f32 (shift, v_f32 (1.0f)), shift);
- shift = vmulq_f32 (shift, d->pi_over_2);
-
- /* Calculate the polynomial approximation.
- Use 2-level Estrin scheme for P(z^2) with deg(P)=7. However,
- a standard implementation using z8 creates spurious underflow
- in the very last fma (when z^8 is small enough).
- Therefore, we split the last fma into a mul and an fma.
- Horner and single-level Estrin have higher errors that exceed
- threshold. */
+ vandq_u32 (pred_xlt0, vreinterpretq_u32_f32 (v_f32 (-1.0f))));
+ float32x4_t shift2 = vreinterpretq_f32_u32 (
+ vandq_u32 (pred_aygtax, vreinterpretq_u32_f32 (v_f32 (0.5f))));
+ shift = vaddq_f32 (shift, shift2);
+
+ /* Calculate the polynomial approximation. */
float32x4_t z2 = vmulq_f32 (z, z);
+ float32x4_t z3 = vmulq_f32 (z2, z);
float32x4_t z4 = vmulq_f32 (z2, z2);
+ float32x4_t z8 = vmulq_f32 (z4, z4);
float32x4_t c1357 = vld1q_f32 (&d->c1);
+
float32x4_t p01 = vfmaq_laneq_f32 (d->c0, z2, c1357, 0);
float32x4_t p23 = vfmaq_laneq_f32 (d->c2, z2, c1357, 1);
float32x4_t p45 = vfmaq_laneq_f32 (d->c4, z2, c1357, 2);
float32x4_t p03 = vfmaq_f32 (p01, z4, p23);
float32x4_t p47 = vfmaq_f32 (p45, z4, p67);
- float32x4_t ret = vfmaq_f32 (p03, z4, vmulq_f32 (z4, p47));
+ float32x4_t poly = vfmaq_f32 (p03, z8, p47);
/* y = shift + z * P(z^2). */
- ret = vaddq_f32 (vfmaq_f32 (z, ret, vmulq_f32 (z2, z)), shift);
+ float32x4_t ret = vfmaq_f32 (z, shift, d->pi);
+ ret = vfmaq_f32 (ret, z3, poly);
if (__glibc_unlikely (v_any_u32 (special_cases)))
{
<https://www.gnu.org/licenses/>. */
#include "sv_math.h"
-#include "poly_sve_f32.h"
static const struct data
{
- float32_t poly[8];
+ float32_t c0, c2, c4, c6;
+ float32_t c1, c3, c5, c7;
float32_t pi_over_2;
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-128, 1.0]. */
- .poly = { -0x1.55555p-2f, 0x1.99935ep-3f, -0x1.24051ep-3f, 0x1.bd7368p-4f,
- -0x1.491f0ep-4f, 0x1.93a2c0p-5f, -0x1.4c3c60p-6f, 0x1.01fd88p-8f },
- .pi_over_2 = 0x1.921fb6p+0f,
+ .c0 = -0x1.5554dcp-2, .c1 = 0x1.9978ecp-3, .c2 = -0x1.230a94p-3,
+ .c3 = 0x1.b4debp-4, .c4 = -0x1.3550dap-4, .c5 = 0x1.61eebp-5,
+ .c6 = -0x1.0c17d4p-6, .c7 = 0x1.7ea694p-9, .pi_over_2 = 0x1.921fb6p+0f,
};
/* Special cases i.e. 0, infinity, nan (fall back to scalar calls). */
/* Fast implementation of SVE atan2f based on atan(x) ~ shift + z + z^3 *
P(z^2) with reduction to [0,1] using z=1/x and shift = pi/2. Maximum
- observed error is 2.95 ULP:
- _ZGVsMxvv_atan2f (0x1.93836cp+6, 0x1.8cae1p+6) got 0x1.967f06p-1
- want 0x1.967f00p-1. */
-svfloat32_t SV_NAME_F2 (atan2) (svfloat32_t y, svfloat32_t x, const svbool_t pg)
+ observed error is 2.21 ULP:
+ _ZGVnN4vv_atan2f (0x1.a04aa8p+6, 0x1.9a274p+6) got 0x1.95ed3ap-1
+ want 0x1.95ed36p-1. */
+svfloat32_t SV_NAME_F2 (atan2) (svfloat32_t y, svfloat32_t x,
+ const svbool_t pg)
{
- const struct data *data_ptr = ptr_barrier (&data);
+ const struct data *d = ptr_barrier (&data);
+ svbool_t ptrue = svptrue_b32 ();
svuint32_t ix = svreinterpret_u32 (x);
svuint32_t iy = svreinterpret_u32 (y);
svbool_t pred_aygtax = svcmpgt (pg, ay, ax);
- /* Set up z for call to atan. */
- svfloat32_t n = svsel (pred_aygtax, svneg_x (pg, ax), ay);
- svfloat32_t d = svsel (pred_aygtax, ay, ax);
- svfloat32_t z = svdiv_x (pg, n, d);
-
- /* Work out the correct shift. */
+ /* Set up z for evaluation of atanf. */
+ svfloat32_t num = svsel (pred_aygtax, svneg_x (pg, ax), ay);
+ svfloat32_t den = svsel (pred_aygtax, ay, ax);
+ svfloat32_t z = svdiv_x (ptrue, num, den);
+
+ /* Work out the correct shift for atan2:
+ Multiplication by pi is done later.
+ -pi when x < 0 and ax < ay
+ -pi/2 when x < 0 and ax > ay
+ 0 when x >= 0 and ax < ay
+ pi/2 when x >= 0 and ax > ay. */
svfloat32_t shift = svreinterpret_f32 (svlsr_x (pg, sign_x, 1));
shift = svsel (pred_aygtax, sv_f32 (1.0), shift);
shift = svreinterpret_f32 (svorr_x (pg, sign_x, svreinterpret_u32 (shift)));
- shift = svmul_x (pg, shift, sv_f32 (data_ptr->pi_over_2));
/* Use pure Estrin scheme for P(z^2) with deg(P)=7. */
- svfloat32_t z2 = svmul_x (pg, z, z);
+ svfloat32_t z2 = svmul_x (ptrue, z, z);
+ svfloat32_t z3 = svmul_x (pg, z2, z);
svfloat32_t z4 = svmul_x (pg, z2, z2);
svfloat32_t z8 = svmul_x (pg, z4, z4);
- svfloat32_t ret = sv_estrin_7_f32_x (pg, z2, z4, z8, data_ptr->poly);
+ svfloat32_t odd_coeffs = svld1rq (ptrue, &d->c1);
- /* ret = shift + z + z^3 * P(z^2). */
- svfloat32_t z3 = svmul_x (pg, z2, z);
- ret = svmla_x (pg, z, z3, ret);
+ svfloat32_t p01 = svmla_lane (sv_f32 (d->c0), z2, odd_coeffs, 0);
+ svfloat32_t p23 = svmla_lane (sv_f32 (d->c2), z2, odd_coeffs, 1);
+ svfloat32_t p45 = svmla_lane (sv_f32 (d->c4), z2, odd_coeffs, 2);
+ svfloat32_t p67 = svmla_lane (sv_f32 (d->c6), z2, odd_coeffs, 3);
- ret = svadd_m (pg, ret, shift);
+ svfloat32_t p03 = svmla_x (pg, p01, z4, p23);
+ svfloat32_t p47 = svmla_x (pg, p45, z4, p67);
+
+ svfloat32_t poly = svmla_x (pg, p03, z8, p47);
+
+ /* ret = shift + z + z^3 * P(z^2). */
+ svfloat32_t ret = svmla_x (pg, z, shift, sv_f32 (d->pi_over_2));
+ ret = svmla_x (pg, ret, z3, poly);
/* Account for the sign of x and y. */
<https://www.gnu.org/licenses/>. */
#include "v_math.h"
-#include "poly_advsimd_f64.h"
static const struct data
{
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-1022, 1.0]. */
- .c0 = V2 (-0x1.5555555555555p-2), .c1 = 0x1.99999999996c1p-3,
- .c2 = V2 (-0x1.2492492478f88p-3), .c3 = 0x1.c71c71bc3951cp-4,
- .c4 = V2 (-0x1.745d160a7e368p-4), .c5 = 0x1.3b139b6a88ba1p-4,
- .c6 = V2 (-0x1.11100ee084227p-4), .c7 = 0x1.e1d0f9696f63bp-5,
- .c8 = V2 (-0x1.aebfe7b418581p-5), .c9 = 0x1.842dbe9b0d916p-5,
- .c10 = V2 (-0x1.5d30140ae5e99p-5), .c11 = 0x1.338e31eb2fbbcp-5,
- .c12 = V2 (-0x1.00e6eece7de8p-5), .c13 = 0x1.860897b29e5efp-6,
- .c14 = V2 (-0x1.0051381722a59p-6), .c15 = 0x1.14e9dc19a4a4ep-7,
- .c16 = V2 (-0x1.d0062b42fe3bfp-9), .c17 = 0x1.17739e210171ap-10,
- .c18 = V2 (-0x1.ab24da7be7402p-13), .c19 = 0x1.358851160a528p-16,
+ .c0 = V2 (-0x1.555555555552ap-2), .c1 = 0x1.9999999995aebp-3,
+ .c2 = V2 (-0x1.24924923923f6p-3), .c3 = 0x1.c71c7184288a2p-4,
+ .c4 = V2 (-0x1.745d11fb3d32bp-4), .c5 = 0x1.3b136a18051b9p-4,
+ .c6 = V2 (-0x1.110e6d985f496p-4), .c7 = 0x1.e1bcf7f08801dp-5,
+ .c8 = V2 (-0x1.ae644e28058c3p-5), .c9 = 0x1.82eeb1fed85c6p-5,
+ .c10 = V2 (-0x1.59d7f901566cbp-5), .c11 = 0x1.2c982855ab069p-5,
+ .c12 = V2 (-0x1.eb49592998177p-6), .c13 = 0x1.69d8b396e3d38p-6,
+ .c14 = V2 (-0x1.ca980345c4204p-7), .c15 = 0x1.dc050eafde0b3p-8,
+ .c16 = V2 (-0x1.7ea70755b8eccp-9), .c17 = 0x1.ba3da3de903e8p-11,
+ .c18 = V2 (-0x1.44a4b059b6f67p-13), .c19 = 0x1.c4a45029e5a91p-17,
.pi_over_2 = V2 (0x1.921fb54442d18p+0),
};
/* Fast implementation of vector atan.
Based on atan(x) ~ shift + z + z^3 * P(z^2) with reduction to [0,1] using
- z=1/x and shift = pi/2. Maximum observed error is 2.27 ulps:
- _ZGVnN2v_atan (0x1.0005af27c23e9p+0) got 0x1.9225645bdd7c1p-1
- want 0x1.9225645bdd7c3p-1. */
+ z=1/x and shift = pi/2. Maximum observed error is 2.45 ulps:
+ _ZGVnN2v_atan (0x1.0008d737eb3e6p+0) got 0x1.92288c551a4c1p-1
+ want 0x1.92288c551a4c3p-1. */
float64x2_t VPCS_ATTR V_NAME_D1 (atan) (float64x2_t x)
{
const struct data *d = ptr_barrier (&data);
y := arctan(x) for x < 1
y := pi/2 + arctan(-1/x) for x > 1
Hence, use z=-1/a if x>=1, otherwise z=a. */
- uint64x2_t red = vcagtq_f64 (x, v_f64 (1.0));
+ uint64x2_t red = vcagtq_f64 (x, v_f64 (-1.0));
/* Avoid dependency in abs(x) in division (and comparison). */
- float64x2_t z = vbslq_f64 (red, vdivq_f64 (v_f64 (1.0), x), x);
+ float64x2_t z = vbslq_f64 (red, vdivq_f64 (v_f64 (-1.0), x), x);
+
float64x2_t shift = vreinterpretq_f64_u64 (
vandq_u64 (red, vreinterpretq_u64_f64 (d->pi_over_2)));
- /* Use absolute value only when needed (odd powers of z). */
- float64x2_t az = vbslq_f64 (
- SignMask, vreinterpretq_f64_u64 (vandq_u64 (SignMask, red)), z);
-
- /* Calculate the polynomial approximation.
- Use split Estrin scheme for P(z^2) with deg(P)=19. Use split instead of
- full scheme to avoid underflow in x^16.
- The order 19 polynomial P approximates
- (atan(sqrt(x))-sqrt(x))/x^(3/2). */
+
+ /* Reinsert sign bit from argument into the shift value. */
+ shift = vreinterpretq_f64_u64 (
+ veorq_u64 (vreinterpretq_u64_f64 (shift), sign));
+
+ /* Calculate polynomial approximation P(z^2) with deg(P)=19. */
float64x2_t z2 = vmulq_f64 (z, z);
- float64x2_t x2 = vmulq_f64 (z2, z2);
- float64x2_t x4 = vmulq_f64 (x2, x2);
- float64x2_t x8 = vmulq_f64 (x4, x4);
+ float64x2_t z4 = vmulq_f64 (z2, z2);
+ float64x2_t z8 = vmulq_f64 (z4, z4);
+ float64x2_t z16 = vmulq_f64 (z8, z8);
- /* estrin_7. */
+ /* Order-7 Estrin. */
float64x2_t p01 = vfmaq_laneq_f64 (d->c0, z2, c13, 0);
float64x2_t p23 = vfmaq_laneq_f64 (d->c2, z2, c13, 1);
- float64x2_t p03 = vfmaq_f64 (p01, x2, p23);
+ float64x2_t p03 = vfmaq_f64 (p01, z4, p23);
float64x2_t p45 = vfmaq_laneq_f64 (d->c4, z2, c57, 0);
float64x2_t p67 = vfmaq_laneq_f64 (d->c6, z2, c57, 1);
- float64x2_t p47 = vfmaq_f64 (p45, x2, p67);
+ float64x2_t p47 = vfmaq_f64 (p45, z4, p67);
- float64x2_t p07 = vfmaq_f64 (p03, x4, p47);
+ float64x2_t p07 = vfmaq_f64 (p03, z8, p47);
- /* estrin_11. */
+ /* Order-11 Estrin. */
float64x2_t p89 = vfmaq_laneq_f64 (d->c8, z2, c911, 0);
float64x2_t p1011 = vfmaq_laneq_f64 (d->c10, z2, c911, 1);
- float64x2_t p811 = vfmaq_f64 (p89, x2, p1011);
+ float64x2_t p811 = vfmaq_f64 (p89, z4, p1011);
float64x2_t p1213 = vfmaq_laneq_f64 (d->c12, z2, c1315, 0);
float64x2_t p1415 = vfmaq_laneq_f64 (d->c14, z2, c1315, 1);
- float64x2_t p1215 = vfmaq_f64 (p1213, x2, p1415);
+ float64x2_t p1215 = vfmaq_f64 (p1213, z4, p1415);
float64x2_t p1617 = vfmaq_laneq_f64 (d->c16, z2, c1719, 0);
float64x2_t p1819 = vfmaq_laneq_f64 (d->c18, z2, c1719, 1);
- float64x2_t p1619 = vfmaq_f64 (p1617, x2, p1819);
+ float64x2_t p1619 = vfmaq_f64 (p1617, z4, p1819);
- float64x2_t p815 = vfmaq_f64 (p811, x4, p1215);
- float64x2_t p819 = vfmaq_f64 (p815, x8, p1619);
+ float64x2_t p815 = vfmaq_f64 (p811, z8, p1215);
+ float64x2_t p819 = vfmaq_f64 (p815, z16, p1619);
- float64x2_t y = vfmaq_f64 (p07, p819, x8);
+ float64x2_t y = vfmaq_f64 (p07, p819, z16);
/* Finalize. y = shift + z + z^3 * P(z^2). */
- y = vfmaq_f64 (az, y, vmulq_f64 (z2, az));
- y = vaddq_f64 (y, shift);
-
- /* y = atan(x) if x>0, -atan(-x) otherwise. */
- y = vreinterpretq_f64_u64 (veorq_u64 (vreinterpretq_u64_f64 (y), sign));
- return y;
+ y = vfmsq_f64 (v_f64 (-1.0), z2, y);
+ return vfmsq_f64 (shift, z, y);
}
<https://www.gnu.org/licenses/>. */
#include "sv_math.h"
-#include "poly_sve_f64.h"
static const struct data
{
- float64_t poly[20];
- float64_t pi_over_2;
+ float64_t c0, c2, c4, c6, c8, c10, c12, c14, c16, c18;
+ float64_t c1, c3, c5, c7, c9, c11, c13, c15, c17, c19;
+ float64_t shift_val, neg_one;
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-1022, 1.0]. */
- .poly = { -0x1.5555555555555p-2, 0x1.99999999996c1p-3, -0x1.2492492478f88p-3,
- 0x1.c71c71bc3951cp-4, -0x1.745d160a7e368p-4, 0x1.3b139b6a88ba1p-4,
- -0x1.11100ee084227p-4, 0x1.e1d0f9696f63bp-5, -0x1.aebfe7b418581p-5,
- 0x1.842dbe9b0d916p-5, -0x1.5d30140ae5e99p-5, 0x1.338e31eb2fbbcp-5,
- -0x1.00e6eece7de8p-5, 0x1.860897b29e5efp-6, -0x1.0051381722a59p-6,
- 0x1.14e9dc19a4a4ep-7, -0x1.d0062b42fe3bfp-9, 0x1.17739e210171ap-10,
- -0x1.ab24da7be7402p-13, 0x1.358851160a528p-16, },
- .pi_over_2 = 0x1.921fb54442d18p+0,
+ .c0 = -0x1.555555555552ap-2, .c1 = 0x1.9999999995aebp-3,
+ .c2 = -0x1.24924923923f6p-3, .c3 = 0x1.c71c7184288a2p-4,
+ .c4 = -0x1.745d11fb3d32bp-4, .c5 = 0x1.3b136a18051b9p-4,
+ .c6 = -0x1.110e6d985f496p-4, .c7 = 0x1.e1bcf7f08801dp-5,
+ .c8 = -0x1.ae644e28058c3p-5, .c9 = 0x1.82eeb1fed85c6p-5,
+ .c10 = -0x1.59d7f901566cbp-5, .c11 = 0x1.2c982855ab069p-5,
+ .c12 = -0x1.eb49592998177p-6, .c13 = 0x1.69d8b396e3d38p-6,
+ .c14 = -0x1.ca980345c4204p-7, .c15 = 0x1.dc050eafde0b3p-8,
+ .c16 = -0x1.7ea70755b8eccp-9, .c17 = 0x1.ba3da3de903e8p-11,
+ .c18 = -0x1.44a4b059b6f67p-13, .c19 = 0x1.c4a45029e5a91p-17,
+ .shift_val = 0x1.490fdaa22168cp+1, .neg_one = -1,
};
/* Useful constants. */
/* Fast implementation of SVE atan.
Based on atan(x) ~ shift + z + z^3 * P(z^2) with reduction to [0,1] using
z=1/x and shift = pi/2. Largest errors are close to 1. The maximum observed
- error is 2.27 ulps:
- _ZGVsMxv_atan (0x1.0005af27c23e9p+0) got 0x1.9225645bdd7c1p-1
- want 0x1.9225645bdd7c3p-1. */
+ error is 2.08 ulps:
+ _ZGVsMxv_atan (0x1.000a7c56975e8p+0) got 0x1.922a3163e15c2p-1
+ want 0x1.922a3163e15c4p-1. */
svfloat64_t SV_NAME_D1 (atan) (svfloat64_t x, const svbool_t pg)
{
const struct data *d = ptr_barrier (&data);
- /* No need to trigger special case. Small cases, infs and nans
- are supported by our approximation technique. */
+ svbool_t ptrue = svptrue_b64 ();
svuint64_t ix = svreinterpret_u64 (x);
svuint64_t sign = svand_x (pg, ix, SignMask);
y := arctan(x) for x < 1
y := pi/2 + arctan(-1/x) for x > 1
Hence, use z=-1/a if x>=1, otherwise z=a. */
- svbool_t red = svacgt (pg, x, 1.0);
- /* Avoid dependency in abs(x) in division (and comparison). */
- svfloat64_t z = svsel (red, svdivr_x (pg, x, 1.0), x);
- /* Use absolute value only when needed (odd powers of z). */
- svfloat64_t az = svabs_x (pg, z);
- az = svneg_m (az, red, az);
+ svbool_t red = svacgt (pg, x, d->neg_one);
+ svfloat64_t z = svsel (red, svdiv_x (pg, sv_f64 (d->neg_one), x), x);
+
+ /* Reuse of -1.0f to reduce constant loads,
+ We need a shift value of 1/2, which is created via -1 + (1 + 1/2). */
+ svfloat64_t shift
+ = svadd_z (red, sv_f64 (d->neg_one), sv_f64 (d->shift_val));
+
+ /* Reinserts the sign bit of the argument to handle the case of x < -1. */
+ shift = svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (shift), sign));
/* Use split Estrin scheme for P(z^2) with deg(P)=19. */
- svfloat64_t z2 = svmul_x (pg, z, z);
- svfloat64_t x2 = svmul_x (pg, z2, z2);
- svfloat64_t x4 = svmul_x (pg, x2, x2);
- svfloat64_t x8 = svmul_x (pg, x4, x4);
+ svfloat64_t z2 = svmul_x (ptrue, z, z);
+ svfloat64_t z4 = svmul_x (ptrue, z2, z2);
+ svfloat64_t z8 = svmul_x (ptrue, z4, z4);
+ svfloat64_t z16 = svmul_x (ptrue, z8, z8);
- svfloat64_t y
- = svmla_x (pg, sv_estrin_7_f64_x (pg, z2, x2, x4, d->poly),
- sv_estrin_11_f64_x (pg, z2, x2, x4, x8, d->poly + 8), x8);
+ /* Order-7 Estrin. */
+ svfloat64_t c13 = svld1rq (ptrue, &d->c1);
+ svfloat64_t c57 = svld1rq (ptrue, &d->c5);
- /* y = shift + z + z^3 * P(z^2). */
- svfloat64_t z3 = svmul_x (pg, z2, az);
- y = svmla_x (pg, az, z3, y);
+ svfloat64_t p01 = svmla_lane (sv_f64 (d->c0), z2, c13, 0);
+ svfloat64_t p23 = svmla_lane (sv_f64 (d->c2), z2, c13, 1);
+ svfloat64_t p45 = svmla_lane (sv_f64 (d->c4), z2, c57, 0);
+ svfloat64_t p67 = svmla_lane (sv_f64 (d->c6), z2, c57, 1);
+
+ svfloat64_t p03 = svmla_x (pg, p01, z4, p23);
+ svfloat64_t p47 = svmla_x (pg, p45, z4, p67);
+ svfloat64_t p07 = svmla_x (pg, p03, z8, p47);
+
+ /* Order-11 Estrin. */
+ svfloat64_t c911 = svld1rq (ptrue, &d->c9);
+ svfloat64_t c1315 = svld1rq (ptrue, &d->c13);
+ svfloat64_t c1719 = svld1rq (ptrue, &d->c17);
- /* Apply shift as indicated by `red` predicate. */
- y = svadd_m (red, y, d->pi_over_2);
+ svfloat64_t p89 = svmla_lane (sv_f64 (d->c8), z2, c911, 0);
+ svfloat64_t p1011 = svmla_lane (sv_f64 (d->c10), z2, c911, 1);
+ svfloat64_t p811 = svmla_x (pg, p89, z4, p1011);
- /* y = atan(x) if x>0, -atan(-x) otherwise. */
- y = svreinterpret_f64 (sveor_x (pg, svreinterpret_u64 (y), sign));
+ svfloat64_t p1213 = svmla_lane (sv_f64 (d->c12), z2, c1315, 0);
+ svfloat64_t p1415 = svmla_lane (sv_f64 (d->c14), z2, c1315, 1);
+ svfloat64_t p1215 = svmla_x (pg, p1213, z4, p1415);
- return y;
+ svfloat64_t p1617 = svmla_lane (sv_f64 (d->c16), z2, c1719, 0);
+ svfloat64_t p1819 = svmla_lane (sv_f64 (d->c18), z2, c1719, 1);
+ svfloat64_t p1619 = svmla_x (pg, p1617, z4, p1819);
+
+ svfloat64_t p815 = svmla_x (pg, p811, z8, p1215);
+ svfloat64_t p819 = svmla_x (pg, p815, z16, p1619);
+
+ svfloat64_t y = svmla_x (pg, p07, z16, p819);
+
+ /* y = shift + z + z^3 * P(z^2). */
+ shift = svadd_m (red, z, shift);
+ y = svmul_x (pg, z2, y);
+ return svmla_x (pg, shift, z, y);
}
static const struct data
{
+ uint32x4_t sign_mask, pi_over_2;
+ float32x4_t neg_one;
+#if WANT_SIMD_EXCEPT
float32x4_t poly[8];
- float32x4_t pi_over_2;
+} data = {
+ .poly = { V4 (-0x1.5554dcp-2), V4 (0x1.9978ecp-3), V4 (-0x1.230a94p-3),
+ V4 (0x1.b4debp-4), V4 (-0x1.3550dap-4), V4 (0x1.61eebp-5),
+ V4 (-0x1.0c17d4p-6), V4 (0x1.7ea694p-9) },
+#else
+ float32x4_t c0, c2, c4, c6;
+ float c1, c3, c5, c7;
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-128, 1.0].
Generated using fpminimax between FLT_MIN and 1. */
- .poly = { V4 (-0x1.55555p-2f), V4 (0x1.99935ep-3f), V4 (-0x1.24051ep-3f),
- V4 (0x1.bd7368p-4f), V4 (-0x1.491f0ep-4f), V4 (0x1.93a2c0p-5f),
- V4 (-0x1.4c3c60p-6f), V4 (0x1.01fd88p-8f) },
- .pi_over_2 = V4 (0x1.921fb6p+0f),
+ .c0 = V4 (-0x1.5554dcp-2), .c1 = 0x1.9978ecp-3,
+ .c2 = V4 (-0x1.230a94p-3), .c3 = 0x1.b4debp-4,
+ .c4 = V4 (-0x1.3550dap-4), .c5 = 0x1.61eebp-5,
+ .c6 = V4 (-0x1.0c17d4p-6), .c7 = 0x1.7ea694p-9,
+#endif
+ .pi_over_2 = V4 (0x3fc90fdb),
+ .neg_one = V4 (-1.0f),
+ .sign_mask = V4 (0x80000000),
};
-#define SignMask v_u32 (0x80000000)
-
-#define P(i) d->poly[i]
-
+#if WANT_SIMD_EXCEPT
#define TinyBound 0x30800000 /* asuint(0x1p-30). */
#define BigBound 0x4e800000 /* asuint(0x1p30). */
-#if WANT_SIMD_EXCEPT
static float32x4_t VPCS_ATTR NOINLINE
special_case (float32x4_t x, float32x4_t y, uint32x4_t special)
{
/* Fast implementation of vector atanf based on
atan(x) ~ shift + z + z^3 * P(z^2) with reduction to [0,1]
- using z=-1/x and shift = pi/2. Maximum observed error is 2.9ulps:
- _ZGVnN4v_atanf (0x1.0468f6p+0) got 0x1.967f06p-1 want 0x1.967fp-1. */
+ using z=-1/x and shift = pi/2. Maximum observed error is 2.02 ulps:
+ _ZGVnN4v_atanf (0x1.03d4cep+0) got 0x1.95ed3ap-1
+ want 0x1.95ed36p-1. */
float32x4_t VPCS_ATTR NOINLINE V_NAME_F1 (atan) (float32x4_t x)
{
const struct data *d = ptr_barrier (&data);
- /* Small cases, infs and nans are supported by our approximation technique,
- but do not set fenv flags correctly. Only trigger special case if we need
- fenv. */
uint32x4_t ix = vreinterpretq_u32_f32 (x);
- uint32x4_t sign = vandq_u32 (ix, SignMask);
+ uint32x4_t sign = vandq_u32 (ix, d->sign_mask);
#if WANT_SIMD_EXCEPT
+ /* Small cases, infs and nans are supported by our approximation technique,
+ but do not set fenv flags correctly. Only trigger special case if we need
+ fenv. */
uint32x4_t ia = vandq_u32 (ix, v_u32 (0x7ff00000));
uint32x4_t special = vcgtq_u32 (vsubq_u32 (ia, v_u32 (TinyBound)),
v_u32 (BigBound - TinyBound));
if (__glibc_unlikely (v_any_u32 (special)))
return special_case (x, x, v_u32 (-1));
#endif
-
/* Argument reduction:
- y := arctan(x) for x < 1
- y := pi/2 + arctan(-1/x) for x > 1
- Hence, use z=-1/a if x>=1, otherwise z=a. */
- uint32x4_t red = vcagtq_f32 (x, v_f32 (1.0));
- /* Avoid dependency in abs(x) in division (and comparison). */
- float32x4_t z = vbslq_f32 (red, vdivq_f32 (v_f32 (1.0f), x), x);
+ y := arctan(x) for |x| < 1
+ y := arctan(-1/x) + pi/2 for x > +1
+ y := arctan(-1/x) - pi/2 for x < -1
+ Hence, use z=-1/a if x>=|-1|, otherwise z=a. */
+ uint32x4_t red = vcagtq_f32 (x, d->neg_one);
+
+ float32x4_t z = vbslq_f32 (red, vdivq_f32 (d->neg_one, x), x);
+
+ /* Shift is calculated as +-pi/2 or 0, depending on the argument case. */
float32x4_t shift = vreinterpretq_f32_u32 (
- vandq_u32 (red, vreinterpretq_u32_f32 (d->pi_over_2)));
- /* Use absolute value only when needed (odd powers of z). */
- float32x4_t az = vbslq_f32 (
- SignMask, vreinterpretq_f32_u32 (vandq_u32 (SignMask, red)), z);
+ vandq_u32 (red, veorq_u32 (d->pi_over_2, sign)));
+
+ float32x4_t z2 = vmulq_f32 (z, z);
+ float32x4_t z3 = vmulq_f32 (z, z2);
+ float32x4_t z4 = vmulq_f32 (z2, z2);
+#if WANT_SIMD_EXCEPT
/* Calculate the polynomial approximation.
Use 2-level Estrin scheme for P(z^2) with deg(P)=7. However,
a standard implementation using z8 creates spurious underflow
in the very last fma (when z^8 is small enough).
- Therefore, we split the last fma into a mul and an fma.
- Horner and single-level Estrin have higher errors that exceed
- threshold. */
- float32x4_t z2 = vmulq_f32 (z, z);
- float32x4_t z4 = vmulq_f32 (z2, z2);
-
+ Therefore, we split the last fma into a mul and an fma. */
float32x4_t y = vfmaq_f32 (
v_pairwise_poly_3_f32 (z2, z4, d->poly), z4,
vmulq_f32 (z4, v_pairwise_poly_3_f32 (z2, z4, d->poly + 4)));
- /* y = shift + z * P(z^2). */
- y = vaddq_f32 (vfmaq_f32 (az, y, vmulq_f32 (z2, az)), shift);
+#else
+ float32x4_t z8 = vmulq_f32 (z4, z4);
+
+ /* Uses an Estrin scheme for polynomial approximation. */
+ float32x4_t odd_coeffs = vld1q_f32 (&d->c1);
+
+ float32x4_t p01 = vfmaq_laneq_f32 (d->c0, z2, odd_coeffs, 0);
+ float32x4_t p23 = vfmaq_laneq_f32 (d->c2, z2, odd_coeffs, 1);
+ float32x4_t p45 = vfmaq_laneq_f32 (d->c4, z2, odd_coeffs, 2);
+ float32x4_t p67 = vfmaq_laneq_f32 (d->c6, z2, odd_coeffs, 3);
- /* y = atan(x) if x>0, -atan(-x) otherwise. */
- y = vreinterpretq_f32_u32 (veorq_u32 (vreinterpretq_u32_f32 (y), sign));
+ float32x4_t p03 = vfmaq_f32 (p01, z4, p23);
+ float32x4_t p47 = vfmaq_f32 (p45, z4, p67);
- return y;
+ float32x4_t y = vfmaq_f32 (p03, z8, p47);
+#endif
+
+ /* y = shift + z * P(z^2). */
+ return vfmaq_f32 (vaddq_f32 (shift, z), z3, y);
}
libmvec_hidden_def (V_NAME_F1 (atan))
HALF_WIDTH_ALIAS_F1 (atan)
<https://www.gnu.org/licenses/>. */
#include "sv_math.h"
-#include "poly_sve_f32.h"
static const struct data
{
- float32_t poly[8];
- float32_t pi_over_2;
+ float32_t c1, c3, c5, c7;
+ float32_t c0, c2, c4, c6;
+ float32_t shift_val, neg_one;
} data = {
/* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on
[2**-128, 1.0]. */
- .poly = { -0x1.55555p-2f, 0x1.99935ep-3f, -0x1.24051ep-3f, 0x1.bd7368p-4f,
- -0x1.491f0ep-4f, 0x1.93a2c0p-5f, -0x1.4c3c60p-6f, 0x1.01fd88p-8f },
- .pi_over_2 = 0x1.921fb6p+0f,
+ .c0 = -0x1.5554dcp-2,
+ .c1 = 0x1.9978ecp-3,
+ .c2 = -0x1.230a94p-3,
+ .c3 = 0x1.b4debp-4,
+ .c4 = -0x1.3550dap-4,
+ .c5 = 0x1.61eebp-5,
+ .c6 = -0x1.0c17d4p-6,
+ .c7 = 0x1.7ea694p-9,
+ /* pi/2, used as a shift value after reduction. */
+ .shift_val = 0x1.921fb54442d18p+0,
+ .neg_one = -1.0f,
};
#define SignMask (0x80000000)
/* Fast implementation of SVE atanf based on
atan(x) ~ shift + z + z^3 * P(z^2) with reduction to [0,1] using
z=-1/x and shift = pi/2.
- Largest observed error is 2.9 ULP, close to +/-1.0:
- _ZGVsMxv_atanf (0x1.0468f6p+0) got -0x1.967f06p-1
- want -0x1.967fp-1. */
+ Largest observed error is 2.12 ULP:
+ _ZGVsMxv_atanf (0x1.03d4cep+0) got 0x1.95ed3ap-1
+ want 0x1.95ed36p-1. */
svfloat32_t SV_NAME_F1 (atan) (svfloat32_t x, const svbool_t pg)
{
const struct data *d = ptr_barrier (&data);
+ svbool_t ptrue = svptrue_b32 ();
/* No need to trigger special case. Small cases, infs and nans
are supported by our approximation technique. */
svuint32_t ix = svreinterpret_u32 (x);
- svuint32_t sign = svand_x (pg, ix, SignMask);
+ svuint32_t sign = svand_x (ptrue, ix, SignMask);
/* Argument reduction:
y := arctan(x) for x < 1
- y := pi/2 + arctan(-1/x) for x > 1
- Hence, use z=-1/a if x>=1, otherwise z=a. */
- svbool_t red = svacgt (pg, x, 1.0f);
- /* Avoid dependency in abs(x) in division (and comparison). */
- svfloat32_t z = svsel (red, svdiv_x (pg, sv_f32 (1.0f), x), x);
- /* Use absolute value only when needed (odd powers of z). */
- svfloat32_t az = svabs_x (pg, z);
- az = svneg_m (az, red, az);
-
- /* Use split Estrin scheme for P(z^2) with deg(P)=7. */
- svfloat32_t z2 = svmul_x (pg, z, z);
- svfloat32_t z4 = svmul_x (pg, z2, z2);
- svfloat32_t z8 = svmul_x (pg, z4, z4);
-
- svfloat32_t y = sv_estrin_7_f32_x (pg, z2, z4, z8, d->poly);
-
- /* y = shift + z + z^3 * P(z^2). */
- svfloat32_t z3 = svmul_x (pg, z2, az);
- y = svmla_x (pg, az, z3, y);
-
- /* Apply shift as indicated by 'red' predicate. */
- y = svadd_m (red, y, sv_f32 (d->pi_over_2));
-
- /* y = atan(x) if x>0, -atan(-x) otherwise. */
- return svreinterpret_f32 (sveor_x (pg, svreinterpret_u32 (y), sign));
+ y := arctan(-1/x) + pi/2 for x > +1
+ y := arctan(-1/x) - pi/2 for x < -1
+ Hence, use z=-1/a if |x|>=|-1|, otherwise z=a. */
+ svbool_t red = svacgt (pg, x, d->neg_one);
+ svfloat32_t z = svsel (red, svdiv_x (pg, sv_f32 (d->neg_one), x), x);
+
+ /* Reinserts the sign bit of the argument to handle the case of x < -1. */
+ svfloat32_t shift = svreinterpret_f32 (
+ sveor_x (red, svreinterpret_u32 (sv_f32 (d->shift_val)), sign));
+
+ svfloat32_t z2 = svmul_x (ptrue, z, z);
+ svfloat32_t z3 = svmul_x (ptrue, z2, z);
+ svfloat32_t z4 = svmul_x (ptrue, z2, z2);
+ svfloat32_t z8 = svmul_x (ptrue, z4, z4);
+
+ svfloat32_t odd_coeffs = svld1rq (ptrue, &d->c1);
+
+ svfloat32_t p01 = svmla_lane (sv_f32 (d->c0), z2, odd_coeffs, 0);
+ svfloat32_t p23 = svmla_lane (sv_f32 (d->c2), z2, odd_coeffs, 1);
+ svfloat32_t p45 = svmla_lane (sv_f32 (d->c4), z2, odd_coeffs, 2);
+ svfloat32_t p67 = svmla_lane (sv_f32 (d->c6), z2, odd_coeffs, 3);
+
+ svfloat32_t p03 = svmla_x (pg, p01, z4, p23);
+ svfloat32_t p47 = svmla_x (pg, p45, z4, p67);
+
+ svfloat32_t y = svmla_x (pg, p03, z8, p47);
+
+ /* shift + z + z^3 * P(z^2). */
+ shift = svadd_m (red, z, shift);
+ return svmla_x (pg, shift, z3, y);
}
projects / thirdparty / glibc.git / commitdiff
