1 /* Copyright (C) 2012-2024 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
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.
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.
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/>. */
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"
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]
37 special_case (uint64_t sbits
, double_t tmp
, uint64_t ki
)
41 if (ki
- (1ull << 16) < 0x80000000)
43 /* The exponent of scale might have overflowed by 1. */
45 scale
= asdouble (sbits
);
46 y
= 2 * (scale
+ scale
* tmp
);
47 return check_oflow (y
);
50 /* n < 0, need special care in the subnormal range. */
51 sbits
+= 1022ull << 52;
52 scale
= asdouble (sbits
);
53 y
= scale
+ scale
* tmp
;
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)
68 /* The underflow exception needs to be signaled explicitly. */
69 math_force_eval (math_opt_barrier (0x1p
-1022) * 0x1p
-1022);
73 return check_uflow (y
);
76 /* Double-precision 10^x approximation. Largest observed error is ~0.513 ULP. */
80 uint64_t ix
= asuint64 (x
);
81 uint32_t abstop
= (ix
>> 52) & 0x7ff;
83 if (__glibc_unlikely (abstop
- SmallTop
>= Thresh
))
85 if (abstop
- SmallTop
>= 0x80000000)
86 /* Avoid spurious underflow for tiny x.
87 Note: 0 is common input. */
90 return ix
== asuint64 (-INFINITY
) ? 0.0 : x
+ 1.0;
92 return __math_oflow (0);
94 return __math_uflow (0);
96 /* Large x is special-cased below. */
100 /* Reduce x: z = x * N / log10(2), k = round(z). */
101 double_t z
= __exp_data
.invlog10_2N
* x
;
106 ki
= converttoint (z
);
108 kd
= math_narrow_eval (z
+ Shift
);
113 /* r = x - k * log10(2), r in [-0.5, 0.5]. */
115 r
= __exp_data
.neglog10_2hiN
* kd
+ r
;
116 r
= __exp_data
.neglog10_2loN
* kd
+ r
;
118 /* exp10(x) = 2^(k/N) * 2^(r/N).
119 Approximate the two components separately. */
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
;
127 double_t tail
= asdouble (__exp_data
.tab
[i
]);
129 /* 2^(r/N) ~= 1 + r * Poly(r). */
131 double_t p
= C (0) + r
* C (1);
132 double_t y
= C (2) + r
* C (3);
137 if (__glibc_unlikely (abstop
== 0))
138 return special_case (sbits
, y
, ki
);
140 /* Assemble components:
141 y = 2^(r/N) * 2^(k/N)
143 double_t s
= asdouble (sbits
);
147 strong_alias (__exp10
, __ieee754_exp10
)
148 libm_alias_finite (__ieee754_exp10
, __exp10
)
150 versioned_symbol (libm
, __exp10
, exp10
, GLIBC_2_39
);
151 libm_alias_double_other (__exp10
, exp10
)
153 libm_alias_double (__exp10
, exp10
)