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