1 /*
2 * Double-precision vector erf(x) function.
3 *
4 * Copyright (c) 2019-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "v_math.h"
9 #include "include/mathlib.h"
10 #include "math_config.h"
11 #include "pl_sig.h"
12 #include "pl_test.h"
13
14 #if V_SUPPORTED
15
16 #define AbsMask v_u64 (0x7fffffffffffffff)
17 #define AbsXMax v_f64 (0x1.8p+2)
18 #define Scale v_f64 (0x1p+3)
19
20 /* Special cases (fall back to scalar calls). */
21 VPCS_ATTR
22 NOINLINE static v_f64_t
specialcase(v_f64_t x,v_f64_t y,v_u64_t cmp)23 specialcase (v_f64_t x, v_f64_t y, v_u64_t cmp)
24 {
25 return v_call_f64 (erf, x, y, cmp);
26 }
27
28 /* A structure to perform look-up in coeffs and other parameter tables. */
29 struct entry
30 {
31 v_f64_t P[V_ERF_NCOEFFS];
32 v_f64_t shift;
33 };
34
35 static inline struct entry
lookup(v_u64_t i)36 lookup (v_u64_t i)
37 {
38 struct entry e;
39 #ifdef SCALAR
40 for (int j = 0; j < V_ERF_NCOEFFS; ++j)
41 e.P[j] = __v_erf_data.coeffs[j][i];
42 e.shift = __v_erf_data.shifts[i];
43 #else
44 for (int j = 0; j < V_ERF_NCOEFFS; ++j)
45 {
46 e.P[j][0] = __v_erf_data.coeffs[j][i[0]];
47 e.P[j][1] = __v_erf_data.coeffs[j][i[1]];
48 }
49 e.shift[0] = __v_erf_data.shifts[i[0]];
50 e.shift[1] = __v_erf_data.shifts[i[1]];
51 #endif
52 return e;
53 }
54
55 /* Optimized double precision vector error function erf. Maximum
56 observed error is 1.75 ULP, in [0.110, 0.111]:
57 verf(0x1.c5e0c2d5d0543p-4) got 0x1.fe0ed62a54987p-4
58 want 0x1.fe0ed62a54985p-4. */
59 VPCS_ATTR
V_NAME(erf)60 v_f64_t V_NAME (erf) (v_f64_t x)
61 {
62 /* Handle both inf/nan as well as small values (|x|<2^-28)
63 If any condition in the lane is true then a loop over
64 scalar calls will be performed. */
65 v_u64_t ix = v_as_u64_f64 (x);
66 v_u64_t atop = (ix >> 48) & v_u64 (0x7fff);
67 v_u64_t special_case
68 = v_cond_u64 (atop - v_u64 (0x3e30) >= v_u64 (0x7ff0 - 0x3e30));
69
70 /* Get sign and absolute value. */
71 v_u64_t sign = v_as_u64_f64 (x) & ~AbsMask;
72 v_f64_t a = v_min_f64 (v_abs_f64 (x), AbsXMax);
73
74 /* Compute index by truncating 8 * a with a=|x| saturated to 6.0. */
75
76 #ifdef SCALAR
77 v_u64_t i = v_trunc_u64 (a * Scale);
78 #else
79 v_u64_t i = vcvtq_n_u64_f64 (a, 3);
80 #endif
81 /* Get polynomial coefficients and shift parameter using lookup. */
82 struct entry dat = lookup (i);
83
84 /* Evaluate polynomial on transformed argument. */
85 v_f64_t z = v_fma_f64 (a, Scale, dat.shift);
86
87 v_f64_t r1 = v_fma_f64 (z, dat.P[1], dat.P[0]);
88 v_f64_t r2 = v_fma_f64 (z, dat.P[3], dat.P[2]);
89 v_f64_t r3 = v_fma_f64 (z, dat.P[5], dat.P[4]);
90 v_f64_t r4 = v_fma_f64 (z, dat.P[7], dat.P[6]);
91 v_f64_t r5 = v_fma_f64 (z, dat.P[9], dat.P[8]);
92
93 v_f64_t z2 = z * z;
94 v_f64_t y = v_fma_f64 (z2, r5, r4);
95 y = v_fma_f64 (z2, y, r3);
96 y = v_fma_f64 (z2, y, r2);
97 y = v_fma_f64 (z2, y, r1);
98
99 /* y=erf(x) if x>0, -erf(-x) otherwise. */
100 y = v_as_f64_u64 (v_as_u64_f64 (y) ^ sign);
101
102 if (unlikely (v_any_u64 (special_case)))
103 return specialcase (x, y, special_case);
104 return y;
105 }
106 VPCS_ALIAS
107
108 PL_SIG (V, D, 1, erf, -6.0, 6.0)
109 PL_TEST_ULP (V_NAME (erf), 1.26)
110 PL_TEST_INTERVAL (V_NAME (erf), 0, 0xffff0000, 10000)
111 PL_TEST_INTERVAL (V_NAME (erf), 0x1p-127, 0x1p-26, 40000)
112 PL_TEST_INTERVAL (V_NAME (erf), -0x1p-127, -0x1p-26, 40000)
113 PL_TEST_INTERVAL (V_NAME (erf), 0x1p-26, 0x1p3, 40000)
114 PL_TEST_INTERVAL (V_NAME (erf), -0x1p-26, -0x1p3, 40000)
115 PL_TEST_INTERVAL (V_NAME (erf), 0, inf, 40000)
116 #endif
117