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