math: new pow

from https://github.com/ARM-software/optimized-routines,
commit 04884bd04eac4b251da4026900010ea7d8850edc

The underflow exception is signaled if the result is in the subnormal
range even if the result is exact.

code size change: +3421 bytes.
benchmark on x86_64 before, after, speedup:

-Os:
   pow rthruput: 102.96 ns/call 33.38 ns/call 3.08x
    pow latency: 144.37 ns/call 54.75 ns/call 2.64x
-O3:
   pow rthruput:  98.91 ns/call 32.79 ns/call 3.02x
    pow latency: 138.74 ns/call 53.78 ns/call 2.58x

3 files changed, 520 insertions, 303 deletions

diff --git a/src/math/pow.c b/src/math/pow.c
index 3ddc1b6f..694c2ef6 100644
--- a/src/math/pow.c
+++ b/src/math/pow.c
@@ -1,328 +1,343 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/e_pow.c */
 /*
- * ====================================================
- * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ * Double-precision x^y function.
  *
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* pow(x,y) return x**y
- *
- *                    n
- * Method:  Let x =  2   * (1+f)
- *      1. Compute and return log2(x) in two pieces:
- *              log2(x) = w1 + w2,
- *         where w1 has 53-24 = 29 bit trailing zeros.
- *      2. Perform y*log2(x) = n+y' by simulating muti-precision
- *         arithmetic, where |y'|<=0.5.
- *      3. Return x**y = 2**n*exp(y'*log2)
- *
- * Special cases:
- *      1.  (anything) ** 0  is 1
- *      2.  1 ** (anything)  is 1
- *      3.  (anything except 1) ** NAN is NAN
- *      4.  NAN ** (anything except 0) is NAN
- *      5.  +-(|x| > 1) **  +INF is +INF
- *      6.  +-(|x| > 1) **  -INF is +0
- *      7.  +-(|x| < 1) **  +INF is +0
- *      8.  +-(|x| < 1) **  -INF is +INF
- *      9.  -1          ** +-INF is 1
- *      10. +0 ** (+anything except 0, NAN)               is +0
- *      11. -0 ** (+anything except 0, NAN, odd integer)  is +0
- *      12. +0 ** (-anything except 0, NAN)               is +INF, raise divbyzero
- *      13. -0 ** (-anything except 0, NAN, odd integer)  is +INF, raise divbyzero
- *      14. -0 ** (+odd integer) is -0
- *      15. -0 ** (-odd integer) is -INF, raise divbyzero
- *      16. +INF ** (+anything except 0,NAN) is +INF
- *      17. +INF ** (-anything except 0,NAN) is +0
- *      18. -INF ** (+odd integer) is -INF
- *      19. -INF ** (anything) = -0 ** (-anything), (anything except odd integer)
- *      20. (anything) ** 1 is (anything)
- *      21. (anything) ** -1 is 1/(anything)
- *      22. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
- *      23. (-anything except 0 and inf) ** (non-integer) is NAN
- *
- * Accuracy:
- *      pow(x,y) returns x**y nearly rounded. In particular
- *                      pow(integer,integer)
- *      always returns the correct integer provided it is
- *      representable.
- *
- * Constants :
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
- * to produce the hexadecimal values shown.
+ * Copyright (c) 2018, Arm Limited.
+ * SPDX-License-Identifier: MIT
  */
 
+#include <math.h>
+#include <stdint.h>
 #include "libm.h"
+#include "exp_data.h"
+#include "pow_data.h"
 
-static const double
-bp[]   = {1.0, 1.5,},
-dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
-dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
-two53  =  9007199254740992.0, /* 0x43400000, 0x00000000 */
-huge   =  1.0e300,
-tiny   =  1.0e-300,
-/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
-L1 =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
-L2 =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
-L3 =  3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
-L4 =  2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
-L5 =  2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
-L6 =  2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
-P1 =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
-P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
-P3 =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
-P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
-P5 =  4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
-lg2     =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
-lg2_h   =  6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
-lg2_l   = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
-ovt     =  8.0085662595372944372e-017, /* -(1024-log2(ovfl+.5ulp)) */
-cp      =  9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
-cp_h    =  9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
-cp_l    = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
-ivln2   =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
-ivln2_h =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
-ivln2_l =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
+/*
+Worst-case error: 0.54 ULP (~= ulperr_exp + 1024*Ln2*relerr_log*2^53)
+relerr_log: 1.3 * 2^-68 (Relative error of log, 1.5 * 2^-68 without fma)
+ulperr_exp: 0.509 ULP (ULP error of exp, 0.511 ULP without fma)
+*/
 
-double pow(double x, double y)
+#define T __pow_log_data.tab
+#define A __pow_log_data.poly
+#define Ln2hi __pow_log_data.ln2hi
+#define Ln2lo __pow_log_data.ln2lo
+#define N (1 << POW_LOG_TABLE_BITS)
+#define OFF 0x3fe6955500000000
+
+/* Top 12 bits of a double (sign and exponent bits).  */
+static inline uint32_t top12(double x)
 {
-	double z,ax,z_h,z_l,p_h,p_l;
-	double y1,t1,t2,r,s,t,u,v,w;
-	int32_t i,j,k,yisint,n;
-	int32_t hx,hy,ix,iy;
-	uint32_t lx,ly;
+	return asuint64(x) >> 52;
+}
 
