From: Dylan Fleming <Dylan.Fleming@arm.com>
Date: Mon, 19 May 2025 11:36:51 +0000 (+0000)
Subject: AArch64: Optimize inverse trig functions
X-Git-Tag: glibc-2.42~237
X-Git-Url: http://git.ipfire.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1e84509e0041c0a83997aba602a585bb3b8285f0;p=thirdparty%2Fglibc.git

AArch64: Optimize inverse trig functions

Improve performance of Inverse trig functions by altering how coefficients are
loaded.

Performance improvement on Neoverse V1:
SVE     acos   14%
AdvSIMD acos   6%

AdvSIMD asin   6%
SVE     asin   5%
AdvSIMD asinf  2%

AdvSIMD atanf  22%
SVE     atanf  20%
SVE     atan   11%
AdvSIMD atan   5%

SVE     atan2  7%
SVE     atan2f 4%
AdvSIMD atan2f 3%
AdvSIMD atan2  2%

Reviewed-by: Wilco Dijkstra  <Wilco.Dijkstra@arm.com>
---

diff --git a/sysdeps/aarch64/fpu/acos_advsimd.c b/sysdeps/aarch64/fpu/acos_advsimd.c
index 7709b5454f..453f780314 100644
--- a/sysdeps/aarch64/fpu/acos_advsimd.c
+++ b/sysdeps/aarch64/fpu/acos_advsimd.c
@@ -18,24 +18,23 @@
    <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),
 };
 
@@ -63,7 +62,7 @@ special_case (float64x2_t x, float64x2_t y, uint64x2_t special)
 
      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.
 
@@ -71,9 +70,9 @@ special_case (float64x2_t x, float64x2_t y, uint64x2_t special)
 
      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);
@@ -99,13 +98,32 @@ float64x2_t VPCS_ATTR V_NAME_D1 (acos) (float64x2_t x)
   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
diff --git a/sysdeps/aarch64/fpu/acos_sve.c b/sysdeps/aarch64/fpu/acos_sve.c
index 74e2f7df0f..104f0d7805 100644
--- a/sysdeps/aarch64/fpu/acos_sve.c
+++ b/sysdeps/aarch64/fpu/acos_sve.c
@@ -18,20 +18,21 @@
    <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,
 };
 
@@ -42,20 +43,21 @@ static const struct data
 
      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);
@@ -70,24 +72,41 @@ svfloat64_t SV_NAME_D1 (acos) (svfloat64_t x, const svbool_t pg)
   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);
 }
diff --git a/sysdeps/aarch64/fpu/asin_advsimd.c b/sysdeps/aarch64/fpu/asin_advsimd.c
index 414211627e..f74141c845 100644
--- a/sysdeps/aarch64/fpu/asin_advsimd.c
+++ b/sysdeps/aarch64/fpu/asin_advsimd.c
@@ -18,24 +18,23 @@
    <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)
@@ -68,8 +67,8 @@ special_case (float64x2_t x, float64x2_t y, uint64x2_t special)
      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);
@@ -86,7 +85,7 @@ float64x2_t VPCS_ATTR V_NAME_D1 (asin) (float64x2_t x)
     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
@@ -99,7 +98,26 @@ float64x2_t VPCS_ATTR V_NAME_D1 (asin) (float64x2_t x)
   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);
diff --git a/sysdeps/aarch64/fpu/asin_sve.c b/sysdeps/aarch64/fpu/asin_sve.c
index 9314466f58..975f408bee 100644
--- a/sysdeps/aarch64/fpu/asin_sve.c
+++ b/sysdeps/aarch64/fpu/asin_sve.c
@@ -18,45 +18,43 @@
    <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);
@@ -70,17 +68,37 @@ svfloat64_t SV_NAME_D1 (asin) (svfloat64_t x, const svbool_t pg)
   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));
 }
diff --git a/sysdeps/aarch64/fpu/asinf_advsimd.c b/sysdeps/aarch64/fpu/asinf_advsimd.c
index 52c7c0ec6e..013936c2c0 100644
--- a/sysdeps/aarch64/fpu/asinf_advsimd.c
+++ b/sysdeps/aarch64/fpu/asinf_advsimd.c
@@ -18,22 +18,21 @@
    <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.  */
 
