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
special_case(uint64_t sbits,double_t tmp,uint64_t ki)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
exp10(double x)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