diff options
Diffstat (limited to 'pl/math/v_log1p_2u5.c')
-rw-r--r-- | pl/math/v_log1p_2u5.c | 120 |
1 files changed, 120 insertions, 0 deletions
diff --git a/pl/math/v_log1p_2u5.c b/pl/math/v_log1p_2u5.c new file mode 100644 index 0000000..e482910 --- /dev/null +++ b/pl/math/v_log1p_2u5.c @@ -0,0 +1,120 @@ +/* + * Double-precision vector log(1+x) function. + * + * Copyright (c) 2022-2023, Arm Limited. + * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception + */ + +#include "v_math.h" +#include "estrin.h" +#include "pl_sig.h" +#include "pl_test.h" + +#if V_SUPPORTED + +#define Ln2Hi v_f64 (0x1.62e42fefa3800p-1) +#define Ln2Lo v_f64 (0x1.ef35793c76730p-45) +#define HfRt2Top 0x3fe6a09e00000000 /* top32(asuint64(sqrt(2)/2)) << 32. */ +#define OneMHfRt2Top \ + 0x00095f6200000000 /* (top32(asuint64(1)) - top32(asuint64(sqrt(2)/2))) \ + << 32. */ +#define OneTop12 0x3ff +#define BottomMask 0xffffffff +#define AbsMask 0x7fffffffffffffff +#define C(i) v_f64 (__log1p_data.coeffs[i]) + +static inline v_f64_t +eval_poly (v_f64_t f) +{ + v_f64_t f2 = f * f; + v_f64_t f4 = f2 * f2; + v_f64_t f8 = f4 * f4; + return ESTRIN_18 (f, f2, f4, f8, f8 * f8, C); +} + +VPCS_ATTR +NOINLINE static v_f64_t +specialcase (v_f64_t x, v_f64_t y, v_u64_t special) +{ + return v_call_f64 (log1p, x, y, special); +} + +/* Vector log1p approximation using polynomial on reduced interval. Routine is a + modification of the algorithm used in scalar log1p, with no shortcut for k=0 + and no narrowing for f and k. Maximum observed error is 2.46 ULP: + __v_log1p(0x1.654a1307242a4p+11) got 0x1.fd5565fb590f4p+2 + want 0x1.fd5565fb590f6p+2 . */ +VPCS_ATTR v_f64_t V_NAME (log1p) (v_f64_t x) +{ + v_u64_t ix = v_as_u64_f64 (x); + v_u64_t ia = ix & AbsMask; + v_u64_t special + = v_cond_u64 ((ia >= v_u64 (0x7ff0000000000000)) + | (ix >= 0xbff0000000000000) | (ix == 0x8000000000000000)); + +#if WANT_SIMD_EXCEPT + if (unlikely (v_any_u64 (special))) + x = v_sel_f64 (special, v_f64 (0), x); +#endif + + /* With x + 1 = t * 2^k (where t = f + 1 and k is chosen such that f + is in [sqrt(2)/2, sqrt(2)]): + log1p(x) = k*log(2) + log1p(f). + + f may not be representable exactly, so we need a correction term: + let m = round(1 + x), c = (1 + x) - m. + c << m: at very small x, log1p(x) ~ x, hence: + log(1+x) - log(m) ~ c/m. + + We therefore calculate log1p(x) by k*log2 + log1p(f) + c/m. */ + + /* Obtain correctly scaled k by manipulation in the exponent. + The scalar algorithm casts down to 32-bit at this point to calculate k and + u_red. We stay in double-width to obtain f and k, using the same constants + as the scalar algorithm but shifted left by 32. */ + v_f64_t m = x + 1; + v_u64_t mi = v_as_u64_f64 (m); + v_u64_t u = mi + OneMHfRt2Top; + + v_s64_t ki = v_as_s64_u64 (u >> 52) - OneTop12; + v_f64_t k = v_to_f64_s64 (ki); + + /* Reduce x to f in [sqrt(2)/2, sqrt(2)]. */ + v_u64_t utop = (u & 0x000fffff00000000) + HfRt2Top; + v_u64_t u_red = utop | (mi & BottomMask); + v_f64_t f = v_as_f64_u64 (u_red) - 1; + + /* Correction term c/m. */ + v_f64_t cm = (x - (m - 1)) / m; + + /* Approximate log1p(x) on the reduced input using a polynomial. Because + log1p(0)=0 we choose an approximation of the form: + x + C0*x^2 + C1*x^3 + C2x^4 + ... + Hence approximation has the form f + f^2 * P(f) + where P(x) = C0 + C1*x + C2x^2 + ... + Assembling this all correctly is dealt with at the final step. */ + v_f64_t p = eval_poly (f); + + v_f64_t ylo = v_fma_f64 (k, Ln2Lo, cm); + v_f64_t yhi = v_fma_f64 (k, Ln2Hi, f); + v_f64_t y = v_fma_f64 (f * f, p, ylo + yhi); + + if (unlikely (v_any_u64 (special))) + return specialcase (v_as_f64_u64 (ix), y, special); + + return y; +} +VPCS_ALIAS + +PL_SIG (V, D, 1, log1p, -0.9, 10.0) +PL_TEST_ULP (V_NAME (log1p), 1.97) +PL_TEST_EXPECT_FENV (V_NAME (log1p), WANT_SIMD_EXCEPT) +PL_TEST_INTERVAL (V_NAME (log1p), -10.0, 10.0, 10000) +PL_TEST_INTERVAL (V_NAME (log1p), 0.0, 0x1p-23, 50000) +PL_TEST_INTERVAL (V_NAME (log1p), 0x1p-23, 0.001, 50000) +PL_TEST_INTERVAL (V_NAME (log1p), 0.001, 1.0, 50000) +PL_TEST_INTERVAL (V_NAME (log1p), 0.0, -0x1p-23, 50000) +PL_TEST_INTERVAL (V_NAME (log1p), -0x1p-23, -0.001, 50000) +PL_TEST_INTERVAL (V_NAME (log1p), -0.001, -1.0, 50000) +PL_TEST_INTERVAL (V_NAME (log1p), -1.0, inf, 5000) +#endif |