1/* Copyright (C) 2012-2024 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
3
4 The GNU C Library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Lesser General Public
6 License as published by the Free Software Foundation; either
7 version 2.1 of the License, or (at your option) any later version.
8
9 The GNU C Library is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public
15 License along with the GNU C Library; if not, see
16 <https://www.gnu.org/licenses/>. */
17
18#include <math.h>
19#include <math-barriers.h>
20#include <math-narrow-eval.h>
21#include <math-svid-compat.h>
22#include <libm-alias-finite.h>
23#include <libm-alias-double.h>
24#include "math_config.h"
25
26#define N (1 << EXP_TABLE_BITS)
27#define IndexMask (N - 1)
28#define OFlowBound 0x1.34413509f79ffp8 /* log10(DBL_MAX). */
29#define UFlowBound -0x1.5ep+8 /* -350. */
30#define SmallTop 0x3c6 /* top12(0x1p-57). */
31#define BigTop 0x407 /* top12(0x1p8). */
32#define Thresh 0x41 /* BigTop - SmallTop. */
33#define Shift __exp_data.shift
34#define C(i) __exp_data.exp10_poly[i]
35
36static double
37special_case (uint64_t sbits, double_t tmp, uint64_t ki)
38{
39 double_t scale, y;
40
41 if (ki - (1ull << 16) < 0x80000000)
42 {
43 /* The exponent of scale might have overflowed by 1. */
44 sbits -= 1ull << 52;
45 scale = asdouble (i: sbits);
46 y = 2 * (scale + scale * tmp);
47 return check_oflow (x: y);
48 }
49
50 /* n < 0, need special care in the subnormal range. */
51 sbits += 1022ull << 52;
52 scale = asdouble (i: sbits);
53 y = scale + scale * tmp;
54
55 if (y < 1.0)
56 {
57 /* Round y to the right precision before scaling it into the subnormal
58 range to avoid double rounding that can cause 0.5+E/2 ulp error where
59 E is the worst-case ulp error outside the subnormal range. So this
60 is only useful if the goal is better than 1 ulp worst-case error. */
61 double_t lo = scale - y + scale * tmp;
62 double_t hi = 1.0 + y;
63 lo = 1.0 - hi + y + lo;
64 y = math_narrow_eval (hi + lo) - 1.0;
65 /* Avoid -0.0 with downward rounding. */
66 if (WANT_ROUNDING && y == 0.0)
67 y = 0.0;
68 /* The underflow exception needs to be signaled explicitly. */
69 math_force_eval (math_opt_barrier (0x1p-1022) * 0x1p-1022);
70 }
71 y = 0x1p-1022 * y;
72
73 return check_uflow (x: y);
74}
75
76/* Double-precision 10^x approximation. Largest observed error is ~0.513 ULP. */
77double
78__exp10 (double x)
79{
80 uint64_t ix = asuint64 (f: x);
81 uint32_t abstop = (ix >> 52) & 0x7ff;
82
83 if (__glibc_unlikely (abstop - SmallTop >= Thresh))
84 {
85 if (abstop - SmallTop >= 0x80000000)
86 /* Avoid spurious underflow for tiny x.
87 Note: 0 is common input. */
88 return x + 1;
89 if (abstop == 0x7ff)
90 return ix == asuint64 (f: -INFINITY) ? 0.0 : x + 1.0;
91 if (x >= OFlowBound)
92 return __math_oflow (0);
93 if (x < UFlowBound)
94 return __math_uflow (0);
95
96 /* Large x is special-cased below. */
97 abstop = 0;
98 }
99
100 /* Reduce x: z = x * N / log10(2), k = round(z). */
101 double_t z = __exp_data.invlog10_2N * x;
102 double_t kd;
103 int64_t ki;
104#if TOINT_INTRINSICS
105 kd = roundtoint (z);
106 ki = converttoint (z);
107#else
108 kd = math_narrow_eval (z + Shift);
109 kd -= Shift;
110 ki = kd;
111#endif
112
113 /* r = x - k * log10(2), r in [-0.5, 0.5]. */
114 double_t r = x;
115 r = __exp_data.neglog10_2hiN * kd + r;
116 r = __exp_data.neglog10_2loN * kd + r;
117
118 /* exp10(x) = 2^(k/N) * 2^(r/N).
119 Approximate the two components separately. */
120
121 /* s = 2^(k/N), using lookup table. */
122 uint64_t e = ki << (52 - EXP_TABLE_BITS);
123 uint64_t i = (ki & IndexMask) * 2;
124 uint64_t u = __exp_data.tab[i + 1];
125 uint64_t sbits = u + e;
126
127 double_t tail = asdouble (i: __exp_data.tab[i]);
128
129 /* 2^(r/N) ~= 1 + r * Poly(r). */
130 double_t r2 = r * r;
131 double_t p = C (0) + r * C (1);
132 double_t y = C (2) + r * C (3);
133 y = y + r2 * C (4);
134 y = p + r2 * y;
135 y = tail + y * r;
136
137 if (__glibc_unlikely (abstop == 0))
138 return special_case (sbits, tmp: y, ki);
139
140 /* Assemble components:
141 y = 2^(r/N) * 2^(k/N)
142 ~= (y + 1) * s. */
143 double_t s = asdouble (i: sbits);
144 return s * y + s;
145}
146
147strong_alias (__exp10, __ieee754_exp10)
148libm_alias_finite (__ieee754_exp10, __exp10)
149#if LIBM_SVID_COMPAT
150versioned_symbol (libm, __exp10, exp10, GLIBC_2_39);
151libm_alias_double_other (__exp10, exp10)
152#else
153libm_alias_double (__exp10, exp10)
154#endif
155

source code of glibc/sysdeps/ieee754/dbl-64/e_exp10.c