[PATCH 7/8] math: Sync pow with Arm Optimized Routines

[Adhemerval Zanella]

Sync with commit 712aa21 for pow.c and 189dfef for pow_log_data.c.
No functional change.

Checked on x86_64-linux-gnu, i686-linux-gnu, and aarch64-linux-gnu.
---
 sysdeps/ieee754/dbl-64/e_pow.c          | 24 +++++++++++++++++++-----
 sysdeps/ieee754/dbl-64/e_pow_log_data.c |  2 +-
 sysdeps/ieee754/dbl-64/math_config.h    |  2 +-
 3 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index 56c393e917..d0f9936ae2 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -73,7 +73,7 @@ log_inline (uint64_t ix, double *tail)
   /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
      |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible.  */
 #ifdef __FP_FAST_FMA
-  r = __builtin_fma (z, invc, -1.0);
+  r = fma (z, invc, -1.0);
 #else
   /* Split z such that rhi, rlo and rhi*rhi are exact and |rlo| <= |r|.  */
   double zhi = asdouble ((iz + (1ULL << 31)) & (-1ULL << 32));
@@ -97,7 +97,7 @@ log_inline (uint64_t ix, double *tail)
   /* k*Ln2 + log(c) + r + A[0]*r*r.  */
 #ifdef __FP_FAST_FMA
   hi = t2 + ar2;
-  lo3 = __builtin_fma (ar, r, -ar2);
+  lo3 = fma (ar, r, -ar2);
   lo4 = t2 - hi + ar2;
 #else
   double arhi = A[0] * rhi;
@@ -107,8 +107,10 @@ log_inline (uint64_t ix, double *tail)
   lo4 = t2 - hi + arhi2;
 #endif
   /* p = log1p(r) - r - A[0]*r*r.  */
+#if POW_LOG_POLY_ORDER == 8
   p = (ar3
        * (A[1] + r * A[2] + ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r * A[6]))));
+#endif
   lo = lo1 + lo2 + lo3 + lo4 + p;
   y = hi + lo;
   *tail = hi - y + lo;
@@ -221,9 +223,13 @@ exp_inline (double x, double xtail, uint32_t sign_bias)
   /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
   z = InvLn2N * x;
 #if TOINT_INTRINSICS
-  /* z - kd is in [-0.5, 0.5] in all rounding modes.  */
   kd = roundtoint (z);
   ki = converttoint (z);
+#elif EXP_USE_TOINT_NARROW
+  /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes.  */
+  kd = eval_as_double (z + Shift);
+  ki = asuint64 (kd) >> 16;
+  kd = (double) (int32_t) ki;
 #else
   /* z - kd is in [-1, 1] in non-nearest rounding modes.  */
   kd = math_narrow_eval (z + Shift);
@@ -244,7 +250,13 @@ exp_inline (double x, double xtail, uint32_t sign_bias)
   r2 = r * r;
   /* Without fma the worst case error is 0.25/N ulp larger.  */
   /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp.  */
+#if EXP_POLY_ORDER == 4
+  tmp = tail + r + r2 * C2 + r * r2 * (C3 + r * C4);
+#elif EXP_POLY_ORDER == 5
   tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
+#elif EXP_POLY_ORDER == 6
+  tmp = tail + r + r2 * (0.5 + r * C3) + r2 * r2 * (C4 + r * C5 + r2 * C6);
+#endif
   if (__glibc_unlikely (abstop == 0))
     return specialcase (tmp, sbits, ki);
   scale = asdouble (sbits);
@@ -360,7 +372,9 @@ __pow (double x, double y)
       if (topx == 0)
 	{
 	  /* Normalize subnormal x so exponent becomes negative.  */
-	  ix = asuint64 (x * 0x1p52);
+	  /* Without the barrier some versions of clang evalutate the mul
+	     unconditionally causing spurious overflow exceptions.  */
+	  ix = asuint64 (math_opt_barrier (x) * 0x1p52);
 	  ix &= 0x7fffffffffffffff;
 	  ix -= 52ULL << 52;
 	}
@@ -371,7 +385,7 @@ __pow (double x, double y)
   double ehi, elo;
 #ifdef __FP_FAST_FMA
   ehi = y * hi;
-  elo = y * lo + __builtin_fma (y, hi, -ehi);
+  elo = y * lo + fma (y, hi, -ehi);
 #else
   double yhi = asdouble (iy & -1ULL << 27);
   double ylo = y - yhi;
diff --git a/sysdeps/ieee754/dbl-64/e_pow_log_data.c b/sysdeps/ieee754/dbl-64/e_pow_log_data.c
index b9fa605ccc..3f651518cf 100644
--- a/sysdeps/ieee754/dbl-64/e_pow_log_data.c
+++ b/sysdeps/ieee754/dbl-64/e_pow_log_data.c
@@ -61,7 +61,7 @@ error and the interval for z is selected such that near x == 1, where log(x)
 is tiny, large cancellation error is avoided in logc + poly(z/c - 1).  */
 .tab = {
 #if N == 128
-#define A(a,b,c) {a,0,b,c},
+#define A(a, b, c) {a, 0, b, c},
 A(0x1.6a00000000000p+0, -0x1.62c82f2b9c800p-2, 0x1.ab42428375680p-48)
 A(0x1.6800000000000p+0, -0x1.5d1bdbf580800p-2, -0x1.ca508d8e0f720p-46)
 A(0x1.6600000000000p+0, -0x1.5767717455800p-2, -0x1.362a4d5b6506dp-45)
diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/math_config.h
index a3506d77eb..8bdd0c3349 100644
--- a/sysdeps/ieee754/dbl-64/math_config.h
+++ b/sysdeps/ieee754/dbl-64/math_config.h
@@ -303,6 +303,6 @@ extern const struct pow_log_data
   /* Note: the pad field is unused, but allows slightly faster indexing.  */
   /* See e_pow_log_data.c for details.  */
   struct {double invc, pad, logc, logctail;} tab[1 << POW_LOG_TABLE_BITS];
-} __pow_log_data attribute_hidden;
+} __pow_log_data __attribute__ ((__aligned__ (16))) attribute_hidden;
 
 #endif
-- 
2.43.0