1*5a02ffc3SAndrew Turner /*
2*5a02ffc3SAndrew Turner * Single-precision vector erf(x) function.
3*5a02ffc3SAndrew Turner *
4*5a02ffc3SAndrew Turner * Copyright (c) 2023, Arm Limited.
5*5a02ffc3SAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*5a02ffc3SAndrew Turner */
7*5a02ffc3SAndrew Turner
8*5a02ffc3SAndrew Turner #include "v_math.h"
9*5a02ffc3SAndrew Turner #include "pl_sig.h"
10*5a02ffc3SAndrew Turner #include "pl_test.h"
11*5a02ffc3SAndrew Turner
12*5a02ffc3SAndrew Turner static const struct data
13*5a02ffc3SAndrew Turner {
14*5a02ffc3SAndrew Turner float32x4_t max, shift, third;
15*5a02ffc3SAndrew Turner #if WANT_SIMD_EXCEPT
16*5a02ffc3SAndrew Turner float32x4_t tiny_bound, scale_minus_one;
17*5a02ffc3SAndrew Turner #endif
18*5a02ffc3SAndrew Turner } data = {
19*5a02ffc3SAndrew Turner .max = V4 (3.9375), /* 4 - 8/128. */
20*5a02ffc3SAndrew Turner .shift = V4 (0x1p16f),
21*5a02ffc3SAndrew Turner .third = V4 (0x1.555556p-2f), /* 1/3. */
22*5a02ffc3SAndrew Turner #if WANT_SIMD_EXCEPT
23*5a02ffc3SAndrew Turner .tiny_bound = V4 (0x1p-62f),
24*5a02ffc3SAndrew Turner .scale_minus_one = V4 (0x1.06eba8p-3f), /* scale - 1.0. */
25*5a02ffc3SAndrew Turner #endif
26*5a02ffc3SAndrew Turner };
27*5a02ffc3SAndrew Turner
28*5a02ffc3SAndrew Turner #define AbsMask 0x7fffffff
29*5a02ffc3SAndrew Turner
30*5a02ffc3SAndrew Turner struct entry
31*5a02ffc3SAndrew Turner {
32*5a02ffc3SAndrew Turner float32x4_t erf;
33*5a02ffc3SAndrew Turner float32x4_t scale;
34*5a02ffc3SAndrew Turner };
35*5a02ffc3SAndrew Turner
36*5a02ffc3SAndrew Turner static inline struct entry
lookup(uint32x4_t i)37*5a02ffc3SAndrew Turner lookup (uint32x4_t i)
38*5a02ffc3SAndrew Turner {
39*5a02ffc3SAndrew Turner struct entry e;
40*5a02ffc3SAndrew Turner float64_t t0 = *((float64_t *) (__erff_data.tab + i[0]));
41*5a02ffc3SAndrew Turner float64_t t1 = *((float64_t *) (__erff_data.tab + i[1]));
42*5a02ffc3SAndrew Turner float64_t t2 = *((float64_t *) (__erff_data.tab + i[2]));
43*5a02ffc3SAndrew Turner float64_t t3 = *((float64_t *) (__erff_data.tab + i[3]));
44*5a02ffc3SAndrew Turner float32x4_t e1 = vreinterpretq_f32_f64 ((float64x2_t){ t0, t1 });
45*5a02ffc3SAndrew Turner float32x4_t e2 = vreinterpretq_f32_f64 ((float64x2_t){ t2, t3 });
46*5a02ffc3SAndrew Turner e.erf = vuzp1q_f32 (e1, e2);
47*5a02ffc3SAndrew Turner e.scale = vuzp2q_f32 (e1, e2);
48*5a02ffc3SAndrew Turner return e;
49*5a02ffc3SAndrew Turner }
50*5a02ffc3SAndrew Turner
51*5a02ffc3SAndrew Turner /* Single-precision implementation of vector erf(x).
52*5a02ffc3SAndrew Turner Approximation based on series expansion near x rounded to
53*5a02ffc3SAndrew Turner nearest multiple of 1/128.
54*5a02ffc3SAndrew Turner Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
55*5a02ffc3SAndrew Turner
56*5a02ffc3SAndrew Turner erf(x) ~ erf(r) + scale * d * [1 - r * d - 1/3 * d^2]
57*5a02ffc3SAndrew Turner
58*5a02ffc3SAndrew Turner Values of erf(r) and scale are read from lookup tables.
59*5a02ffc3SAndrew Turner For |x| > 3.9375, erf(|x|) rounds to 1.0f.
60*5a02ffc3SAndrew Turner
61*5a02ffc3SAndrew Turner Maximum error: 1.93 ULP
62*5a02ffc3SAndrew Turner _ZGVnN4v_erff(0x1.c373e6p-9) got 0x1.fd686cp-9
63*5a02ffc3SAndrew Turner want 0x1.fd6868p-9. */
V_NAME_F1(erf)64*5a02ffc3SAndrew Turner float32x4_t VPCS_ATTR V_NAME_F1 (erf) (float32x4_t x)
65*5a02ffc3SAndrew Turner {
66*5a02ffc3SAndrew Turner const struct data *dat = ptr_barrier (&data);
67*5a02ffc3SAndrew Turner
68*5a02ffc3SAndrew Turner #if WANT_SIMD_EXCEPT
69*5a02ffc3SAndrew Turner /* |x| < 2^-62. */
70*5a02ffc3SAndrew Turner uint32x4_t cmp = vcaltq_f32 (x, dat->tiny_bound);
71*5a02ffc3SAndrew Turner float32x4_t xm = x;
72*5a02ffc3SAndrew Turner /* If any lanes are special, mask them with 1 and retain a copy of x to allow
73*5a02ffc3SAndrew Turner special case handler to fix special lanes later. This is only necessary if
74*5a02ffc3SAndrew Turner fenv exceptions are to be triggered correctly. */
75*5a02ffc3SAndrew Turner if (unlikely (v_any_u32 (cmp)))
76*5a02ffc3SAndrew Turner x = vbslq_f32 (cmp, v_f32 (1), x);
77*5a02ffc3SAndrew Turner #endif
78*5a02ffc3SAndrew Turner
79*5a02ffc3SAndrew Turner float32x4_t a = vabsq_f32 (x);
80*5a02ffc3SAndrew Turner uint32x4_t a_gt_max = vcgtq_f32 (a, dat->max);
81*5a02ffc3SAndrew Turner
82*5a02ffc3SAndrew Turner /* Lookup erf(r) and scale(r) in tables, e.g. set erf(r) to 0 and scale to
83*5a02ffc3SAndrew Turner 2/sqrt(pi), when x reduced to r = 0. */
84*5a02ffc3SAndrew Turner float32x4_t shift = dat->shift;
85*5a02ffc3SAndrew Turner float32x4_t z = vaddq_f32 (a, shift);
86*5a02ffc3SAndrew Turner
87*5a02ffc3SAndrew Turner uint32x4_t i
88*5a02ffc3SAndrew Turner = vsubq_u32 (vreinterpretq_u32_f32 (z), vreinterpretq_u32_f32 (shift));
89*5a02ffc3SAndrew Turner i = vminq_u32 (i, v_u32 (512));
90*5a02ffc3SAndrew Turner struct entry e = lookup (i);
91*5a02ffc3SAndrew Turner
92*5a02ffc3SAndrew Turner float32x4_t r = vsubq_f32 (z, shift);
93*5a02ffc3SAndrew Turner
94*5a02ffc3SAndrew Turner /* erf(x) ~ erf(r) + scale * d * (1 - r * d - 1/3 * d^2). */
95*5a02ffc3SAndrew Turner float32x4_t d = vsubq_f32 (a, r);
96*5a02ffc3SAndrew Turner float32x4_t d2 = vmulq_f32 (d, d);
97*5a02ffc3SAndrew Turner float32x4_t y = vfmaq_f32 (r, dat->third, d);
98*5a02ffc3SAndrew Turner y = vfmaq_f32 (e.erf, e.scale, vfmsq_f32 (d, d2, y));
99*5a02ffc3SAndrew Turner
100*5a02ffc3SAndrew Turner /* Solves the |x| = inf case. */
101*5a02ffc3SAndrew Turner y = vbslq_f32 (a_gt_max, v_f32 (1.0f), y);
102*5a02ffc3SAndrew Turner
103*5a02ffc3SAndrew Turner /* Copy sign. */
104*5a02ffc3SAndrew Turner y = vbslq_f32 (v_u32 (AbsMask), y, x);
105*5a02ffc3SAndrew Turner
106*5a02ffc3SAndrew Turner #if WANT_SIMD_EXCEPT
107*5a02ffc3SAndrew Turner if (unlikely (v_any_u32 (cmp)))
108*5a02ffc3SAndrew Turner return vbslq_f32 (cmp, vfmaq_f32 (xm, dat->scale_minus_one, xm), y);
109*5a02ffc3SAndrew Turner #endif
110*5a02ffc3SAndrew Turner return y;
111*5a02ffc3SAndrew Turner }
112*5a02ffc3SAndrew Turner
113*5a02ffc3SAndrew Turner PL_SIG (V, F, 1, erf, -4.0, 4.0)
114*5a02ffc3SAndrew Turner PL_TEST_ULP (V_NAME_F1 (erf), 1.43)
115*5a02ffc3SAndrew Turner PL_TEST_EXPECT_FENV (V_NAME_F1 (erf), WANT_SIMD_EXCEPT)
116*5a02ffc3SAndrew Turner PL_TEST_SYM_INTERVAL (V_NAME_F1 (erf), 0, 3.9375, 40000)
117*5a02ffc3SAndrew Turner PL_TEST_SYM_INTERVAL (V_NAME_F1 (erf), 3.9375, inf, 40000)
118*5a02ffc3SAndrew Turner PL_TEST_SYM_INTERVAL (V_NAME_F1 (erf), 0, inf, 40000)
119