1*0928368fSKristof Beyls /*
2*0928368fSKristof Beyls * Single-precision e^x function.
3*0928368fSKristof Beyls *
4*0928368fSKristof Beyls * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
5*0928368fSKristof Beyls * See https://llvm.org/LICENSE.txt for license information.
6*0928368fSKristof Beyls * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7*0928368fSKristof Beyls */
8*0928368fSKristof Beyls
9*0928368fSKristof Beyls #include <math.h>
10*0928368fSKristof Beyls #include <stdint.h>
11*0928368fSKristof Beyls #include "math_config.h"
12*0928368fSKristof Beyls
13*0928368fSKristof Beyls /*
14*0928368fSKristof Beyls EXP2F_TABLE_BITS = 5
15*0928368fSKristof Beyls EXP2F_POLY_ORDER = 3
16*0928368fSKristof Beyls
17*0928368fSKristof Beyls ULP error: 0.502 (nearest rounding.)
18*0928368fSKristof Beyls Relative error: 1.69 * 2^-34 in [-ln2/64, ln2/64] (before rounding.)
19*0928368fSKristof Beyls Wrong count: 170635 (all nearest rounding wrong results with fma.)
20*0928368fSKristof Beyls Non-nearest ULP error: 1 (rounded ULP error)
21*0928368fSKristof Beyls */
22*0928368fSKristof Beyls
23*0928368fSKristof Beyls #define N (1 << EXP2F_TABLE_BITS)
24*0928368fSKristof Beyls #define InvLn2N __exp2f_data.invln2_scaled
25*0928368fSKristof Beyls #define T __exp2f_data.tab
26*0928368fSKristof Beyls #define C __exp2f_data.poly_scaled
27*0928368fSKristof Beyls
28*0928368fSKristof Beyls static inline uint32_t
top12(float x)29*0928368fSKristof Beyls top12 (float x)
30*0928368fSKristof Beyls {
31*0928368fSKristof Beyls return asuint (x) >> 20;
32*0928368fSKristof Beyls }
33*0928368fSKristof Beyls
34*0928368fSKristof Beyls float
expf(float x)35*0928368fSKristof Beyls expf (float x)
36*0928368fSKristof Beyls {
37*0928368fSKristof Beyls uint32_t abstop;
38*0928368fSKristof Beyls uint64_t ki, t;
39*0928368fSKristof Beyls /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */
40*0928368fSKristof Beyls double_t kd, xd, z, r, r2, y, s;
41*0928368fSKristof Beyls
42*0928368fSKristof Beyls xd = (double_t) x;
43*0928368fSKristof Beyls abstop = top12 (x) & 0x7ff;
44*0928368fSKristof Beyls if (unlikely (abstop >= top12 (88.0f)))
45*0928368fSKristof Beyls {
46*0928368fSKristof Beyls /* |x| >= 88 or x is nan. */
47*0928368fSKristof Beyls if (asuint (x) == asuint (-INFINITY))
48*0928368fSKristof Beyls return 0.0f;
49*0928368fSKristof Beyls if (abstop >= top12 (INFINITY))
50*0928368fSKristof Beyls return x + x;
51*0928368fSKristof Beyls if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */
52*0928368fSKristof Beyls return __math_oflowf (0);
53*0928368fSKristof Beyls if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */
54*0928368fSKristof Beyls return __math_uflowf (0);
55*0928368fSKristof Beyls #if WANT_ERRNO_UFLOW
56*0928368fSKristof Beyls if (x < -0x1.9d1d9ep6f) /* x < log(0x1p-149) ~= -103.28 */
57*0928368fSKristof Beyls return __math_may_uflowf (0);
58*0928368fSKristof Beyls #endif
59*0928368fSKristof Beyls }
60*0928368fSKristof Beyls
61*0928368fSKristof Beyls /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */
62*0928368fSKristof Beyls z = InvLn2N * xd;
63*0928368fSKristof Beyls
64*0928368fSKristof Beyls /* Round and convert z to int, the result is in [-150*N, 128*N] and
65*0928368fSKristof Beyls ideally nearest int is used, otherwise the magnitude of r can be
66*0928368fSKristof Beyls bigger which gives larger approximation error. */
67*0928368fSKristof Beyls #if TOINT_INTRINSICS
68*0928368fSKristof Beyls kd = roundtoint (z);
69*0928368fSKristof Beyls ki = converttoint (z);
70*0928368fSKristof Beyls #else
71*0928368fSKristof Beyls # define SHIFT __exp2f_data.shift
72*0928368fSKristof Beyls kd = eval_as_double (z + SHIFT);
73*0928368fSKristof Beyls ki = asuint64 (kd);
74*0928368fSKristof Beyls kd -= SHIFT;
75*0928368fSKristof Beyls #endif
76*0928368fSKristof Beyls r = z - kd;
77*0928368fSKristof Beyls
78*0928368fSKristof Beyls /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */
79*0928368fSKristof Beyls t = T[ki % N];
80*0928368fSKristof Beyls t += ki << (52 - EXP2F_TABLE_BITS);
81*0928368fSKristof Beyls s = asdouble (t);
82*0928368fSKristof Beyls z = C[0] * r + C[1];
83*0928368fSKristof Beyls r2 = r * r;
84*0928368fSKristof Beyls y = C[2] * r + 1;
85*0928368fSKristof Beyls y = z * r2 + y;
86*0928368fSKristof Beyls y = y * s;
87*0928368fSKristof Beyls return eval_as_float (y);
88*0928368fSKristof Beyls }
89*0928368fSKristof Beyls #if USE_GLIBC_ABI
90*0928368fSKristof Beyls strong_alias (expf, __expf_finite)
91*0928368fSKristof Beyls hidden_alias (expf, __ieee754_expf)
92*0928368fSKristof Beyls #endif
93