-	EXTRACT_WORDS(hx, lx, x);
-	EXTRACT_WORDS(hy, ly, y);
-	ix = hx & 0x7fffffff;
-	iy = hy & 0x7fffffff;
+/* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
+   additional 15 bits precision.  IX is the bit representation of x, but
+   normalized in the subnormal range using the sign bit for the exponent.  */
+static inline double_t log_inline(uint64_t ix, double_t *tail)
+{
+	/* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
+	double_t z, r, y, invc, logc, logctail, kd, hi, t1, t2, lo, lo1, lo2, p;
+	uint64_t iz, tmp;
+	int k, i;
 
-	/* x**0 = 1, even if x is NaN */
-	if ((iy|ly) == 0)
-		return 1.0;
-	/* 1**y = 1, even if y is NaN */
-	if (hx == 0x3ff00000 && lx == 0)
-		return 1.0;
-	/* NaN if either arg is NaN */
-	if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0) ||
-	    iy > 0x7ff00000 || (iy == 0x7ff00000 && ly != 0))
-		return x + y;
+	/* x = 2^k z; where z is in range [OFF,2*OFF) and exact.
+	   The range is split into N subintervals.
+	   The ith subinterval contains z and c is near its center.  */
+	tmp = ix - OFF;
+	i = (tmp >> (52 - POW_LOG_TABLE_BITS)) % N;
+	k = (int64_t)tmp >> 52; /* arithmetic shift */
+	iz = ix - (tmp & 0xfffULL << 52);
+	z = asdouble(iz);
+	kd = (double_t)k;
 
-	/* determine if y is an odd int when x < 0
-	 * yisint = 0       ... y is not an integer
-	 * yisint = 1       ... y is an odd int
-	 * yisint = 2       ... y is an even int
-	 */
-	yisint = 0;
-	if (hx < 0) {
-		if (iy >= 0x43400000)
-			yisint = 2; /* even integer y */
-		else if (iy >= 0x3ff00000) {
-			k = (iy>>20) - 0x3ff;  /* exponent */
-			if (k > 20) {
-				uint32_t j = ly>>(52-k);
-				if ((j<<(52-k)) == ly)
-					yisint = 2 - (j&1);
-			} else if (ly == 0) {
-				uint32_t j = iy>>(20-k);
-				if ((j<<(20-k)) == iy)
-					yisint = 2 - (j&1);
-			}
-		}
-	}
+	/* log(x) = k*Ln2 + log(c) + log1p(z/c-1).  */
+	invc = T[i].invc;
+	logc = T[i].logc;
+	logctail = T[i].logctail;
 
-	/* special value of y */
-	if (ly == 0) {
-		if (iy == 0x7ff00000) {  /* y is +-inf */
-			if (((ix-0x3ff00000)|lx) == 0)  /* (-1)**+-inf is 1 */
-				return 1.0;
-			else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
-				return hy >= 0 ? y : 0.0;
-			else                       /* (|x|<1)**+-inf = 0,inf */
-				return hy >= 0 ? 0.0 : -y;
-		}
-		if (iy == 0x3ff00000) {    /* y is +-1 */
-			if (hy >= 0)
-				return x;
-			y = 1/x;
-#if FLT_EVAL_METHOD!=0
-			{
-				union {double f; uint64_t i;} u = {y};
-				uint64_t i = u.i & -1ULL/2;
-				if (i>>52 == 0 && (i&(i-1)))
-					FORCE_EVAL((float)y);
-			}
+	/* 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.  */
+#if __FP_FAST_FMA
+	r = __builtin_fma(z, invc, -1.0);
+#else
+	/* Split z such that rhi, rlo and rhi*rhi are exact and |rlo| <= |r|.  */
+	double_t zhi = asdouble((iz + (1ULL << 31)) & (-1ULL << 32));
+	double_t zlo = z - zhi;
+	double_t rhi = zhi * invc - 1.0;
+	double_t rlo = zlo * invc;
+	r = rhi + rlo;
 #endif
-			return y;
-		}
-		if (hy == 0x40000000)    /* y is 2 */
-			return x*x;
-		if (hy == 0x3fe00000) {  /* y is 0.5 */
-			if (hx >= 0)     /* x >= +0 */
-				return sqrt(x);
-		}
+
+	/* k*Ln2 + log(c) + r.  */
+	t1 = kd * Ln2hi + logc;
+	t2 = t1 + r;
+	lo1 = kd * Ln2lo + logctail;
+	lo2 = t1 - t2 + r;
+
+	/* Evaluation is optimized assuming superscalar pipelined execution.  */
+	double_t ar, ar2, ar3, lo3, lo4;
+	ar = A[0] * r; /* A[0] = -0.5.  */
+	ar2 = r * ar;
+	ar3 = r * ar2;
+	/* k*Ln2 + log(c) + r + A[0]*r*r.  */
+#if __FP_FAST_FMA
+	hi = t2 + ar2;
+	lo3 = __builtin_fma(ar, r, -ar2);
+	lo4 = t2 - hi + ar2;
+#else
+	double_t arhi = A[0] * rhi;
+	double_t arhi2 = rhi * arhi;
+	hi = t2 + arhi2;
+	lo3 = rlo * (ar + arhi);
+	lo4 = t2 - hi + arhi2;
+#endif
+	/* p = log1p(r) - r - A[0]*r*r.  */
+	p = (ar3 * (A[1] + r * A[2] +
+		    ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r * A[6]))));
+	lo = lo1 + lo2 + lo3 + lo4 + p;
+	y = hi + lo;
+	*tail = hi - y + lo;
+	return y;
+}
+
+#undef N
+#undef T
+#define N (1 << EXP_TABLE_BITS)
+#define InvLn2N __exp_data.invln2N
+#define NegLn2hiN __exp_data.negln2hiN
+#define NegLn2loN __exp_data.negln2loN
+#define Shift __exp_data.shift
+#define T __exp_data.tab
+#define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
+#define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
+#define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
+#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
+#define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
+
+/* Handle cases that may overflow or underflow when computing the result that
+   is scale*(1+TMP) without intermediate rounding.  The bit representation of
+   scale is in SBITS, however it has a computed exponent that may have
+   overflown into the sign bit so that needs to be adjusted before using it as
+   a double.  (int32_t)KI is the k used in the argument reduction and exponent
+   adjustment of scale, positive k here means the result may overflow and
+   negative k means the result may underflow.  */
+static inline double specialcase(double_t tmp, uint64_t sbits, uint64_t ki)
+{
+	double_t scale, y;
+
+	if ((ki & 0x80000000) == 0) {
+		/* k > 0, the exponent of scale might have overflowed by <= 460.  */
+		sbits -= 1009ull << 52;
+		scale = asdouble(sbits);
+		y = 0x1p1009 * (scale + scale * tmp);
+		return eval_as_double(y);
+	}
+	/* k < 0, need special care in the subnormal range.  */
+	sbits += 1022ull << 52;
+	/* Note: sbits is signed scale.  */
+	scale = asdouble(sbits);
+	y = scale + scale * tmp;
+	if (fabs(y) < 1.0) {
+		/* Round y to the right precision before scaling it into the subnormal
+		   range to avoid double rounding that can cause 0.5+E/2 ulp error where
+		   E is the worst-case ulp error outside the subnormal range.  So this
+		   is only useful if the goal is better than 1 ulp worst-case error.  */
+		double_t hi, lo, one = 1.0;
+		if (y < 0.0)
+			one = -1.0;
+		lo = scale - y + scale * tmp;
+		hi = one + y;
+		lo = one - hi + y + lo;
+		y = eval_as_double(hi + lo) - one;
+		/* Fix the sign of 0.  */
+		if (y == 0.0)
+			y = asdouble(sbits & 0x8000000000000000);
+		/* The underflow exception needs to be signaled explicitly.  */
+		fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022);
 	}
