1/* 2 * Double-precision x^y function. 3 * 4 * Copyright (c) 2018-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8#include "math_config.h" 9 10/* Scalar version of pow used for fallbacks in vector implementations. */ 11 12/* Data is defined in v_pow_log_data.c. */ 13#define N_LOG (1 << V_POW_LOG_TABLE_BITS) 14#define Off 0x3fe6955500000000 15#define As __v_pow_log_data.poly 16 17/* Data is defined in v_pow_exp_data.c. */ 18#define N_EXP (1 << V_POW_EXP_TABLE_BITS) 19#define SignBias (0x800 << V_POW_EXP_TABLE_BITS) 20#define SmallExp 0x3c9 /* top12(0x1p-54). */ 21#define BigExp 0x408 /* top12(512.0). */ 22#define ThresExp 0x03f /* BigExp - SmallExp. */ 23#define InvLn2N __v_pow_exp_data.n_over_ln2 24#define Ln2HiN __v_pow_exp_data.ln2_over_n_hi 25#define Ln2LoN __v_pow_exp_data.ln2_over_n_lo 26#define SBits __v_pow_exp_data.sbits 27#define Cs __v_pow_exp_data.poly 28 29/* Constants associated with pow. */ 30#define SmallPowX 0x001 /* top12(0x1p-126). */ 31#define BigPowX 0x7ff /* top12(INFINITY). */ 32#define ThresPowX 0x7fe /* BigPowX - SmallPowX. */ 33#define SmallPowY 0x3be /* top12(0x1.e7b6p-65). */ 34#define BigPowY 0x43e /* top12(0x1.749p62). */ 35#define ThresPowY 0x080 /* BigPowY - SmallPowY. */ 36 37/* Top 12 bits of a double (sign and exponent bits). */ 38static inline uint32_t 39top12 (double x) 40{ 41 return asuint64 (x) >> 52; 42} 43 44/* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about 45 additional 15 bits precision. IX is the bit representation of x, but 46 normalized in the subnormal range using the sign bit for the exponent. */ 47static inline double 48log_inline (uint64_t ix, double *tail) 49{ 50 /* x = 2^k z; where z is in range [Off,2*Off) and exact. 51 The range is split into N subintervals. 52 The ith subinterval contains z and c is near its center. */ 53 uint64_t tmp = ix - Off; 54 int i = (tmp >> (52 - V_POW_LOG_TABLE_BITS)) & (N_LOG - 1); 55 int k = (int64_t) tmp >> 52; /* arithmetic shift. */ 56 uint64_t iz = ix - (tmp & 0xfffULL << 52); 57 double z = asdouble (iz); 58 double kd = (double) k; 59 60 /* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */ 61 double invc = __v_pow_log_data.invc[i]; 62 double logc = __v_pow_log_data.logc[i]; 63 double logctail = __v_pow_log_data.logctail[i]; 64 65 /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and 66 |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */ 67 double r = fma (z, invc, -1.0); 68 69 /* k*Ln2 + log(c) + r. */ 70 double t1 = kd * __v_pow_log_data.ln2_hi + logc; 71 double t2 = t1 + r; 72 double lo1 = kd * __v_pow_log_data.ln2_lo + logctail; 73 double lo2 = t1 - t2 + r; 74 75 /* Evaluation is optimized assuming superscalar pipelined execution. */ 76 double ar = As[0] * r; 77 double ar2 = r * ar; 78 double ar3 = r * ar2; 79 /* k*Ln2 + log(c) + r + A[0]*r*r. */ 80 double hi = t2 + ar2; 81 double lo3 = fma (ar, r, -ar2); 82 double lo4 = t2 - hi + ar2; 83 /* p = log1p(r) - r - A[0]*r*r. */ 84 double p = (ar3 85 * (As[1] + r * As[2] 86 + ar2 * (As[3] + r * As[4] + ar2 * (As[5] + r * As[6])))); 87 double lo = lo1 + lo2 + lo3 + lo4 + p; 88 double y = hi + lo; 89 *tail = hi - y + lo; 90 return y; 91} 92 93/* Handle cases that may overflow or underflow when computing the result that 94 is scale*(1+TMP) without intermediate rounding. The bit representation of 95 scale is in SBITS, however it has a computed exponent that may have 96 overflown into the sign bit so that needs to be adjusted before using it as 97 a double. (int32_t)KI is the k used in the argument reduction and exponent 98 adjustment of scale, positive k here means the result may overflow and 99 negative k means the result may underflow. */ 100static inline double 101special_case (double tmp, uint64_t sbits, uint64_t ki) 102{ 103 double scale, y; 104 105 if ((ki & 0x80000000) == 0) 106 { 107 /* k > 0, the exponent of scale might have overflowed by <= 460. */ 108 sbits -= 1009ull << 52; 109 scale = asdouble (sbits); 110 y = 0x1p1009 * (scale + scale * tmp); 111 return check_oflow (eval_as_double (y)); 112 } 113 /* k < 0, need special care in the subnormal range. */ 114 sbits += 1022ull << 52; 115 /* Note: sbits is signed scale. */ 116 scale = asdouble (sbits); 117 y = scale + scale * tmp; 118#if WANT_SIMD_EXCEPT 119 if (fabs (y) < 1.0) 120 { 121 /* Round y to the right precision before scaling it into the subnormal 122 range to avoid double rounding that can cause 0.5+E/2 ulp error where 123 E is the worst-case ulp error outside the subnormal range. So this 124 is only useful if the goal is better than 1 ulp worst-case error. */ 125 double hi, lo, one = 1.0; 126 if (y < 0.0) 127 one = -1.0; 128 lo = scale - y + scale * tmp; 129 hi = one + y; 130 lo = one - hi + y + lo; 131 y = eval_as_double (hi + lo) - one; 132 /* Fix the sign of 0. */ 133 if (y == 0.0) 134 y = asdouble (sbits & 0x8000000000000000); 135 /* The underflow exception needs to be signaled explicitly. */ 136 force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022); 137 } 138#endif 139 y = 0x1p-1022 * y; 140 return check_uflow (eval_as_double (y)); 141} 142 143/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|. 144 The sign_bias argument is SignBias or 0 and sets the sign to -1 or 1. */ 145static inline double 146exp_inline (double x, double xtail, uint32_t sign_bias) 147{ 148 uint32_t abstop = top12 (x) & 0x7ff; 149 if (unlikely (abstop - SmallExp >= ThresExp)) 150 { 151 if (abstop - SmallExp >= 0x80000000) 152 { 153 /* Avoid spurious underflow for tiny x. */ 154 /* Note: 0 is common input. */ 155 return sign_bias ? -1.0 : 1.0; 156 } 157 if (abstop >= top12 (1024.0)) 158 { 159 /* Note: inf and nan are already handled. */ 160 /* Skip errno handling. */ 161#if WANT_SIMD_EXCEPT 162 return asuint64 (x) >> 63 ? __math_uflow (sign_bias) 163 : __math_oflow (sign_bias); 164#else 165 double res_uoflow = asuint64 (x) >> 63 ? 0.0 : INFINITY; 166 return sign_bias ? -res_uoflow : res_uoflow; 167#endif 168 } 169 /* Large x is special cased below. */ 170 abstop = 0; 171 } 172 173 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ 174 /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */ 175 double z = InvLn2N * x; 176 double kd = round (z); 177 uint64_t ki = lround (z); 178 double r = x - kd * Ln2HiN - kd * Ln2LoN; 179 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */ 180 r += xtail; 181 /* 2^(k/N) ~= scale. */ 182 uint64_t idx = ki & (N_EXP - 1); 183 uint64_t top = (ki + sign_bias) << (52 - V_POW_EXP_TABLE_BITS); 184 /* This is only a valid scale when -1023*N < k < 1024*N. */ 185 uint64_t sbits = SBits[idx] + top; 186 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */ 187 /* Evaluation is optimized assuming superscalar pipelined execution. */ 188 double r2 = r * r; 189 double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]); 190 if (unlikely (abstop == 0)) 191 return special_case (tmp, sbits, ki); 192 double scale = asdouble (sbits); 193 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there 194 is no spurious underflow here even without fma. */ 195 return eval_as_double (scale + scale * tmp); 196} 197 198/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|. 199 A version of exp_inline that is not inlined and for which sign_bias is 200 equal to 0. */ 201static double NOINLINE 202exp_nosignbias (double x, double xtail) 203{ 204 uint32_t abstop = top12 (x) & 0x7ff; 205 if (unlikely (abstop - SmallExp >= ThresExp)) 206 { 207 /* Avoid spurious underflow for tiny x. */ 208 if (abstop - SmallExp >= 0x80000000) 209 return 1.0; 210 /* Note: inf and nan are already handled. */ 211 if (abstop >= top12 (1024.0)) 212#if WANT_SIMD_EXCEPT 213 return asuint64 (x) >> 63 ? __math_uflow (0) : __math_oflow (0); 214#else 215 return asuint64 (x) >> 63 ? 0.0 : INFINITY; 216#endif 217 /* Large x is special cased below. */ 218 abstop = 0; 219 } 220 221 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ 222 /* x = ln2/N*k + r, with k integer and r in [-ln2/2N, ln2/2N]. */ 223 double z = InvLn2N * x; 224 double kd = round (z); 225 uint64_t ki = lround (z); 226 double r = x - kd * Ln2HiN - kd * Ln2LoN; 227 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */ 228 r += xtail; 229 /* 2^(k/N) ~= scale. */ 230 uint64_t idx = ki & (N_EXP - 1); 231 uint64_t top = ki << (52 - V_POW_EXP_TABLE_BITS); 232 /* This is only a valid scale when -1023*N < k < 1024*N. */ 233 uint64_t sbits = SBits[idx] + top; 234 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */ 235 double r2 = r * r; 236 double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]); 237 if (unlikely (abstop == 0)) 238 return special_case (tmp, sbits, ki); 239 double scale = asdouble (sbits); 240 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there 241 is no spurious underflow here even without fma. */ 242 return eval_as_double (scale + scale * tmp); 243} 244 245/* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is 246 the bit representation of a non-zero finite floating-point value. */ 247static inline int 248checkint (uint64_t iy) 249{ 250 int e = iy >> 52 & 0x7ff; 251 if (e < 0x3ff) 252 return 0; 253 if (e > 0x3ff + 52) 254 return 2; 255 if (iy & ((1ULL << (0x3ff + 52 - e)) - 1)) 256 return 0; 257 if (iy & (1ULL << (0x3ff + 52 - e))) 258 return 1; 259 return 2; 260} 261 262/* Returns 1 if input is the bit representation of 0, infinity or nan. */ 263static inline int 264zeroinfnan (uint64_t i) 265{ 266 return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1; 267} 268 269static double NOINLINE 270__pl_finite_pow (double x, double y) 271{ 272 uint32_t sign_bias = 0; 273 uint64_t ix, iy; 274 uint32_t topx, topy; 275 276 ix = asuint64 (x); 277 iy = asuint64 (y); 278 topx = top12 (x); 279 topy = top12 (y); 280 if (unlikely (topx - SmallPowX >= ThresPowX 281 || (topy & 0x7ff) - SmallPowY >= ThresPowY)) 282 { 283 /* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0 284 and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */ 285 /* Special cases: (x < 0x1p-126 or inf or nan) or 286 (|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */ 287 if (unlikely (zeroinfnan (iy))) 288 { 289 if (2 * iy == 0) 290 return issignaling_inline (x) ? x + y : 1.0; 291 if (ix == asuint64 (1.0)) 292 return issignaling_inline (y) ? x + y : 1.0; 293 if (2 * ix > 2 * asuint64 (INFINITY) 294 || 2 * iy > 2 * asuint64 (INFINITY)) 295 return x + y; 296 if (2 * ix == 2 * asuint64 (1.0)) 297 return 1.0; 298 if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63)) 299 return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ 300 return y * y; 301 } 302 if (unlikely (zeroinfnan (ix))) 303 { 304 double x2 = x * x; 305 if (ix >> 63 && checkint (iy) == 1) 306 { 307 x2 = -x2; 308 sign_bias = 1; 309 } 310#if WANT_SIMD_EXCEPT 311 if (2 * ix == 0 && iy >> 63) 312 return __math_divzero (sign_bias); 313#endif 314 /* Without the barrier some versions of clang hoist the 1/x2 and 315 thus division by zero exception can be signaled spuriously. */ 316 return iy >> 63 ? opt_barrier_double (1 / x2) : x2; 317 } 318 /* Here x and y are non-zero finite. */ 319 if (ix >> 63) 320 { 321 /* Finite x < 0. */ 322 int yint = checkint (iy); 323 if (yint == 0) 324#if WANT_SIMD_EXCEPT 325 return __math_invalid (x); 326#else 327 return __builtin_nan (""); 328#endif 329 if (yint == 1) 330 sign_bias = SignBias; 331 ix &= 0x7fffffffffffffff; 332 topx &= 0x7ff; 333 } 334 if ((topy & 0x7ff) - SmallPowY >= ThresPowY) 335 { 336 /* Note: sign_bias == 0 here because y is not odd. */ 337 if (ix == asuint64 (1.0)) 338 return 1.0; 339 /* |y| < 2^-65, x^y ~= 1 + y*log(x). */ 340 if ((topy & 0x7ff) < SmallPowY) 341 return 1.0; 342#if WANT_SIMD_EXCEPT 343 return (ix > asuint64 (1.0)) == (topy < 0x800) ? __math_oflow (0) 344 : __math_uflow (0); 345#else 346 return (ix > asuint64 (1.0)) == (topy < 0x800) ? INFINITY : 0; 347#endif 348 } 349 if (topx == 0) 350 { 351 /* Normalize subnormal x so exponent becomes negative. */ 352 /* Without the barrier some versions of clang evalutate the mul 353 unconditionally causing spurious overflow exceptions. */ 354 ix = asuint64 (opt_barrier_double (x) * 0x1p52); 355 ix &= 0x7fffffffffffffff; 356 ix -= 52ULL << 52; 357 } 358 } 359 360 double lo; 361 double hi = log_inline (ix, &lo); 362 double ehi = y * hi; 363 double elo = y * lo + fma (y, hi, -ehi); 364 return exp_inline (ehi, elo, sign_bias); 365} 366