pl/math: Use Neon indexed fmla/mul with constants stored in neon registers.

Replaced various vfmaq instructions with vfmaq_laneq variants to reduce constant vector loads.

9 files changed, 69 insertions, 70 deletions

diff --git a/pl/math/v_cosh_2u.c b/pl/math/v_cosh_2u.c
index ccb67d7..060d7ec 100644
--- a/pl/math/v_cosh_2u.c
+++ b/pl/math/v_cosh_2u.c
@@ -12,7 +12,7 @@
 static const struct data
 {
   float64x2_t poly[3];
-  float64x2_t inv_ln2, ln2_hi, ln2_lo, shift, thres;
+  float64x2_t inv_ln2, ln2, shift, thres;
   uint64x2_t index_mask, special_bound;
 } data = {
   .poly = { V2 (0x1.fffffffffffd4p-2), V2 (0x1.5555571d6b68cp-3),
@@ -20,8 +20,7 @@ static const struct data
 
   .inv_ln2 = V2 (0x1.71547652b82fep8), /* N/ln2.  */
   /* -ln2/N.  */
-  .ln2_hi = V2 (-0x1.62e42fefa39efp-9),
-  .ln2_lo = V2 (-0x1.abc9e3b39803f3p-64),
+  .ln2 = {-0x1.62e42fefa39efp-9, -0x1.abc9e3b39803f3p-64},
   .shift = V2 (0x1.8p+52),
   .thres = V2 (704.0),
 
@@ -49,8 +48,8 @@ exp_inline (float64x2_t x)
   float64x2_t n = vsubq_f64 (z, d->shift);
 
   /* r = x - n*ln2/N.  */
-  float64x2_t r = vfmaq_f64 (x, d->ln2_hi, n);
-  r = vfmaq_f64 (r, d->ln2_lo, n);
+  float64x2_t r = vfmaq_laneq_f64 (x, n, d->ln2, 0);
+  r = vfmaq_laneq_f64 (r, n, d->ln2, 1);
 
   uint64x2_t e = vshlq_n_u64 (u, 52 - V_EXP_TAIL_TABLE_BITS);
   uint64x2_t i = vandq_u64 (u, d->index_mask);
diff --git a/pl/math/v_exp10f_2u4.c b/pl/math/v_exp10f_2u4.c
index 8a80fbd..0d37dcb 100644
--- a/pl/math/v_exp10f_2u4.c
+++ b/pl/math/v_exp10f_2u4.c
@@ -16,7 +16,8 @@
 static const struct data
 {
   float32x4_t poly[5];
-  float32x4_t shift, log10_2, log2_10_hi, log2_10_lo;
+  float32x4_t log10_2_and_inv, shift;
+
 #if !WANT_SIMD_EXCEPT
   float32x4_t scale_thresh;
 #endif
@@ -29,9 +30,9 @@ static const struct data
   .poly = { V4 (0x1.26bb16p+1f), V4 (0x1.5350d2p+1f), V4 (0x1.04744ap+1f),
 	    V4 (0x1.2d8176p+0f), V4 (0x1.12b41ap-1f) },
   .shift = V4 (0x1.8p23f),
-  .log10_2 = V4 (0x1.a934fp+1),
-  .log2_10_hi = V4 (0x1.344136p-2),
-  .log2_10_lo = V4 (-0x1.ec10cp-27),
+
+  /* Stores constants 1/log10(2), log10(2)_high, log10(2)_low, 0.  */
+  .log10_2_and_inv = { 0x1.a934fp+1, 0x1.344136p-2, -0x1.ec10cp-27, 0 },
 #if !WANT_SIMD_EXCEPT
   .scale_thresh = V4 (ScaleBound)
 #endif
@@ -56,9 +57,9 @@ special_case (float32x4_t x, float32x4_t y, uint32x4_t cmp)
 
 #else
 
-# define SpecialBound 126.0f /* rint (log2 (2^127 / (1 + sqrt (2)))).  */
-# define SpecialOffset v_u32 (0x82000000)
-# define SpecialBias v_u32 (0x7f000000)
+#  define SpecialBound 126.0f /* rint (log2 (2^127 / (1 + sqrt (2)))).  */
+#  define SpecialOffset v_u32 (0x82000000)
+#  define SpecialBias v_u32 (0x7f000000)
 
 static float32x4_t VPCS_ATTR NOINLINE
 special_case (float32x4_t poly, float32x4_t n, uint32x4_t e, uint32x4_t cmp1,
@@ -103,10 +104,10 @@ float32x4_t VPCS_ATTR V_NAME_F1 (exp10) (float32x4_t x)
   /* exp10(x) = 2^n * 10^r = 2^n * (1 + poly (r)),
      with poly(r) in [1/sqrt(2), sqrt(2)] and
      x = r + n * log10 (2), with r in [-log10(2)/2, log10(2)/2].  */
-  float32x4_t z = vfmaq_f32 (d->shift, x, d->log10_2);
+  float32x4_t z = vfmaq_laneq_f32 (d->shift, x, d->log10_2_and_inv, 0);
   float32x4_t n = vsubq_f32 (z, d->shift);
-  float32x4_t r = vfmsq_f32 (x, n, d->log2_10_hi);
-  r = vfmsq_f32 (r, n, d->log2_10_lo);
+  float32x4_t r = vfmsq_laneq_f32 (x, n, d->log10_2_and_inv, 1);
+  r = vfmsq_laneq_f32 (r, n, d->log10_2_and_inv, 2);
   uint32x4_t e = vshlq_n_u32 (vreinterpretq_u32_f32 (z), 23);
 
   float32x4_t scale = vreinterpretq_f32_u32 (vaddq_u32 (e, ExponentBias));
diff --git a/pl/math/v_expf_inline.h b/pl/math/v_expf_inline.h
index 682cbc5..1666837 100644
--- a/pl/math/v_expf_inline.h
+++ b/pl/math/v_expf_inline.h
@@ -14,7 +14,7 @@
 struct v_expf_data
 {
   float32x4_t poly[5];
-  float32x4_t shift, invln2, ln2_hi, ln2_lo;
+  float32x4_t shift, invln2_and_ln2;
 };
 
 /* maxerr: 1.45358 +0.5 ulp.  */
@@ -22,9 +22,8 @@ struct v_expf_data
   {                                                                           \
     .poly = { V4 (0x1.0e4020p-7f), V4 (0x1.573e2ep-5f), V4 (0x1.555e66p-3f),  \
 	      V4 (0x1.fffdb6p-2f), V4 (0x1.ffffecp-1f) },                     \
-                                                                              \
-    .shift = V4 (0x1.8p23f), .invln2 = V4 (0x1.715476p+0f),                   \
-    .ln2_hi = V4 (0x1.62e4p-1f), .ln2_lo = V4 (0x1.7f7d1cp-20f),              \
+    .shift = V4 (0x1.8p23f),                                                  \
+    .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 },   \
   }
 
 #define ExponentBias v_u32 (0x3f800000) /* asuint(1.0f).  */
@@ -40,10 +39,10 @@ v_expf_inline (float32x4_t x, const struct v_expf_data *d)
   /* exp(x) = 2^n (1 + poly(r)), with 1 + poly(r) in [1/sqrt(2),sqrt(2)]
      x = ln2*n + r, with r in [-ln2/2, ln2/2].  */
   float32x4_t n, r, z;
-  z = vfmaq_f32 (d->shift, x, d->invln2);
+  z = vfmaq_laneq_f32 (d->shift, x, d->invln2_and_ln2, 0);
   n = vsubq_f32 (z, d->shift);
-  r = vfmsq_f32 (x, n, d->ln2_hi);
-  r = vfmsq_f32 (r, n, d->ln2_lo);
+  r = vfmsq_laneq_f32 (x, n, d->invln2_and_ln2, 1);
+  r = vfmsq_laneq_f32 (r, n, d->invln2_and_ln2, 2);
   uint32x4_t e = vshlq_n_u32 (vreinterpretq_u32_f32 (z), 23);
   float32x4_t scale = vreinterpretq_f32_u32 (vaddq_u32 (e, ExponentBias));
 
diff --git a/pl/math/v_expm1_2u5.c b/pl/math/v_expm1_2u5.c
index 0c2aa07..f867361 100644
--- a/pl/math/v_expm1_2u5.c
+++ b/pl/math/v_expm1_2u5.c
@@ -13,7 +13,7 @@
 static const struct data
 {
   float64x2_t poly[11];
-  float64x2_t invln2, ln2_lo, ln2_hi, shift;
+  float64x2_t invln2, ln2, shift;
   int64x2_t exponent_bias;
 #if WANT_SIMD_EXCEPT
   uint64x2_t thresh, tiny_bound;
@@ -28,8 +28,7 @@ static const struct data
 	    V2 (0x1.71ddf82db5bb4p-19), V2 (0x1.27e517fc0d54bp-22),
 	    V2 (0x1.af5eedae67435p-26), V2 (0x1.1f143d060a28ap-29) },
   .invln2 = V2 (0x1.71547652b82fep0),
-  .ln2_hi = V2 (0x1.62e42fefa39efp-1),
-  .ln2_lo = V2 (0x1.abc9e3b39803fp-56),
+  .ln2 = { 0x1.62e42fefa39efp-1, 0x1.abc9e3b39803fp-56 },
   .shift = V2 (0x1.8p52),
   .exponent_bias = V2 (0x3ff0000000000000),
 #if WANT_SIMD_EXCEPT
@@ -83,8 +82,8 @@ float64x2_t VPCS_ATTR V_NAME_D1 (expm1) (float64x2_t x)
      where 2^i is exact because i is an integer.  */
   float64x2_t n = vsubq_f64 (vfmaq_f64 (d->shift, d->invln2, x), d->shift);
   int64x2_t i = vcvtq_s64_f64 (n);
-  float64x2_t f = vfmsq_f64 (x, n, d->ln2_hi);
-  f = vfmsq_f64 (f, n, d->ln2_lo);
+  float64x2_t f = vfmsq_laneq_f64 (x, n, d->ln2, 0);
+  f = vfmsq_laneq_f64 (f, n, d->ln2, 1);
 
   /* Approximate expm1(f) using polynomial.
      Taylor expansion for expm1(x) has the form:
diff --git a/pl/math/v_expm1f_1u6.c b/pl/math/v_expm1f_1u6.c
index b755b30..ea6f255 100644
--- a/pl/math/v_expm1f_1u6.c
+++ b/pl/math/v_expm1f_1u6.c
@@ -13,7 +13,8 @@
 static const struct data
 {
   float32x4_t poly[5];
-  float32x4_t invln2, ln2_lo, ln2_hi, shift;
+  float32x4_t invln2_and_ln2;
+  float32x4_t shift;
   int32x4_t exponent_bias;
 #if WANT_SIMD_EXCEPT
   uint32x4_t thresh;
@@ -24,9 +25,8 @@ static const struct data
   /* Generated using fpminimax with degree=5 in [-log(2)/2, log(2)/2].  */
   .poly = { V4 (0x1.fffffep-2), V4 (0x1.5554aep-3), V4 (0x1.555736p-5),
 	    V4 (0x1.12287cp-7), V4 (0x1.6b55a2p-10) },
-  .invln2 = V4 (0x1.715476p+0f),
-  .ln2_hi = V4 (0x1.62e4p-1f),
-  .ln2_lo = V4 (0x1.7f7d1cp-20f),
+  /* Stores constants: invln2, ln2_hi, ln2_lo, 0.  */
+  .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 },
   .shift = V4 (0x1.8p23f),
   .exponent_bias = V4 (0x3f800000),
 #if !WANT_SIMD_EXCEPT
@@ -78,10 +78,11 @@ float32x4_t VPCS_ATTR V_NAME_F1 (expm1) (float32x4_t x)
      and f = x - i * ln2, then f is in [-ln2/2, ln2/2].
      exp(x) - 1 = 2^i * (expm1(f) + 1) - 1
      where 2^i is exact because i is an integer.  */
-  float32x4_t j = vsubq_f32 (vfmaq_f32 (d->shift, d->invln2, x), d->shift);
+  float32x4_t j = vsubq_f32 (
+      vfmaq_laneq_f32 (d->shift, x, d->invln2_and_ln2, 0), d->shift);
   int32x4_t i = vcvtq_s32_f32 (j);
-  float32x4_t f = vfmsq_f32 (x, j, d->ln2_hi);
-  f = vfmsq_f32 (f, j, d->ln2_lo);
+  float32x4_t f = vfmsq_laneq_f32 (x, j, d->invln2_and_ln2, 1);
+  f = vfmsq_laneq_f32 (f, j, d->invln2_and_ln2, 2);
 
   /* Approximate expm1(f) using polynomial.
      Taylor expansion for expm1(x) has the form:
diff --git a/pl/math/v_expm1f_inline.h b/pl/math/v_expm1f_inline.h
index 5a706f5..6ae94c4 100644
--- a/pl/math/v_expm1f_inline.h
+++ b/pl/math/v_expm1f_inline.h
@@ -16,20 +16,19 @@
 struct v_expm1f_data
 {
   float32x4_t poly[5];
-  float32x4_t invln2, ln2_lo, ln2_hi, shift;
+  float32x4_t invln2_and_ln2, shift;
   int32x4_t exponent_bias;
 };
 
 /* Coefficients generated using fpminimax with degree=5 in [-log(2)/2,
- * log(2)/2]. Exponent bias is asuint(1.0f).  */
+   log(2)/2]. Exponent bias is asuint(1.0f).
+   invln2_and_ln2 Stores constants: invln2, ln2_lo, ln2_hi, 0.  */
 #define V_EXPM1F_DATA                                                         \
   {                                                                           \
     .poly = { V4 (0x1.fffffep-2), V4 (0x1.5554aep-3), V4 (0x1.555736p-5),     \
 	      V4 (0x1.12287cp-7), V4 (0x1.6b55a2p-10) },                      \
-                                                                              \
-    .invln2 = V4 (0x1.715476p+0f), .ln2_hi = V4 (0x1.62e4p-1f),               \
-    .ln2_lo = V4 (0x1.7f7d1cp-20f), .shift = V4 (0x1.8p23f),                  \
-    .exponent_bias = V4 (0x3f800000),                                         \
+    .shift = V4 (0x1.8p23f), .exponent_bias = V4 (0x3f800000),                \
+    .invln2_and_ln2 = { 0x1.715476p+0f, 0x1.62e4p-1f, 0x1.7f7d1cp-20f, 0 },   \
   }
 
 static inline float32x4_t
@@ -40,10 +39,11 @@ expm1f_inline (float32x4_t x, const struct v_expm1f_data *d)
      calling routine should handle special values if required.  */
 
   /* Reduce argument: f in [-ln2/2, ln2/2], i is exact.  */
-  float32x4_t j = vsubq_f32 (vfmaq_f32 (d->shift, d->invln2, x), d->shift);
+  float32x4_t j = vsubq_f32 (
+      vfmaq_laneq_f32 (d->shift, x, d->invln2_and_ln2, 0), d->shift);
   int32x4_t i = vcvtq_s32_f32 (j);
-  float32x4_t f = vfmsq_f32 (x, j, d->ln2_hi);
-  f = vfmsq_f32 (f, j, d->ln2_lo);
+  float32x4_t f = vfmsq_laneq_f32 (x, j, d->invln2_and_ln2, 1);
+  f = vfmsq_laneq_f32 (f, j, d->invln2_and_ln2, 2);
 
   /* Approximate expm1(f) with polynomial P, expm1(f) ~= f + f^2 * P(f).
      Uses Estrin scheme, where the main _ZGVnN4v_expm1f routine uses
diff --git a/pl/math/v_sinh_3u.c b/pl/math/v_sinh_3u.c
index b99d6b9..393efa7 100644
--- a/pl/math/v_sinh_3u.c
+++ b/pl/math/v_sinh_3u.c
@@ -13,7 +13,7 @@
 static const struct data
 {
   float64x2_t poly[11];
-  float64x2_t inv_ln2, m_ln2_hi, m_ln2_lo, shift;
+  float64x2_t inv_ln2, m_ln2, shift;
   uint64x2_t halff;
   int64x2_t onef;
 #if WANT_SIMD_EXCEPT
@@ -30,8 +30,7 @@ static const struct data
 	    V2 (0x1.af5eedae67435p-26), V2 (0x1.1f143d060a28ap-29), },
 
   .inv_ln2 = V2 (0x1.71547652b82fep0),
-  .m_ln2_hi = V2 (-0x1.62e42fefa39efp-1),
-  .m_ln2_lo = V2 (-0x1.abc9e3b39803fp-56),
+  .m_ln2 = (float64x2_t) {-0x1.62e42fefa39efp-1, -0x1.abc9e3b39803fp-56},
   .shift = V2 (0x1.8p52),
 
   .halff = V2 (0x3fe0000000000000),
@@ -58,8 +57,8 @@ expm1_inline (float64x2_t x)
      and   f = x - i * ln2 (f in [-ln2/2, ln2/2]).  */
   float64x2_t j = vsubq_f64 (vfmaq_f64 (d->shift, d->inv_ln2, x), d->shift);
   int64x2_t i = vcvtq_s64_f64 (j);
-  float64x2_t f = vfmaq_f64 (x, j, d->m_ln2_hi);
-  f = vfmaq_f64 (f, j, d->m_ln2_lo);
+  float64x2_t f = vfmaq_laneq_f64 (x, j, d->m_ln2, 0);
+  f = vfmaq_laneq_f64 (f, j, d->m_ln2, 1);
   /* Approximate expm1(f) using polynomial.  */
   float64x2_t f2 = vmulq_f64 (f, f);
   float64x2_t f4 = vmulq_f64 (f2, f2);
diff --git a/pl/math/v_tan_3u5.c b/pl/math/v_tan_3u5.c
index 90612f1..b390c92 100644
--- a/pl/math/v_tan_3u5.c
+++ b/pl/math/v_tan_3u5.c
@@ -13,7 +13,7 @@
 static const struct data
 {
   float64x2_t poly[9];
-  float64x2_t half_pi_hi, half_pi_lo, two_over_pi, shift;
+  float64x2_t half_pi, two_over_pi, shift;
 #if !WANT_SIMD_EXCEPT
   float64x2_t range_val;
 #endif
@@ -24,8 +24,7 @@ static const struct data
 	    V2 (0x1.226e5e5ecdfa3p-7), V2 (0x1.d6c7ddbf87047p-9),
 	    V2 (0x1.7ea75d05b583ep-10), V2 (0x1.289f22964a03cp-11),
 	    V2 (0x1.4e4fd14147622p-12) },
-  .half_pi_hi = V2 (0x1.921fb54442d18p0),
-  .half_pi_lo = V2 (0x1.1a62633145c07p-54),
+  .half_pi = { 0x1.921fb54442d18p0, 0x1.1a62633145c07p-54 },
   .two_over_pi = V2 (0x1.45f306dc9c883p-1),
   .shift = V2 (0x1.8p52),
 #if !WANT_SIMD_EXCEPT
@@ -51,10 +50,10 @@ special_case (float64x2_t x)
 float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x)
 {
   const struct data *dat = ptr_barrier (&data);
-  /* Our argument reduction cannot calculate q with sufficient accuracy for very
-     large inputs. Fall back to scalar routine for all lanes if any are too
-     large, or Inf/NaN. If fenv exceptions are expected, also fall back for tiny
-     input to avoid underflow.  */
+  /* Our argument reduction cannot calculate q with sufficient accuracy for
+     very large inputs. Fall back to scalar routine for all lanes if any are
+     too large, or Inf/NaN. If fenv exceptions are expected, also fall back for
+     tiny input to avoid underflow.  */
 #if WANT_SIMD_EXCEPT
   uint64x2_t iax = vreinterpretq_u64_f64 (vabsq_f64 (x));
   /* iax - tiny_bound > range_val - tiny_bound.  */
@@ -72,8 +71,8 @@ float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x)
   /* Use q to reduce x to r in [-pi/4, pi/4], by:
      r = x - q * pi/2, in extended precision.  */
   float64x2_t r = x;