+	y = 0x1p-1022 * y;
+	return eval_as_double(y);
+}
 
-	ax = fabs(x);
-	/* special value of x */
-	if (lx == 0) {
-		if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { /* x is +-0,+-inf,+-1 */
-			z = ax;
-			if (hy < 0)   /* z = (1/|x|) */
-				z = 1.0/z;
-			if (hx < 0) {
-				if (((ix-0x3ff00000)|yisint) == 0) {
-					z = (z-z)/(z-z); /* (-1)**non-int is NaN */
-				} else if (yisint == 1)
-					z = -z;          /* (x<0)**odd = -(|x|**odd) */
-			}
-			return z;
+#define SIGN_BIAS (0x800 << EXP_TABLE_BITS)
+
+/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+   The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1.  */
+static inline double exp_inline(double_t x, double_t xtail, uint32_t sign_bias)
+{
+	uint32_t abstop;
+	uint64_t ki, idx, top, sbits;
+	/* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
+	double_t kd, z, r, r2, scale, tail, tmp;
+
+	abstop = top12(x) & 0x7ff;
+	if (predict_false(abstop - top12(0x1p-54) >=
+			  top12(512.0) - top12(0x1p-54))) {
+		if (abstop - top12(0x1p-54) >= 0x80000000) {
+			/* Avoid spurious underflow for tiny x.  */
+			/* Note: 0 is common input.  */
+			double_t one = WANT_ROUNDING ? 1.0 + x : 1.0;
+			return sign_bias ? -one : one;
+		}
+		if (abstop >= top12(1024.0)) {
+			/* Note: inf and nan are already handled.  */
+			if (asuint64(x) >> 63)
+				return __math_uflow(sign_bias);
+			else
+				return __math_oflow(sign_bias);
 		}
+		/* Large x is special cased below.  */
+		abstop = 0;
 	}
 
-	s = 1.0; /* sign of result */
-	if (hx < 0) {
-		if (yisint == 0) /* (x<0)**(non-int) is NaN */
-			return (x-x)/(x-x);
-		if (yisint == 1) /* (x<0)**(odd int) */
-			s = -1.0;
-	}
+	/* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
+	/* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
+	z = InvLn2N * x;
+#if TOINT_INTRINSICS
+	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_t)(int32_t)ki;
+#else
+	/* z - kd is in [-1, 1] in non-nearest rounding modes.  */
+	kd = eval_as_double(z + Shift);
+	ki = asuint64(kd);
+	kd -= Shift;
+#endif
+	r = x + kd * NegLn2hiN + kd * NegLn2loN;
+	/* The code assumes 2^-200 < |xtail| < 2^-8/N.  */
+	r += xtail;
+	/* 2^(k/N) ~= scale * (1 + tail).  */
+	idx = 2 * (ki % N);
+	top = (ki + sign_bias) << (52 - EXP_TABLE_BITS);
+	tail = asdouble(T[idx]);
+	/* This is only a valid scale when -1023*N < k < 1024*N.  */
+	sbits = T[idx + 1] + top;
+	/* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
+	/* Evaluation is optimized assuming superscalar pipelined execution.  */
+	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.  */
+	tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
+	if (predict_false(abstop == 0))
+		return specialcase(tmp, sbits, ki);
+	scale = asdouble(sbits);
+	/* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+	   is no spurious underflow here even without fma.  */
+	return eval_as_double(scale + scale * tmp);
+}
 
