Add new pow implementation

The algorithm is exp(y * log(x)), where log(x) is computed with about
1.8*2^-66 relative error, returning the result in two doubles, and the
exp part uses the same algorithm (and lookup tables) as exp, but takes
the input as two doubles and a sign (to handle negative bases with odd
integer exponent).

There is separate code path when fma is not available but the worst
case error is about 0.67 ULP in both cases.  The lookup table and
consts for log are 4224 bytes, the code is 1196 bytes.  The non-nearest
rounding error is less than 1 ULP.

Improvements on Cortex-A72 compared to current glibc master:
latency: 1.8x
thruput: 2.5x

1 files changed, 20 insertions, 0 deletions

diff --git a/math/math_config.h b/math/math_config.h
index 28b7d26..3625285 100644
--- a/math/math_config.h
+++ b/math/math_config.h
@@ -139,6 +139,15 @@ issignalingf_inline (float x)
   return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000;
 }
 
+static inline int
+issignaling_inline (double x)
+{
+  uint64_t ix = asuint64 (x);
+  if (!IEEE_754_2008_SNAN)
+    return (ix & 0x7ff8000000000000) == 0x7ff8000000000000;
+  return 2 * (ix ^ 0x0008000000000000) > 2 * 0x7ff8000000000000ULL;
+}
+
 /* Force the evaluation of a floating-point expression for its side-effect.  */
 #if __aarch64__ && __GNUC__
 static inline void
@@ -323,4 +332,15 @@ extern const struct log2_data {
 #endif
 } __log2_data HIDDEN;
 
+#define POW_LOG_TABLE_BITS 8
+#define POW_LOG_POLY_ORDER 7
+#define POW_LOG_POLY1_ORDER 9
+extern const struct pow_log_data {
+  double ln2hi;
+  double ln2lo;
+  double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1.  */
+  double poly1[POW_LOG_POLY1_ORDER - 1];
+  struct {double invc, logc;} tab[1 << POW_LOG_TABLE_BITS];
+} __pow_log_data HIDDEN;
+
 #endif