@@ -47,11 +46,8 @@ special_case (float32x4_t x, float32x4_t y, uint32x4_t special)
 
 /* 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,
@@ -80,24 +76,31 @@ float32x4_t VPCS_ATTR NOINLINE V_NAME_F1 (asin) (float32x4_t x)
 #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);
diff --git a/sysdeps/aarch64/fpu/atan2_advsimd.c b/sysdeps/aarch64/fpu/atan2_advsimd.c
index 00b4a4f083..a31d52f3ac 100644
--- a/sysdeps/aarch64/fpu/atan2_advsimd.c
+++ b/sysdeps/aarch64/fpu/atan2_advsimd.c
@@ -19,40 +19,38 @@
 
 #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)
@@ -77,10 +75,9 @@ zeroinfnan (uint64x2_t i, const struct data *d)
 }
 
 /* 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);
@@ -101,26 +98,29 @@ float64x2_t VPCS_ATTR V_NAME_D2 (atan2) (float64x2_t y, float64x2_t x)
   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);
@@ -128,45 +128,43 @@ float64x2_t VPCS_ATTR V_NAME_D2 (atan2) (float64x2_t y, float64x2_t x)
   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;
 }
diff --git a/sysdeps/aarch64/fpu/atan2_sve.c b/sysdeps/aarch64/fpu/atan2_sve.c
index 163f61308b..9e2dd249d4 100644
--- a/sysdeps/aarch64/fpu/atan2_sve.c
+++ b/sysdeps/aarch64/fpu/atan2_sve.c
@@ -19,25 +19,25 @@
 
 #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,
@@ -56,15 +56,17 @@ zeroinfnan (svuint64_t i, const svbool_t pg)
 }
 
 /* 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);
@@ -81,32 +83,67 @@ svfloat64_t SV_NAME_D2 (atan2) (svfloat64_t y, svfloat64_t x, const svbool_t 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)))
diff --git a/sysdeps/aarch64/fpu/atan2f_advsimd.c b/sysdeps/aarch64/fpu/atan2f_advsimd.c
index e65406f492..75d873897a 100644
--- a/sysdeps/aarch64/fpu/atan2f_advsimd.c
+++ b/sysdeps/aarch64/fpu/atan2f_advsimd.c
@@ -18,22 +18,22 @@
    <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)
@@ -54,13 +54,13 @@ static inline uint32x4_t
 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);
@@ -81,28 +81,31 @@ float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (atan2) (float32x4_t y, float32x4_t x)
   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);
@@ -110,10 +113,11 @@ float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (atan2) (float32x4_t y, float32x4_t x)
   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)))
     {
diff --git a/sysdeps/aarch64/fpu/atan2f_sve.c b/sysdeps/aarch64/fpu/atan2f_sve.c
index 5f26e2a365..4d9341952d 100644
--- a/sysdeps/aarch64/fpu/atan2f_sve.c
+++ b/sysdeps/aarch64/fpu/atan2f_sve.c
@@ -18,18 +18,18 @@
    <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).  */
@@ -51,12 +51,14 @@ zeroinfnan (svuint32_t i, const svbool_t pg)
 
 /* 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);
@@ -76,29 +78,42 @@ svfloat32_t SV_NAME_F2 (atan2) (svfloat32_t y, svfloat32_t x, const svbool_t pg)
 
   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.  */
 
diff --git a/sysdeps/aarch64/fpu/atan_advsimd.c b/sysdeps/aarch64/fpu/atan_advsimd.c
index f024fd1d74..da0d3715df 100644
--- a/sysdeps/aarch64/fpu/atan_advsimd.c
+++ b/sysdeps/aarch64/fpu/atan_advsimd.c
@@ -18,7 +18,6 @@
    <https://www.gnu.org/licenses/>.  */
 
 #include "v_math.h"
-#include "poly_advsimd_f64.h"
 
 static const struct data
 {
@@ -28,16 +27,16 @@ 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),
 };
 
@@ -47,9 +46,9 @@ static const struct data
 
 /* 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);
@@ -78,59 +77,53 @@ float64x2_t VPCS_ATTR V_NAME_D1 (atan) (float64x2_t x)
      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);
 }
diff --git a/sysdeps/aarch64/fpu/atan_sve.c b/sysdeps/aarch64/fpu/atan_sve.c
index 3880cedff4..a6b0489cf6 100644
--- a/sysdeps/aarch64/fpu/atan_sve.c
+++ b/sysdeps/aarch64/fpu/atan_sve.c
@@ -18,23 +18,26 @@
    <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.  */
@@ -43,15 +46,14 @@ static const struct data
 /* 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);
 
@@ -59,32 +61,60 @@ svfloat64_t SV_NAME_D1 (atan) (svfloat64_t x, const svbool_t pg)
      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);
 }
diff --git a/sysdeps/aarch64/fpu/atanf_advsimd.c b/sysdeps/aarch64/fpu/atanf_advsimd.c
index 472865ed74..817a47ef3e 100644
--- a/sysdeps/aarch64/fpu/atanf_advsimd.c
+++ b/sysdeps/aarch64/fpu/atanf_advsimd.c
@@ -22,26 +22,35 @@
 
 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)
 {
@@ -51,19 +60,20 @@ 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));
@@ -71,41 +81,52 @@ float32x4_t VPCS_ATTR NOINLINE V_NAME_F1 (atan) (float32x4_t x)
   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)
diff --git a/sysdeps/aarch64/fpu/atanf_sve.c b/sysdeps/aarch64/fpu/atanf_sve.c
index 3a98d70c50..6558223e41 100644
--- a/sysdeps/aarch64/fpu/atanf_sve.c
+++ b/sysdeps/aarch64/fpu/atanf_sve.c
@@ -18,18 +18,26 @@
    <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)
@@ -37,43 +45,49 @@ static const struct data
 /* 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);
 }