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