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