2  * Double-precision x^y function.
4  * Copyright (c) 2018-2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
15 Worst-case error: 0.54 ULP (~= ulperr_exp + 1024*Ln2*relerr_log*2^53)
16 relerr_log: 1.3 * 2^-68 (Relative error of log, 1.5 * 2^-68 without fma)
27 /* Top 12 bits of a double (sign and exponent bits).  */
29 top12 (double x)  in top12()
36    normalized in the subnormal range using the sign bit for the exponent.  */
45   /* x = 2^k z; where z is in range [OFF,2*OFF) and exact.  in log_inline()
46      The range is split into N subintervals.  in log_inline()
48   tmp = ix - OFF;  in log_inline()
49   i = (tmp >> (52 - POW_LOG_TABLE_BITS)) % N;  in log_inline()
51   iz = ix - (tmp & 0xfffULL << 52);  in log_inline()
55   /* log(x) = k*Ln2 + log(c) + log1p(z/c-1).  */  in log_inline()
61      |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible.  */  in log_inline()
63   r = fma (z, invc, -1.0);  in log_inline()
66   double_t zhi = asdouble ((iz + (1ULL << 31)) & (-1ULL << 32));  in log_inline()
67   double_t zlo = z - zhi;  in log_inline()
68   double_t rhi = zhi * invc - 1.0;  in log_inline()
77   lo2 = t1 - t2 + r;  in log_inline()
81   ar = A[0] * r; /* A[0] = -0.5.  */  in log_inline()
87   lo3 = fma (ar, r, -ar2);  in log_inline()
88   lo4 = t2 - hi + ar2;  in log_inline()
94   lo4 = t2 - hi + arhi2;  in log_inline()
96   /* p = log1p(r) - r - A[0]*r*r.  */  in log_inline()
103   *tail = hi - y + lo;  in log_inline()
115 #define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
116 #define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
117 #define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
118 #define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
119 #define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
125    a double.  (int32_t)KI is the k used in the argument reduction and exponent
128 static inline double
136       sbits -= 1009ull << 52;  in specialcase()
141   /* k < 0, need special care in the subnormal range.  */  in specialcase()
149 	 range to avoid double rounding that can cause 0.5+E/2 ulp error where  in specialcase()
150 	 E is the worst-case ulp error outside the subnormal range.  So this  in specialcase()
151 	 is only useful if the goal is better than 1 ulp worst-case error.  */  in specialcase()
154 	one = -1.0;  in specialcase()
155       lo = scale - y + scale * tmp;  in specialcase()
157       lo = one - hi + y + lo;  in specialcase()
158       y = eval_as_double (hi + lo) - one;  in specialcase()
163       force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);  in specialcase()
165   y = 0x1p-1022 * y;  in specialcase()
171 /* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
172    The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1.  */
173 static inline double
182   if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))  in exp_inline()
184       if (abstop - top12 (0x1p-54) >= 0x80000000)  in exp_inline()
189 	  return sign_bias ? -one : one;  in exp_inline()
203   /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */  in exp_inline()
204   /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */  in exp_inline()
210   /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes.  */  in exp_inline()
215   /* z - kd is in [-1, 1] in non-nearest rounding modes.  */  in exp_inline()
218   kd -= Shift;  in exp_inline()
221   /* The code assumes 2^-200 < |xtail| < 2^-8/N.  */  in exp_inline()
225   top = (ki + sign_bias) << (52 - EXP_TABLE_BITS);  in exp_inline()
227   /* This is only a valid scale when -1023*N < k < 1024*N.  */  in exp_inline()
229   /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */  in exp_inline()
244   /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there  in exp_inline()
250    the bit representation of a non-zero finite floating-point value.  */
259   if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))  in checkint()
261   if (iy & (1ULL << (0x3ff + 52 - e)))  in checkint()
270   return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1;  in zeroinfnan()
273 double
274 pow (double x, double y)  in pow()
284   if (unlikely (topx - 0x001 >= 0x7ff - 0x001  in pow()
285 		|| (topy & 0x7ff) - 0x3be >= 0x43e - 0x3be))  in pow()
288 	 and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1.  */  in pow()
289       /* Special cases: (x < 0x1p-126 or inf or nan) or  in pow()
290 	 (|y| < 0x1p-65 or |y| >= 0x1p63 or nan).  */  in pow()
303 	    return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf.  */  in pow()
311 	      x2 = -x2;  in pow()
320       /* Here x and y are non-zero finite.  */  in pow()
332       if ((topy & 0x7ff) - 0x3be >= 0x43e - 0x3be)  in pow()
339 	      /* |y| < 2^-65, x^y ~= 1 + y*log(x).  */  in pow()
341 		return ix > asuint64 (1.0) ? 1.0 + y : 1.0 - y;  in pow()
355 	  ix -= 52ULL << 52;  in pow()
364   elo = y * lo + fma (y, hi, -ehi);  in pow()
366   double_t yhi = asdouble (iy & -1ULL << 27);  in pow()
367   double_t ylo = y - yhi;  in pow()
368   double_t lhi = asdouble (asuint64 (hi) & -1ULL << 27);  in pow()
369   double_t llo = hi - lhi + lo;  in pow()
371   elo = ylo * lhi + y * llo; /* |elo| < |ehi| * 2^-25.  */  in pow()
379 long double powl (long double x, long double y) { return pow (x, y); }  in strong_alias()
386 TEST_INTERVAL2 (pow, -0.5, -2.0, 0, inf, 20000)
387 TEST_INTERVAL2 (pow, 0.5, 2.0, -0, -inf, 20000)
388 TEST_INTERVAL2 (pow, -0.5, -2.0, -0, -inf, 20000)
389 TEST_INTERVAL2 (pow, 0.5, 2.0, 0x1p-10, 0x1p10, 40000)
390 TEST_INTERVAL2 (pow, 0.5, 2.0, -0x1p-10, -0x1p10, 40000)
392 TEST_INTERVAL2 (pow, 0, inf, -0.5, -2.0, 80000)
393 TEST_INTERVAL2 (pow, 0x1.fp-1, 0x1.08p0, 0x1p8, 0x1p17, 80000)
394 TEST_INTERVAL2 (pow, 0x1.fp-1, 0x1.08p0, -0x1p8, -0x1p17, 80000)
395 TEST_INTERVAL2 (pow, 0, 0x1p-1000, 0, 1.0, 50000)
397 TEST_INTERVAL2 (pow, 0x1.ffffffffffff0p-1, 0x1.0000000000008p0, 0x1p60, 0x1p68,
399 TEST_INTERVAL2 (pow, 0x1.ffffffffff000p-1, 0x1p0, 0x1p50, 0x1p52, 50000)
400 TEST_INTERVAL2 (pow, -0x1.ffffffffff000p-1, -0x1p0, 0x1p50, 0x1p52, 50000)