2  * Double-precision e^x function.
4  * Copyright (c) 2018-2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
21 #define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
22 #define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
23 #define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
24 #define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
25 #define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
31    a double.  (int32_t)KI is the k used in the argument reduction and exponent
34 static inline double
42       sbits -= 1009ull << 52;  in specialcase()
47   /* k < 0, need special care in the subnormal range.  */  in specialcase()
54 	 range to avoid double rounding that can cause 0.5+E/2 ulp error where  in specialcase()
55 	 E is the worst-case ulp error outside the subnormal range.  So this  in specialcase()
56 	 is only useful if the goal is better than 1 ulp worst-case error.  */  in specialcase()
58       lo = scale - y + scale * tmp;  in specialcase()
60       lo = 1.0 - hi + y + lo;  in specialcase()
61       y = eval_as_double (hi + lo) - 1.0;  in specialcase()
62       /* Avoid -0.0 with downward rounding.  */  in specialcase()
66       force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);  in specialcase()
68   y = 0x1p-1022 * y;  in specialcase()
72 /* Top 12 bits of a double (sign and exponent bits).  */
74 top12 (double x)  in top12()
79 /* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
81 static inline double
82 exp_inline (double x, double xtail)  in exp_inline()
90   if (unlikely (abstop - top12 (0x1p-54) >= top12 (512.0) - top12 (0x1p-54)))  in exp_inline()
92       if (abstop - top12 (0x1p-54) >= 0x80000000)  in exp_inline()
98 	  if (asuint64 (x) == asuint64 (-INFINITY))  in exp_inline()
111   /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */  in exp_inline()
112   /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */  in exp_inline()
118   /* z - kd is in [-0.5-2^-16, 0.5] in all rounding modes.  */  in exp_inline()
123   /* z - kd is in [-1, 1] in non-nearest rounding modes.  */  in exp_inline()
126   kd -= Shift;  in exp_inline()
129   /* The code assumes 2^-200 < |xtail| < 2^-8/N.  */  in exp_inline()
134   top = ki << (52 - EXP_TABLE_BITS);  in exp_inline()
136   /* This is only a valid scale when -1023*N < k < 1024*N.  */  in exp_inline()
138   /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */  in exp_inline()
153   /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there  in exp_inline()
158 double
159 exp (double x)  in exp()
168 long double expl (long double x) { return exp (x); }  in strong_alias()
172 TEST_SIG (S, D, 1, exp, -9.9, 9.9)
176 TEST_SYM_INTERVAL (exp, 0x1p-6, 0x1p6, 400000)