-	/* |y| is huge */
-	if (iy > 0x41e00000) { /* if |y| > 2**31 */
-		if (iy > 0x43f00000) {  /* if |y| > 2**64, must o/uflow */
-			if (ix <= 0x3fefffff)
-				return hy < 0 ? huge*huge : tiny*tiny;
-			if (ix >= 0x3ff00000)
-				return hy > 0 ? huge*huge : tiny*tiny;
+/* Returns 0 if not int, 1 if odd int, 2 if even int.  The argument is
+   the bit representation of a non-zero finite floating-point value.  */
+static inline int checkint(uint64_t iy)
+{
+	int e = iy >> 52 & 0x7ff;
+	if (e < 0x3ff)
+		return 0;
+	if (e > 0x3ff + 52)
+		return 2;
+	if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
+		return 0;
+	if (iy & (1ULL << (0x3ff + 52 - e)))
+		return 1;
+	return 2;
+}
+
+/* Returns 1 if input is the bit representation of 0, infinity or nan.  */
+static inline int zeroinfnan(uint64_t i)
+{
+	return 2 * i - 1 >= 2 * asuint64(INFINITY) - 1;
+}
+
+double pow(double x, double y)
+{
+	uint32_t sign_bias = 0;
+	uint64_t ix, iy;
+	uint32_t topx, topy;
+
+	ix = asuint64(x);
+	iy = asuint64(y);
+	topx = top12(x);
+	topy = top12(y);
+	if (predict_false(topx - 0x001 >= 0x7ff - 0x001 ||
+			  (topy & 0x7ff) - 0x3be >= 0x43e - 0x3be)) {
+		/* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
+		   and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1.  */
+		/* Special cases: (x < 0x1p-126 or inf or nan) or
+		   (|y| < 0x1p-65 or |y| >= 0x1p63 or nan).  */
+		if (predict_false(zeroinfnan(iy))) {
+			if (2 * iy == 0)
+				return issignaling_inline(x) ? x + y : 1.0;
+			if (ix == asuint64(1.0))
+				return issignaling_inline(y) ? x + y : 1.0;
+			if (2 * ix > 2 * asuint64(INFINITY) ||
+			    2 * iy > 2 * asuint64(INFINITY))
+				return x + y;
+			if (2 * ix == 2 * asuint64(1.0))
+				return 1.0;
+			if ((2 * ix < 2 * asuint64(1.0)) == !(iy >> 63))
+				return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf.  */
+			return y * y;
 		}
-		/* over/underflow if x is not close to one */
-		if (ix < 0x3fefffff)
-			return hy < 0 ? s*huge*huge : s*tiny*tiny;
-		if (ix > 0x3ff00000)
-			return hy > 0 ? s*huge*huge : s*tiny*tiny;
-		/* now |1-x| is tiny <= 2**-20, suffice to compute
-		   log(x) by x-x^2/2+x^3/3-x^4/4 */
-		t = ax - 1.0;       /* t has 20 trailing zeros */
-		w = (t*t)*(0.5 - t*(0.3333333333333333333333-t*0.25));
-		u = ivln2_h*t;      /* ivln2_h has 21 sig. bits */
-		v = t*ivln2_l - w*ivln2;
-		t1 = u + v;
-		SET_LOW_WORD(t1, 0);
-		t2 = v - (t1-u);
-	} else {
-		double ss,s2,s_h,s_l,t_h,t_l;
-		n = 0;
-		/* take care subnormal number */
-		if (ix < 0x00100000) {
-			ax *= two53;
-			n -= 53;
-			GET_HIGH_WORD(ix,ax);
+		if (predict_false(zeroinfnan(ix))) {
+			double_t x2 = x * x;
+			if (ix >> 63 && checkint(iy) == 1)
+				x2 = -x2;
+			/* Without the barrier some versions of clang hoist the 1/x2 and
+			   thus division by zero exception can be signaled spuriously.  */
+			return iy >> 63 ? fp_barrier(1 / x2) : x2;
 		}
-		n += ((ix)>>20) - 0x3ff;
-		j = ix & 0x000fffff;
-		/* determine interval */
-		ix = j | 0x3ff00000;   /* normalize ix */
-		if (j <= 0x3988E)      /* |x|<sqrt(3/2) */
-			k = 0;
-		else if (j < 0xBB67A)  /* |x|<sqrt(3)   */
-			k = 1;
-		else {
-			k = 0;
-			n += 1;
-			ix -= 0x00100000;
+		/* Here x and y are non-zero finite.  */
+		if (ix >> 63) {
+			/* Finite x < 0.  */
+			int yint = checkint(iy);
+			if (yint == 0)
+				return __math_invalid(x);
+			if (yint == 1)
+				sign_bias = SIGN_BIAS;
+			ix &= 0x7fffffffffffffff;
+			topx &= 0x7ff;
+		}
+		if ((topy & 0x7ff) - 0x3be >= 0x43e - 0x3be) {
+			/* Note: sign_bias == 0 here because y is not odd.  */
+			if (ix == asuint64(1.0))
+				return 1.0;
+			if ((topy & 0x7ff) < 0x3be) {
+				/* |y| < 2^-65, x^y ~= 1 + y*log(x).  */
+				if (WANT_ROUNDING)
+					return ix > asuint64(1.0) ? 1.0 + y :
+								    1.0 - y;
+				else
+					return 1.0;
+			}
+			return (ix > asuint64(1.0)) == (topy < 0x800) ?
+				       __math_oflow(0) :
+				       __math_uflow(0);
+		}
+		if (topx == 0) {
+			/* Normalize subnormal x so exponent becomes negative.  */
+			ix = asuint64(x * 0x1p52);
+			ix &= 0x7fffffffffffffff;
+			ix -= 52ULL << 52;
 		}
-		SET_HIGH_WORD(ax, ix);
-
-		/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
-		u = ax - bp[k];        /* bp[0]=1.0, bp[1]=1.5 */
-		v = 1.0/(ax+bp[k]);
-		ss = u*v;
-		s_h = ss;
-		SET_LOW_WORD(s_h, 0);
-		/* t_h=ax+bp[k] High */
-		t_h = 0.0;
-		SET_HIGH_WORD(t_h, ((ix>>1)|0x20000000) + 0x00080000 + (k<<18));
-		t_l = ax - (t_h-bp[k]);
-		s_l = v*((u-s_h*t_h)-s_h*t_l);
-		/* compute log(ax) */
-		s2 = ss*ss;
-		r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
-		r += s_l*(s_h+ss);
-		s2 = s_h*s_h;
-		t_h = 3.0 + s2 + r;
-		SET_LOW_WORD(t_h, 0);
-		t_l = r - ((t_h-3.0)-s2);
-		/* u+v = ss*(1+...) */
-		u = s_h*t_h;
-		v = s_l*t_h + t_l*ss;
-		/* 2/(3log2)*(ss+...) */
-		p_h = u + v;
-		SET_LOW_WORD(p_h, 0);
-		p_l = v - (p_h-u);
-		z_h = cp_h*p_h;        /* cp_h+cp_l = 2/(3*log2) */
-		z_l = cp_l*p_h+p_l*cp + dp_l[k];
-		/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
-		t = (double)n;
-		t1 = ((z_h + z_l) + dp_h[k]) + t;
-		SET_LOW_WORD(t1, 0);
-		t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
 	}
 
-	/* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
-	y1 = y;
-	SET_LOW_WORD(y1, 0);
-	p_l = (y-y1)*t1 + y*t2;
-	p_h = y1*t1;
-	z = p_l + p_h;
-	EXTRACT_WORDS(j, i, z);
-	if (j >= 0x40900000) {                      /* z >= 1024 */
-		if (((j-0x40900000)|i) != 0)        /* if z > 1024 */
-			return s*huge*huge;         /* overflow */
-		if (p_l + ovt > z - p_h)
-			return s*huge*huge;         /* overflow */
-	} else if ((j&0x7fffffff) >= 0x4090cc00) {  /* z <= -1075 */  // FIXME: instead of abs(j) use unsigned j
-		if (((j-0xc090cc00)|i) != 0)        /* z < -1075 */
-			return s*tiny*tiny;         /* underflow */
-		if (p_l <= z - p_h)
-			return s*tiny*tiny;         /* underflow */
-	}
-	/*
-	 * compute 2**(p_h+p_l)
-	 */
-	i = j & 0x7fffffff;
-	k = (i>>20) - 0x3ff;
-	n = 0;
-	if (i > 0x3fe00000) {  /* if |z| > 0.5, set n = [z+0.5] */
-		n = j + (0x00100000>>(k+1));
-		k = ((n&0x7fffffff)>>20) - 0x3ff;  /* new k for n */
-		t = 0.0;
-		SET_HIGH_WORD(t, n & ~(0x000fffff>>k));
-		n = ((n&0x000fffff)|0x00100000)>>(20-k);
-		if (j < 0)
-			n = -n;
-		p_h -= t;
-	}
-	t = p_l + p_h;
-	SET_LOW_WORD(t, 0);
-	u = t*lg2_h;
-	v = (p_l-(t-p_h))*lg2 + t*lg2_l;
-	z = u + v;
-	w = v - (z-u);
-	t = z*z;
-	t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
-	r = (z*t1)/(t1-2.0) - (w + z*w);
-	z = 1.0 - (r-z);
-	GET_HIGH_WORD(j, z);
-	j += n<<20;
-	if ((j>>20) <= 0)  /* subnormal output */
-		z = scalbn(z,n);
-	else
-		SET_HIGH_WORD(z, j);
-	return s*z;
+	double_t lo;
+	double_t hi = log_inline(ix, &lo);
+	double_t ehi, elo;
+#if __FP_FAST_FMA
+	ehi = y * hi;
+	elo = y * lo + __builtin_fma(y, hi, -ehi);
+#else
+	double_t yhi = asdouble(iy & -1ULL << 27);
+	double_t ylo = y - yhi;
+	double_t lhi = asdouble(asuint64(hi) & -1ULL << 27);
+	double_t llo = hi - lhi + lo;
+	ehi = yhi * lhi;
+	elo = ylo * lhi + y * llo; /* |elo| < |ehi| * 2^-25.  */
+#endif
+	return exp_inline(ehi, elo, sign_bias);
 }
diff --git a/src/math/pow_data.c b/src/math/pow_data.c
new file mode 100644
index 00000000..81e760de
--- /dev/null
+++ b/src/math/pow_data.c
@@ -0,0 +1,180 @@
+/*
+ * Data for the log part of pow.
+ *
+ * Copyright (c) 2018, Arm Limited.
+ * SPDX-License-Identifier: MIT
+ */
+
+#include "pow_data.h"
+
+#define N (1 << POW_LOG_TABLE_BITS)
+
+const struct pow_log_data __pow_log_data = {
+.ln2hi = 0x1.62e42fefa3800p-1,
+.ln2lo = 0x1.ef35793c76730p-45,
+.poly = {
+// relative error: 0x1.11922ap-70
+// in -0x1.6bp-8 0x1.6bp-8
+// Coefficients are scaled to match the scaling during evaluation.
+-0x1p-1,
+0x1.555555555556p-2 * -2,
+-0x1.0000000000006p-2 * -2,
+0x1.999999959554ep-3 * 4,
+-0x1.555555529a47ap-3 * 4,
+0x1.2495b9b4845e9p-3 * -8,
+-0x1.0002b8b263fc3p-3 * -8,
+},
+/* Algorithm:
+
+	x = 2^k z
+	log(x) = k ln2 + log(c) + log(z/c)
+	log(z/c) = poly(z/c - 1)
+
+where z is in [0x1.69555p-1; 0x1.69555p0] which is split into N subintervals
+and z falls into the ith one, then table entries are computed as
+
+	tab[i].invc = 1/c
+	tab[i].logc = round(0x1p43*log(c))/0x1p43
+	tab[i].logctail = (double)(log(c) - logc)
+
+where c is chosen near the center of the subinterval such that 1/c has only a
+few precision bits so z/c - 1 is exactly representible as double:
+
+	1/c = center < 1 ? round(N/center)/N : round(2*N/center)/N/2
+
+Note: |z/c - 1| < 1/N for the chosen c, |log(c) - logc - logctail| < 0x1p-97,
+the last few bits of logc are rounded away so k*ln2hi + logc has no rounding
+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 = {
+#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)
+A(0x1.6400000000000p+0, -0x1.51aad872df800p-2, -0x1.684e49eb067d5p-49)
+A(0x1.6200000000000p+0, -0x1.4be5f95777800p-2, -0x1.41b6993293ee0p-47)
+A(0x1.6000000000000p+0, -0x1.4618bc21c6000p-2, 0x1.3d82f484c84ccp-46)
+A(0x1.5e00000000000p+0, -0x1.404308686a800p-2, 0x1.c42f3ed820b3ap-50)
+A(0x1.5c00000000000p+0, -0x1.3a64c55694800p-2, 0x1.0b1c686519460p-45)
+A(0x1.5a00000000000p+0, -0x1.347dd9a988000p-2, 0x1.5594dd4c58092p-45)
+A(0x1.5800000000000p+0, -0x1.2e8e2bae12000p-2, 0x1.67b1e99b72bd8p-45)
+A(0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46)
+A(0x1.5600000000000p+0, -0x1.2895a13de8800p-2, 0x1.5ca14b6cfb03fp-46)
+A(0x1.5400000000000p+0, -0x1.22941fbcf7800p-2, -0x1.65a242853da76p-46)
+A(0x1.5200000000000p+0, -0x1.1c898c1699800p-2, -0x1.fafbc68e75404p-46)
+A(0x1.5000000000000p+0, -0x1.1675cababa800p-2, 0x1.f1fc63382a8f0p-46)
+A(0x1.4e00000000000p+0, -0x1.1058bf9ae4800p-2, -0x1.6a8c4fd055a66p-45)
+A(0x1.4c00000000000p+0, -0x1.0a324e2739000p-2, -0x1.c6bee7ef4030ep-47)
+A(0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48)
+A(0x1.4a00000000000p+0, -0x1.0402594b4d000p-2, -0x1.036b89ef42d7fp-48)
+A(0x1.4800000000000p+0, -0x1.fb9186d5e4000p-3, 0x1.d572aab993c87p-47)
+A(0x1.4600000000000p+0, -0x1.ef0adcbdc6000p-3, 0x1.b26b79c86af24p-45)
+A(0x1.4400000000000p+0, -0x1.e27076e2af000p-3, -0x1.72f4f543fff10p-46)
+A(0x1.4200000000000p+0, -0x1.d5c216b4fc000p-3, 0x1.1ba91bbca681bp-45)
+A(0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45)
+A(0x1.4000000000000p+0, -0x1.c8ff7c79aa000p-3, 0x1.7794f689f8434p-45)
+A(0x1.3e00000000000p+0, -0x1.bc286742d9000p-3, 0x1.94eb0318bb78fp-46)
+A(0x1.3c00000000000p+0, -0x1.af3c94e80c000p-3, 0x1.a4e633fcd9066p-52)
+A(0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45)
+A(0x1.3a00000000000p+0, -0x1.a23bc1fe2b000p-3, -0x1.58c64dc46c1eap-45)
+A(0x1.3800000000000p+0, -0x1.9525a9cf45000p-3, -0x1.ad1d904c1d4e3p-45)
+A(0x1.3600000000000p+0, -0x1.87fa06520d000p-3, 0x1.bbdbf7fdbfa09p-45)
+A(0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45)
+A(0x1.3400000000000p+0, -0x1.7ab890210e000p-3, 0x1.bdb9072534a58p-45)
+A(0x1.3200000000000p+0, -0x1.6d60fe719d000p-3, -0x1.0e46aa3b2e266p-46)
+A(0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46)
+A(0x1.3000000000000p+0, -0x1.5ff3070a79000p-3, -0x1.e9e439f105039p-46)
+A(0x1.2e00000000000p+0, -0x1.526e5e3a1b000p-3, -0x1.0de8b90075b8fp-45)
+A(0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46)
+A(0x1.2c00000000000p+0, -0x1.44d2b6ccb8000p-3, 0x1.70cc16135783cp-46)
+A(0x1.2a00000000000p+0, -0x1.371fc201e9000p-3, 0x1.178864d27543ap-48)
+A(0x1.2800000000000p+0, -0x1.29552f81ff000p-3, -0x1.48d301771c408p-45)
+A(0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45)
+A(0x1.2600000000000p+0, -0x1.1b72ad52f6000p-3, -0x1.e80a41811a396p-45)
+A(0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47)
+A(0x1.2400000000000p+0, -0x1.0d77e7cd09000p-3, 0x1.a699688e85bf4p-47)
+A(0x1.2200000000000p+0, -0x1.fec9131dbe000p-4, -0x1.575545ca333f2p-45)
+A(0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45)
+A(0x1.2000000000000p+0, -0x1.e27076e2b0000p-4, 0x1.a342c2af0003cp-45)
+A(0x1.1e00000000000p+0, -0x1.c5e548f5bc000p-4, -0x1.d0c57585fbe06p-46)
+A(0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45)
+A(0x1.1c00000000000p+0, -0x1.a926d3a4ae000p-4, 0x1.53935e85baac8p-45)
+A(0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46)
+A(0x1.1a00000000000p+0, -0x1.8c345d631a000p-4, 0x1.37c294d2f5668p-46)
+A(0x1.1800000000000p+0, -0x1.6f0d28ae56000p-4, -0x1.69737c93373dap-45)
+A(0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46)
+A(0x1.1600000000000p+0, -0x1.51b073f062000p-4, 0x1.f025b61c65e57p-46)
+A(0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45)
+A(0x1.1400000000000p+0, -0x1.341d7961be000p-4, 0x1.c5edaccf913dfp-45)
+A(0x1.1200000000000p+0, -0x1.16536eea38000p-4, 0x1.47c5e768fa309p-46)
+A(0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45)
+A(0x1.1000000000000p+0, -0x1.f0a30c0118000p-5, 0x1.d599e83368e91p-45)
+A(0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46)
+A(0x1.0e00000000000p+0, -0x1.b42dd71198000p-5, 0x1.c827ae5d6704cp-46)
+A(0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45)
+A(0x1.0c00000000000p+0, -0x1.77458f632c000p-5, -0x1.cfc4634f2a1eep-45)
+A(0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48)
+A(0x1.0a00000000000p+0, -0x1.39e87b9fec000p-5, 0x1.502b7f526feaap-48)
+A(0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45)
+A(0x1.0800000000000p+0, -0x1.f829b0e780000p-6, -0x1.980267c7e09e4p-45)
+A(0x1.0600000000000p+0, -0x1.7b91b07d58000p-6, -0x1.88d5493faa639p-45)
+A(0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50)
+A(0x1.0400000000000p+0, -0x1.fc0a8b0fc0000p-7, -0x1.f1e7cf6d3a69cp-50)
+A(0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46)
+A(0x1.0200000000000p+0, -0x1.fe02a6b100000p-8, -0x1.9e23f0dda40e4p-46)
+A(0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0)
+A(0x1.0000000000000p+0, 0x0.0000000000000p+0, 0x0.0000000000000p+0)
+A(0x1.fc00000000000p-1, 0x1.0101575890000p-7, -0x1.0c76b999d2be8p-46)
+A(0x1.f800000000000p-1, 0x1.0205658938000p-6, -0x1.3dc5b06e2f7d2p-45)
+A(0x1.f400000000000p-1, 0x1.8492528c90000p-6, -0x1.aa0ba325a0c34p-45)
+A(0x1.f000000000000p-1, 0x1.0415d89e74000p-5, 0x1.111c05cf1d753p-47)
+A(0x1.ec00000000000p-1, 0x1.466aed42e0000p-5, -0x1.c167375bdfd28p-45)
+A(0x1.e800000000000p-1, 0x1.894aa149fc000p-5, -0x1.97995d05a267dp-46)
+A(0x1.e400000000000p-1, 0x1.ccb73cdddc000p-5, -0x1.a68f247d82807p-46)
+A(0x1.e200000000000p-1, 0x1.eea31c006c000p-5, -0x1.e113e4fc93b7bp-47)
+A(0x1.de00000000000p-1, 0x1.1973bd1466000p-4, -0x1.5325d560d9e9bp-45)
+A(0x1.da00000000000p-1, 0x1.3bdf5a7d1e000p-4, 0x1.cc85ea5db4ed7p-45)
+A(0x1.d600000000000p-1, 0x1.5e95a4d97a000p-4, -0x1.c69063c5d1d1ep-45)
+A(0x1.d400000000000p-1, 0x1.700d30aeac000p-4, 0x1.c1e8da99ded32p-49)
+A(0x1.d000000000000p-1, 0x1.9335e5d594000p-4, 0x1.3115c3abd47dap-45)
+A(0x1.cc00000000000p-1, 0x1.b6ac88dad6000p-4, -0x1.390802bf768e5p-46)
+A(0x1.ca00000000000p-1, 0x1.c885801bc4000p-4, 0x1.646d1c65aacd3p-45)
+A(0x1.c600000000000p-1, 0x1.ec739830a2000p-4, -0x1.dc068afe645e0p-45)
+A(0x1.c400000000000p-1, 0x1.fe89139dbe000p-4, -0x1.534d64fa10afdp-45)
+A(0x1.c000000000000p-1, 0x1.1178e8227e000p-3, 0x1.1ef78ce2d07f2p-45)
+A(0x1.be00000000000p-1, 0x1.1aa2b7e23f000p-3, 0x1.ca78e44389934p-45)
+A(0x1.ba00000000000p-1, 0x1.2d1610c868000p-3, 0x1.39d6ccb81b4a1p-47)
+A(0x1.b800000000000p-1, 0x1.365fcb0159000p-3, 0x1.62fa8234b7289p-51)
+A(0x1.b400000000000p-1, 0x1.4913d8333b000p-3, 0x1.5837954fdb678p-45)
+A(0x1.b200000000000p-1, 0x1.527e5e4a1b000p-3, 0x1.633e8e5697dc7p-45)
+A(0x1.ae00000000000p-1, 0x1.6574ebe8c1000p-3, 0x1.9cf8b2c3c2e78p-46)
+A(0x1.ac00000000000p-1, 0x1.6f0128b757000p-3, -0x1.5118de59c21e1p-45)
+A(0x1.aa00000000000p-1, 0x1.7898d85445000p-3, -0x1.c661070914305p-46)
+A(0x1.a600000000000p-1, 0x1.8beafeb390000p-3, -0x1.73d54aae92cd1p-47)
+A(0x1.a400000000000p-1, 0x1.95a5adcf70000p-3, 0x1.7f22858a0ff6fp-47)
+A(0x1.a000000000000p-1, 0x1.a93ed3c8ae000p-3, -0x1.8724350562169p-45)
+A(0x1.9e00000000000p-1, 0x1.b31d8575bd000p-3, -0x1.c358d4eace1aap-47)
+A(0x1.9c00000000000p-1, 0x1.bd087383be000p-3, -0x1.d4bc4595412b6p-45)
+A(0x1.9a00000000000p-1, 0x1.c6ffbc6f01000p-3, -0x1.1ec72c5962bd2p-48)
+A(0x1.9600000000000p-1, 0x1.db13db0d49000p-3, -0x1.aff2af715b035p-45)
+A(0x1.9400000000000p-1, 0x1.e530effe71000p-3, 0x1.212276041f430p-51)
+A(0x1.9200000000000p-1, 0x1.ef5ade4dd0000p-3, -0x1.a211565bb8e11p-51)
+A(0x1.9000000000000p-1, 0x1.f991c6cb3b000p-3, 0x1.bcbecca0cdf30p-46)
+A(0x1.8c00000000000p-1, 0x1.07138604d5800p-2, 0x1.89cdb16ed4e91p-48)
+A(0x1.8a00000000000p-1, 0x1.0c42d67616000p-2, 0x1.7188b163ceae9p-45)
+A(0x1.8800000000000p-1, 0x1.1178e8227e800p-2, -0x1.c210e63a5f01cp-45)
+A(0x1.8600000000000p-1, 0x1.16b5ccbacf800p-2, 0x1.b9acdf7a51681p-45)
+A(0x1.8400000000000p-1, 0x1.1bf99635a6800p-2, 0x1.ca6ed5147bdb7p-45)
+A(0x1.8200000000000p-1, 0x1.214456d0eb800p-2, 0x1.a87deba46baeap-47)
+A(0x1.7e00000000000p-1, 0x1.2bef07cdc9000p-2, 0x1.a9cfa4a5004f4p-45)
+A(0x1.7c00000000000p-1, 0x1.314f1e1d36000p-2, -0x1.8e27ad3213cb8p-45)
+A(0x1.7a00000000000p-1, 0x1.36b6776be1000p-2, 0x1.16ecdb0f177c8p-46)
+A(0x1.7800000000000p-1, 0x1.3c25277333000p-2, 0x1.83b54b606bd5cp-46)
+A(0x1.7600000000000p-1, 0x1.419b423d5e800p-2, 0x1.8e436ec90e09dp-47)
+A(0x1.7400000000000p-1, 0x1.4718dc271c800p-2, -0x1.f27ce0967d675p-45)
+A(0x1.7200000000000p-1, 0x1.4c9e09e173000p-2, -0x1.e20891b0ad8a4p-45)
+A(0x1.7000000000000p-1, 0x1.522ae0738a000p-2, 0x1.ebe708164c759p-45)
+A(0x1.6e00000000000p-1, 0x1.57bf753c8d000p-2, 0x1.fadedee5d40efp-46)
+A(0x1.6c00000000000p-1, 0x1.5d5bddf596000p-2, -0x1.a0b2a08a465dcp-47)
+},
+};
diff --git a/src/math/pow_data.h b/src/math/pow_data.h
new file mode 100644
index 00000000..5d609ae8
--- /dev/null
+++ b/src/math/pow_data.h
@@ -0,0 +1,22 @@
+/*
+ * Copyright (c) 2018, Arm Limited.
+ * SPDX-License-Identifier: MIT
+ */
+#ifndef _POW_DATA_H
+#define _POW_DATA_H
+
+#include <features.h>
+
+#define POW_LOG_TABLE_BITS 7
+#define POW_LOG_POLY_ORDER 8
+extern hidden const struct pow_log_data {
+	double ln2hi;
+	double ln2lo;
+	double poly[POW_LOG_POLY_ORDER - 1]; /* First coefficient is 1.  */
+	/* Note: the pad field is unused, but allows slightly faster indexing.  */
+	struct {
+		double invc, pad, logc, logctail;
+	} tab[1 << POW_LOG_TABLE_BITS];
+} __pow_log_data;
+
+#endif