diff options
Diffstat (limited to 'pl/math/expf.c')
-rw-r--r-- | pl/math/expf.c | 76 |
1 files changed, 76 insertions, 0 deletions
diff --git a/pl/math/expf.c b/pl/math/expf.c new file mode 100644 index 0000000..c325e45 --- /dev/null +++ b/pl/math/expf.c @@ -0,0 +1,76 @@ +/* + * Single-precision e^x function. + * + * Copyright (c) 2017-2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include <math.h> +#include <stdint.h> +#include "math_config.h" + +/* +EXPF_TABLE_BITS = 5 +EXPF_POLY_ORDER = 3 + +ULP error: 0.502 (nearest rounding.) +Relative error: 1.69 * 2^-34 in [-ln2/64, ln2/64] (before rounding.) +Wrong count: 170635 (all nearest rounding wrong results with fma.) +Non-nearest ULP error: 1 (rounded ULP error) +*/ + +#define N (1 << EXPF_TABLE_BITS) +#define InvLn2N __expf_data.invln2_scaled +#define T __expf_data.tab +#define C __expf_data.poly_scaled + +static inline uint32_t +top12 (float x) +{ + return asuint (x) >> 20; +} + +float +optr_aor_exp_f32 (float x) +{ + uint32_t abstop; + uint64_t ki, t; + /* double_t for better performance on targets with FLT_EVAL_METHOD==2. */ + double_t kd, xd, z, r, r2, y, s; + + xd = (double_t) x; + abstop = top12 (x) & 0x7ff; + if (unlikely (abstop >= top12 (88.0f))) + { + /* |x| >= 88 or x is nan. */ + if (asuint (x) == asuint (-INFINITY)) + return 0.0f; + if (abstop >= top12 (INFINITY)) + return x + x; + if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */ + return __math_oflowf (0); + if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */ + return __math_uflowf (0); + } + + /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ + z = InvLn2N * xd; + + /* Round and convert z to int, the result is in [-150*N, 128*N] and + ideally nearest int is used, otherwise the magnitude of r can be + bigger which gives larger approximation error. */ + kd = roundtoint (z); + ki = converttoint (z); + r = z - kd; + + /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = T[ki % N]; + t += ki << (52 - EXPF_TABLE_BITS); + s = asdouble (t); + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return eval_as_float (y); +} |