1 /*
2 * e_expf.c - single-precision exp function
3 *
4 * Copyright (c) 2009-2018, Arm Limited.
5 * SPDX-License-Identifier: MIT
6 */
7
8 /*
9 * Algorithm was once taken from Cody & Waite, but has been munged
10 * out of all recognition by SGT.
11 */
12
13 #include <math.h>
14 #include <errno.h>
15 #include "math_private.h"
16
17 float
expf(float X)18 expf(float X)
19 {
20 int N; float XN, g, Rg, Result;
21 unsigned ix = fai(X), edgecaseflag = 0;
22
23 /*
24 * Handle infinities, NaNs and big numbers.
25 */
26 if (__builtin_expect((ix << 1) - 0x67000000 > 0x85500000 - 0x67000000, 0)) {
27 if (!(0x7f800000 & ~ix)) {
28 if (ix == 0xff800000)
29 return 0.0f;
30 else
31 return FLOAT_INFNAN(X);/* do the right thing with both kinds of NaN and with +inf */
32 } else if ((ix << 1) < 0x67000000) {
33 return 1.0f; /* magnitude so small the answer can't be distinguished from 1 */
34 } else if ((ix << 1) > 0x85a00000) {
35 __set_errno(ERANGE);
36 if (ix & 0x80000000) {
37 return FLOAT_UNDERFLOW;
38 } else {
39 return FLOAT_OVERFLOW;
40 }
41 } else {
42 edgecaseflag = 1;
43 }
44 }
45
46 /*
47 * Split the input into an integer multiple of log(2)/4, and a
48 * fractional part.
49 */
50 XN = X * (4.0f*1.4426950408889634074f);
51 #ifdef __TARGET_FPU_SOFTVFP
52 XN = _frnd(XN);
53 N = (int)XN;
54 #else
55 N = (int)(XN + (ix & 0x80000000 ? -0.5f : 0.5f));
56 XN = N;
57 #endif
58 g = (X - XN * 0x1.62ep-3F) - XN * 0x1.0bfbe8p-17F; /* prec-and-a-half representation of log(2)/4 */
59
60 /*
61 * Now we compute exp(X) in, conceptually, three parts:
62 * - a pure power of two which we get from N>>2
63 * - exp(g) for g in [-log(2)/8,+log(2)/8], which we compute
64 * using a Remez-generated polynomial approximation
65 * - exp(k*log(2)/4) (aka 2^(k/4)) for k in [0..3], which we
66 * get from a lookup table in precision-and-a-half and
67 * multiply by g.
68 *
69 * We gain a bit of extra precision by the fact that actually
70 * our polynomial approximation gives us exp(g)-1, and we add
71 * the 1 back on by tweaking the prec-and-a-half multiplication
72 * step.
73 *
74 * Coefficients generated by the command
75
76 ./auxiliary/remez.jl --variable=g --suffix=f -- '-log(BigFloat(2))/8' '+log(BigFloat(2))/8' 3 0 '(expm1(x))/x'
77
78 */
79 Rg = g * (
80 9.999999412829185331953781321128516523408059996430919985217971370689774264850229e-01f+g*(4.999999608551332693833317084753864837160947932961832943901913087652889900683833e-01f+g*(1.667292360203016574303631953046104769969440903672618034272397630620346717392378e-01f+g*(4.168230895653321517750133783431970715648192153539929404872173693978116154823859e-02f)))
81 );
82
83 /*
84 * Do the table lookup and combine it with Rg, to get our final
85 * answer apart from the exponent.
86 */
87 {
88 static const float twotokover4top[4] = { 0x1p+0F, 0x1.306p+0F, 0x1.6ap+0F, 0x1.ae8p+0F };
89 static const float twotokover4bot[4] = { 0x0p+0F, 0x1.fc1464p-13F, 0x1.3cccfep-13F, 0x1.3f32b6p-13F };
90 static const float twotokover4all[4] = { 0x1p+0F, 0x1.306fep+0F, 0x1.6a09e6p+0F, 0x1.ae89fap+0F };
91 int index = (N & 3);
92 Rg = twotokover4top[index] + (twotokover4bot[index] + twotokover4all[index]*Rg);
93 N >>= 2;
94 }
95
96 /*
97 * Combine the output exponent and mantissa, and return.
98 */
99 if (__builtin_expect(edgecaseflag, 0)) {
100 Result = fhex(((N/2) << 23) + 0x3f800000);
101 Result *= Rg;
102 Result *= fhex(((N-N/2) << 23) + 0x3f800000);
103 /*
104 * Step not mentioned in C&W: set errno reliably.
105 */
106 if (fai(Result) == 0)
107 return MATHERR_EXPF_UFL(Result);
108 if (fai(Result) == 0x7f800000)
109 return MATHERR_EXPF_OFL(Result);
110 return FLOAT_CHECKDENORM(Result);
111 } else {
112 Result = fhex(N * 8388608.0f + (float)0x3f800000);
113 Result *= Rg;
114 }
115
116 return Result;
117 }
118