]>
Commit | Line | Data |
---|---|---|
72aa6233 | 1 | /* Single-precision e^x function. |
6d7e8eda | 2 | Copyright (C) 2017-2023 Free Software Foundation, Inc. |
dc30f461 | 3 | This file is part of the GNU C Library. |
f7eac6eb | 4 | |
dc30f461 | 5 | The GNU C Library is free software; you can redistribute it and/or |
41bdb6e2 AJ |
6 | modify it under the terms of the GNU Lesser General Public |
7 | License as published by the Free Software Foundation; either | |
8 | version 2.1 of the License, or (at your option) any later version. | |
f7eac6eb | 9 | |
dc30f461 UD |
10 | The GNU C Library is distributed in the hope that it will be useful, |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
41bdb6e2 | 13 | Lesser General Public License for more details. |
f7eac6eb | 14 | |
41bdb6e2 | 15 | You should have received a copy of the GNU Lesser General Public |
59ba27a6 | 16 | License along with the GNU C Library; if not, see |
5a82c748 | 17 | <https://www.gnu.org/licenses/>. */ |
dc30f461 | 18 | |
f7a0b063 SN |
19 | #ifdef __expf |
20 | # undef libm_hidden_proto | |
21 | # define libm_hidden_proto(ignored) | |
22 | #endif | |
23 | ||
dc30f461 | 24 | #include <math.h> |
aaee3cd8 | 25 | #include <math-narrow-eval.h> |
72aa6233 | 26 | #include <stdint.h> |
220622dd | 27 | #include <libm-alias-finite.h> |
24b6515d | 28 | #include <libm-alias-float.h> |
72aa6233 SN |
29 | #include "math_config.h" |
30 | ||
31 | /* | |
32 | EXP2F_TABLE_BITS = 5 | |
33 | EXP2F_POLY_ORDER = 3 | |
34 | ||
35 | ULP error: 0.502 (nearest rounding.) | |
36 | Relative error: 1.69 * 2^-34 in [-ln2/64, ln2/64] (before rounding.) | |
37 | Wrong count: 170635 (all nearest rounding wrong results with fma.) | |
38 | Non-nearest ULP error: 1 (rounded ULP error) | |
39 | */ | |
40 | ||
41 | #define N (1 << EXP2F_TABLE_BITS) | |
42 | #define InvLn2N __exp2f_data.invln2_scaled | |
43 | #define T __exp2f_data.tab | |
44 | #define C __exp2f_data.poly_scaled | |
45 | ||
46 | static inline uint32_t | |
47 | top12 (float x) | |
48 | { | |
49 | return asuint (x) >> 20; | |
50 | } | |
dc30f461 UD |
51 | |
52 | float | |
f7a0b063 | 53 | __expf (float x) |
f7eac6eb | 54 | { |
72aa6233 SN |
55 | uint32_t abstop; |
56 | uint64_t ki, t; | |
57 | /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ | |
58 | double_t kd, xd, z, r, r2, y, s; | |
59 | ||
60 | xd = (double_t) x; | |
61 | abstop = top12 (x) & 0x7ff; | |
62 | if (__glibc_unlikely (abstop >= top12 (88.0f))) | |
dc30f461 | 63 | { |
72aa6233 SN |
64 | /* |x| >= 88 or x is nan. */ |
65 | if (asuint (x) == asuint (-INFINITY)) | |
66 | return 0.0f; | |
67 | if (abstop >= top12 (INFINITY)) | |
68 | return x + x; | |
69 | if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */ | |
70 | return __math_oflowf (0); | |
71 | if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */ | |
72 | return __math_uflowf (0); | |
73 | #if WANT_ERRNO_UFLOW | |
74 | if (x < -0x1.9d1d9ep6f) /* x < log(0x1p-149) ~= -103.28 */ | |
75 | return __math_may_uflowf (0); | |
76 | #endif | |
dc30f461 | 77 | } |
72aa6233 SN |
78 | |
79 | /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ | |
80 | z = InvLn2N * xd; | |
81 | ||
82 | /* Round and convert z to int, the result is in [-150*N, 128*N] and | |
83 | ideally ties-to-even rule is used, otherwise the magnitude of r | |
84 | can be bigger which gives larger approximation error. */ | |
85 | #if TOINT_INTRINSICS | |
86 | kd = roundtoint (z); | |
87 | ki = converttoint (z); | |
43cfdf8f | 88 | #else |
72aa6233 SN |
89 | # define SHIFT __exp2f_data.shift |
90 | kd = math_narrow_eval ((double) (z + SHIFT)); /* Needs to be double. */ | |
91 | ki = asuint64 (kd); | |
92 | kd -= SHIFT; | |
93 | #endif | |
94 | r = z - kd; | |
95 | ||
96 | /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ | |
97 | t = T[ki % N]; | |
98 | t += ki << (52 - EXP2F_TABLE_BITS); | |
99 | s = asdouble (t); | |
100 | z = C[0] * r + C[1]; | |
101 | r2 = r * r; | |
102 | y = C[2] * r + 1; | |
103 | y = z * r2 + y; | |
104 | y = y * s; | |
105 | return (float) y; | |
f7eac6eb | 106 | } |
f7a0b063 SN |
107 | |
108 | #ifndef __expf | |
109 | hidden_def (__expf) | |
110 | strong_alias (__expf, __ieee754_expf) | |
220622dd | 111 | libm_alias_finite (__ieee754_expf, __expf) |
f7a0b063 | 112 | versioned_symbol (libm, __expf, expf, GLIBC_2_27); |
24b6515d | 113 | libm_alias_float_other (__exp, exp) |
f7a0b063 | 114 | #endif |