diff options
Diffstat (limited to 'pl/math/log1pf_2u1.c')
| -rw-r--r-- | pl/math/log1pf_2u1.c | 165 |
1 files changed, 165 insertions, 0 deletions
diff --git a/pl/math/log1pf_2u1.c b/pl/math/log1pf_2u1.c new file mode 100644 index 0000000..fcfd05a --- /dev/null +++ b/pl/math/log1pf_2u1.c @@ -0,0 +1,165 @@ +/* + * Single-precision log(1+x) function. + * + * Copyright (c) 2022-2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "hornerf.h" +#include "math_config.h" +#include "pl_sig.h" +#include "pl_test.h" + +#define Ln2 (0x1.62e43p-1f) +#define SignMask (0x80000000) + +/* Biased exponent of the largest float m for which m^8 underflows. */ +#define M8UFLOW_BOUND_BEXP 112 +/* Biased exponent of the largest float for which we just return x. */ +#define TINY_BOUND_BEXP 103 + +#define C(i) __log1pf_data.coeffs[i] + +static inline float +eval_poly (float m, uint32_t e) +{ +#ifdef LOG1PF_2U5 + + /* 2.5 ulp variant. Approximate log(1+m) on [-0.25, 0.5] using + slightly modified Estrin scheme (no x^0 term, and x term is just x). */ + float p_12 = fmaf (m, C (1), C (0)); + float p_34 = fmaf (m, C (3), C (2)); + float p_56 = fmaf (m, C (5), C (4)); + float p_78 = fmaf (m, C (7), C (6)); + + float m2 = m * m; + float p_02 = fmaf (m2, p_12, m); + float p_36 = fmaf (m2, p_56, p_34); + float p_79 = fmaf (m2, C (8), p_78); + + float m4 = m2 * m2; + float p_06 = fmaf (m4, p_36, p_02); + + if (unlikely (e < M8UFLOW_BOUND_BEXP)) + return p_06; + + float m8 = m4 * m4; + return fmaf (m8, p_79, p_06); + +#elif defined(LOG1PF_1U3) + + /* 1.3 ulp variant. Approximate log(1+m) on [-0.25, 0.5] using Horner + scheme. Our polynomial approximation for log1p has the form + x + C1 * x^2 + C2 * x^3 + C3 * x^4 + ... + Hence approximation has the form m + m^2 * P(m) + where P(x) = C1 + C2 * x + C3 * x^2 + ... . */ + return fmaf (m, m * HORNER_8 (m, C), m); + +#else +#error No log1pf approximation exists with the requested precision. Options are 13 or 25. +#endif +} + +static inline uint32_t +biased_exponent (uint32_t ix) +{ + return (ix & 0x7f800000) >> 23; +} + +/* log1pf approximation using polynomial on reduced interval. Worst-case error + when using Estrin is roughly 2.02 ULP: + log1pf(0x1.21e13ap-2) got 0x1.fe8028p-3 want 0x1.fe802cp-3. */ +float +log1pf (float x) +{ + uint32_t ix = asuint (x); + uint32_t ia = ix & ~SignMask; + uint32_t ia12 = ia >> 20; + uint32_t e = biased_exponent (ix); + + /* Handle special cases first. */ + if (unlikely (ia12 >= 0x7f8 || ix >= 0xbf800000 || ix == 0x80000000 + || e <= TINY_BOUND_BEXP)) + { + if (ix == 0xff800000) + { + /* x == -Inf => log1pf(x) = NaN. */ + return NAN; + } + if ((ix == 0x7f800000 || e <= TINY_BOUND_BEXP) && ia12 <= 0x7f8) + { + /* |x| < TinyBound => log1p(x) = x. + x == Inf => log1pf(x) = Inf. */ + return x; + } + if (ix == 0xbf800000) + { + /* x == -1.0 => log1pf(x) = -Inf. */ + return __math_divzerof (-1); + } + if (ia12 >= 0x7f8) + { + /* x == +/-NaN => log1pf(x) = NaN. */ + return __math_invalidf (asfloat (ia)); + } + /* x < -1.0 => log1pf(x) = NaN. */ + return __math_invalidf (x); + } + + /* With x + 1 = t * 2^k (where t = m + 1 and k is chosen such that m + is in [-0.25, 0.5]): + log1p(x) = log(t) + log(2^k) = log1p(m) + k*log(2). + + We approximate log1p(m) with a polynomial, then scale by + k*log(2). Instead of doing this directly, we use an intermediate + scale factor s = 4*k*log(2) to ensure the scale is representable + as a normalised fp32 number. */ + + if (ix <= 0x3f000000 || ia <= 0x3e800000) + { + /* If x is in [-0.25, 0.5] then we can shortcut all the logic + below, as k = 0 and m = x. All we need is to return the + polynomial. */ + return eval_poly (x, e); + } + + float m = x + 1.0f; + + /* k is used scale the input. 0x3f400000 is chosen as we are trying to + reduce x to the range [-0.25, 0.5]. Inside this range, k is 0. + Outside this range, if k is reinterpreted as (NOT CONVERTED TO) float: + let k = sign * 2^p where sign = -1 if x < 0 + 1 otherwise + and p is a negative integer whose magnitude increases with the + magnitude of x. */ + int k = (asuint (m) - 0x3f400000) & 0xff800000; + + /* By using integer arithmetic, we obtain the necessary scaling by + subtracting the unbiased exponent of k from the exponent of x. */ + float m_scale = asfloat (asuint (x) - k); + + /* Scale up to ensure that the scale factor is representable as normalised + fp32 number (s in [2**-126,2**26]), and scale m down accordingly. */ + float s = asfloat (asuint (4.0f) - k); + m_scale = m_scale + fmaf (0.25f, s, -1.0f); + + float p = eval_poly (m_scale, biased_exponent (asuint (m_scale))); + + /* The scale factor to be applied back at the end - by multiplying float(k) + by 2^-23 we get the unbiased exponent of k. */ + float scale_back = (float) k * 0x1.0p-23f; + + /* Apply the scaling back. */ + return fmaf (scale_back, Ln2, p); +} + +PL_SIG (S, F, 1, log1p, -0.9, 10.0) +PL_TEST_ULP (log1pf, 1.52) +PL_TEST_INTERVAL (log1pf, -10.0, 10.0, 10000) +PL_TEST_INTERVAL (log1pf, 0.0, 0x1p-23, 50000) +PL_TEST_INTERVAL (log1pf, 0x1p-23, 0.001, 50000) +PL_TEST_INTERVAL (log1pf, 0.001, 1.0, 50000) +PL_TEST_INTERVAL (log1pf, 0.0, -0x1p-23, 50000) +PL_TEST_INTERVAL (log1pf, -0x1p-23, -0.001, 50000) +PL_TEST_INTERVAL (log1pf, -0.001, -1.0, 50000) +PL_TEST_INTERVAL (log1pf, -1.0, inf, 5000) |