1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision vector log(1+x) function.
3*f3087befSAndrew Turner *
4*f3087befSAndrew Turner * Copyright (c) 2022-2024, Arm Limited.
5*f3087befSAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*f3087befSAndrew Turner */
7*f3087befSAndrew Turner
8*f3087befSAndrew Turner #include "v_math.h"
9*f3087befSAndrew Turner #include "test_sig.h"
10*f3087befSAndrew Turner #include "test_defs.h"
11*f3087befSAndrew Turner
12*f3087befSAndrew Turner #define WANT_V_LOG1P_K0_SHORTCUT 0
13*f3087befSAndrew Turner #include "v_log1p_inline.h"
14*f3087befSAndrew Turner
15*f3087befSAndrew Turner const static struct data
16*f3087befSAndrew Turner {
17*f3087befSAndrew Turner struct v_log1p_data d;
18*f3087befSAndrew Turner uint64x2_t inf, minus_one;
19*f3087befSAndrew Turner } data = { .d = V_LOG1P_CONSTANTS_TABLE,
20*f3087befSAndrew Turner .inf = V2 (0x7ff0000000000000),
21*f3087befSAndrew Turner .minus_one = V2 (0xbff0000000000000) };
22*f3087befSAndrew Turner
23*f3087befSAndrew Turner #define BottomMask v_u64 (0xffffffff)
24*f3087befSAndrew Turner
25*f3087befSAndrew Turner static float64x2_t NOINLINE VPCS_ATTR
special_case(float64x2_t x,uint64x2_t cmp,const struct data * d)26*f3087befSAndrew Turner special_case (float64x2_t x, uint64x2_t cmp, const struct data *d)
27*f3087befSAndrew Turner {
28*f3087befSAndrew Turner /* Side-step special lanes so fenv exceptions are not triggered
29*f3087befSAndrew Turner inadvertently. */
30*f3087befSAndrew Turner float64x2_t x_nospecial = v_zerofy_f64 (x, cmp);
31*f3087befSAndrew Turner return v_call_f64 (log1p, x, log1p_inline (x_nospecial, &d->d), cmp);
32*f3087befSAndrew Turner }
33*f3087befSAndrew Turner
34*f3087befSAndrew Turner /* Vector log1p approximation using polynomial on reduced interval. Routine is
35*f3087befSAndrew Turner a modification of the algorithm used in scalar log1p, with no shortcut for
36*f3087befSAndrew Turner k=0 and no narrowing for f and k. Maximum observed error is 2.45 ULP:
37*f3087befSAndrew Turner _ZGVnN2v_log1p(0x1.658f7035c4014p+11) got 0x1.fd61d0727429dp+2
38*f3087befSAndrew Turner want 0x1.fd61d0727429fp+2 . */
V_NAME_D1(log1p)39*f3087befSAndrew Turner VPCS_ATTR float64x2_t V_NAME_D1 (log1p) (float64x2_t x)
40*f3087befSAndrew Turner {
41*f3087befSAndrew Turner const struct data *d = ptr_barrier (&data);
42*f3087befSAndrew Turner uint64x2_t ix = vreinterpretq_u64_f64 (x);
43*f3087befSAndrew Turner uint64x2_t ia = vreinterpretq_u64_f64 (vabsq_f64 (x));
44*f3087befSAndrew Turner
45*f3087befSAndrew Turner uint64x2_t special_cases
46*f3087befSAndrew Turner = vorrq_u64 (vcgeq_u64 (ia, d->inf), vcgeq_u64 (ix, d->minus_one));
47*f3087befSAndrew Turner
48*f3087befSAndrew Turner if (unlikely (v_any_u64 (special_cases)))
49*f3087befSAndrew Turner return special_case (x, special_cases, d);
50*f3087befSAndrew Turner
51*f3087befSAndrew Turner return log1p_inline (x, &d->d);
52*f3087befSAndrew Turner }
53*f3087befSAndrew Turner
54*f3087befSAndrew Turner TEST_SIG (V, D, 1, log1p, -0.9, 10.0)
55*f3087befSAndrew Turner TEST_ULP (V_NAME_D1 (log1p), 1.95)
56*f3087befSAndrew Turner TEST_DISABLE_FENV_IF_NOT (V_NAME_D1 (log1p), WANT_SIMD_EXCEPT)
57*f3087befSAndrew Turner TEST_SYM_INTERVAL (V_NAME_D1 (log1p), 0.0, 0x1p-23, 50000)
58*f3087befSAndrew Turner TEST_SYM_INTERVAL (V_NAME_D1 (log1p), 0x1p-23, 0.001, 50000)
59*f3087befSAndrew Turner TEST_SYM_INTERVAL (V_NAME_D1 (log1p), 0.001, 1.0, 50000)
60*f3087befSAndrew Turner TEST_INTERVAL (V_NAME_D1 (log1p), 1, inf, 40000)
61*f3087befSAndrew Turner TEST_INTERVAL (V_NAME_D1 (log1p), -1.0, -inf, 500)
62