1 /* 2 * Double-precision 10^x function. 3 * 4 * Copyright (c) 2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "math_config.h" 9 10 #define N (1 << EXP_TABLE_BITS) 11 #define IndexMask (N - 1) 12 #define OFlowBound 0x1.34413509f79ffp8 /* log10(DBL_MAX). */ 13 #define UFlowBound -0x1.5ep+8 /* -350. */ 14 #define SmallTop 0x3c6 /* top12(0x1p-57). */ 15 #define BigTop 0x407 /* top12(0x1p8). */ 16 #define Thresh 0x41 /* BigTop - SmallTop. */ 17 #define Shift __exp_data.shift 18 #define C(i) __exp_data.exp10_poly[i] 19 20 static double 21 special_case (uint64_t sbits, double_t tmp, uint64_t ki) 22 { 23 double_t scale, y; 24 25 if (ki - (1ull << 16) < 0x80000000) 26 { 27 /* The exponent of scale might have overflowed by 1. */ 28 sbits -= 1ull << 52; 29 scale = asdouble (sbits); 30 y = 2 * (scale + scale * tmp); 31 return check_oflow (eval_as_double (y)); 32 } 33 34 /* n < 0, need special care in the subnormal range. */ 35 sbits += 1022ull << 52; 36 scale = asdouble (sbits); 37 y = scale + scale * tmp; 38 39 if (y < 1.0) 40 { 41 /* Round y to the right precision before scaling it into the subnormal 42 range to avoid double rounding that can cause 0.5+E/2 ulp error where 43 E is the worst-case ulp error outside the subnormal range. So this 44 is only useful if the goal is better than 1 ulp worst-case error. */ 45 double_t lo = scale - y + scale * tmp; 46 double_t hi = 1.0 + y; 47 lo = 1.0 - hi + y + lo; 48 y = eval_as_double (hi + lo) - 1.0; 49 /* Avoid -0.0 with downward rounding. */ 50 if (WANT_ROUNDING && y == 0.0) 51 y = 0.0; 52 /* The underflow exception needs to be signaled explicitly. */ 53 force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022); 54 } 55 y = 0x1p-1022 * y; 56 57 return check_uflow (y); 58 } 59 60 /* Double-precision 10^x approximation. Largest observed error is ~0.513 ULP. */ 61 double 62 exp10 (double x) 63 { 64 uint64_t ix = asuint64 (x); 65 uint32_t abstop = (ix >> 52) & 0x7ff; 66 67 if (unlikely (abstop - SmallTop >= Thresh)) 68 { 69 if (abstop - SmallTop >= 0x80000000) 70 /* Avoid spurious underflow for tiny x. 71 Note: 0 is common input. */ 72 return x + 1; 73 if (abstop == 0x7ff) 74 return ix == asuint64 (-INFINITY) ? 0.0 : x + 1.0; 75 if (x >= OFlowBound) 76 return __math_oflow (0); 77 if (x < UFlowBound) 78 return __math_uflow (0); 79 80 /* Large x is special-cased below. */ 81 abstop = 0; 82 } 83 84 /* Reduce x: z = x * N / log10(2), k = round(z). */ 85 double_t z = __exp_data.invlog10_2N * x; 86 double_t kd; 87 int64_t ki; 88 #if TOINT_INTRINSICS 89 kd = roundtoint (z); 90 ki = converttoint (z); 91 #else 92 kd = eval_as_double (z + Shift); 93 kd -= Shift; 94 ki = kd; 95 #endif 96 97 /* r = x - k * log10(2), r in [-0.5, 0.5]. */ 98 double_t r = x; 99 r = __exp_data.neglog10_2hiN * kd + r; 100 r = __exp_data.neglog10_2loN * kd + r; 101 102 /* exp10(x) = 2^(k/N) * 2^(r/N). 103 Approximate the two components separately. */ 104 105 /* s = 2^(k/N), using lookup table. */ 106 uint64_t e = ki << (52 - EXP_TABLE_BITS); 107 uint64_t i = (ki & IndexMask) * 2; 108 uint64_t u = __exp_data.tab[i + 1]; 109 uint64_t sbits = u + e; 110 111 double_t tail = asdouble (__exp_data.tab[i]); 112 113 /* 2^(r/N) ~= 1 + r * Poly(r). */ 114 double_t r2 = r * r; 115 double_t p = C (0) + r * C (1); 116 double_t y = C (2) + r * C (3); 117 y = y + r2 * C (4); 118 y = p + r2 * y; 119 y = tail + y * r; 120 121 if (unlikely (abstop == 0)) 122 return special_case (sbits, y, ki); 123 124 /* Assemble components: 125 y = 2^(r/N) * 2^(k/N) 126 ~= (y + 1) * s. */ 127 double_t s = asdouble (sbits); 128 return eval_as_double (s * y + s); 129 } 130