-  r = vfmsq_f64 (r, q, dat->half_pi_hi);
-  r = vfmsq_f64 (r, q, dat->half_pi_lo);
+  r = vfmsq_laneq_f64 (r, q, dat->half_pi, 0);
+  r = vfmsq_laneq_f64 (r, q, dat->half_pi, 1);
   /* Further reduce r to [-pi/8, pi/8], to be reconstructed using double angle
      formula.  */
   r = vmulq_n_f64 (r, 0.5);
@@ -96,7 +95,8 @@ float64x2_t VPCS_ATTR V_NAME_D1 (tan) (float64x2_t x)
      and reciprocity around pi/2:
      tan(x) = 1 / (tan(pi/2 - x))
      to assemble result using change-of-sign and conditional selection of
-     numerator/denominator, dependent on odd/even-ness of q (hence quadrant). */
+     numerator/denominator, dependent on odd/even-ness of q (hence quadrant).
+   */
   float64x2_t n = vfmaq_f64 (v_f64 (-1), p, p);
   float64x2_t d = vaddq_f64 (p, p);
 
diff --git a/pl/math/v_tanf_3u5.c b/pl/math/v_tanf_3u5.c
index e1f0a30..c3dd616 100644
--- a/pl/math/v_tanf_3u5.c
+++ b/pl/math/v_tanf_3u5.c
@@ -13,7 +13,8 @@
 static const struct data
 {
   float32x4_t poly[6];
-  float32x4_t neg_half_pi_1, neg_half_pi_2, neg_half_pi_3, two_over_pi, shift;
+  float32x4_t pi_consts;
+  float32x4_t shift;
 #if !WANT_SIMD_EXCEPT
   float32x4_t range_val;
 #endif
@@ -21,10 +22,9 @@ static const struct data
   /* Coefficients generated using FPMinimax.  */
   .poly = { V4 (0x1.55555p-2f), V4 (0x1.11166p-3f), V4 (0x1.b88a78p-5f),
 	    V4 (0x1.7b5756p-6f), V4 (0x1.4ef4cep-8f), V4 (0x1.0e1e74p-7f) },
-  .neg_half_pi_1 = V4 (-0x1.921fb6p+0f),
-  .neg_half_pi_2 = V4 (0x1.777a5cp-25f),
-  .neg_half_pi_3 = V4 (0x1.ee59dap-50f),
-  .two_over_pi = V4 (0x1.45f306p-1f),
+  /* Stores constants: (-pi/2)_high, (-pi/2)_mid, (-pi/2)_low, and 2/pi.  */
+  .pi_consts
+  = { -0x1.921fb6p+0f, 0x1.777a5cp-25f, 0x1.ee59dap-50f, 0x1.45f306p-1f },
   .shift = V4 (0x1.8p+23f),
 #if !WANT_SIMD_EXCEPT
   .range_val = V4 (0x1p15f),
@@ -48,10 +48,11 @@ eval_poly (float32x4_t z, const struct data *d)
 {
   float32x4_t z2 = vmulq_f32 (z, z);
 #if WANT_SIMD_EXCEPT
-  /* Tiny z (<= 0x1p-31) will underflow when calculating z^4. If fp exceptions
-     are to be triggered correctly, sidestep this by fixing such lanes to 0.  */
+  /* Tiny z (<= 0x1p-31) will underflow when calculating z^4.
+     If fp exceptions are to be triggered correctly,
+     sidestep this by fixing such lanes to 0.  */
   uint32x4_t will_uflow
-    = vcleq_u32 (vreinterpretq_u32_f32 (vabsq_f32 (z)), TinyBound);
+      = vcleq_u32 (vreinterpretq_u32_f32 (vabsq_f32 (z)), TinyBound);
   if (unlikely (v_any_u32 (will_uflow)))
     z2 = vbslq_f32 (will_uflow, v_f32 (0), z2);
 #endif
@@ -84,16 +85,16 @@ float32x4_t VPCS_ATTR V_NAME_F1 (tan) (float32x4_t x)
 #endif
 
   /* n = rint(x/(pi/2)).  */
-  float32x4_t q = vfmaq_f32 (d->shift, d->two_over_pi, x);
+  float32x4_t q = vfmaq_laneq_f32 (d->shift, x, d->pi_consts, 3);
   float32x4_t n = vsubq_f32 (q, d->shift);
   /* Determine if x lives in an interval, where |tan(x)| grows to infinity.  */
   uint32x4_t pred_alt = vtstq_u32 (vreinterpretq_u32_f32 (q), v_u32 (1));
 
   /* r = x - n * (pi/2)  (range reduction into -pi./4 .. pi/4).  */
   float32x4_t r;
-  r = vfmaq_f32 (x, d->neg_half_pi_1, n);
-  r = vfmaq_f32 (r, d->neg_half_pi_2, n);
-  r = vfmaq_f32 (r, d->neg_half_pi_3, n);
+  r = vfmaq_laneq_f32 (x, n, d->pi_consts, 0);
+  r = vfmaq_laneq_f32 (r, n, d->pi_consts, 1);
+  r = vfmaq_laneq_f32 (r, n, d->pi_consts, 2);
 
   /* If x lives in an interval, where |tan(x)|
      - is finite, then use a polynomial approximation of the form