1 /*
2  * Single-precision 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 "math_config.h"
9 #include "pl_sig.h"
10 #include "pl_test.h"
11 
12 #define TwoOverSqrtPiMinusOne 0x1.06eba8p-3f
13 #define Shift 0x1p16f
14 #define OneThird 0x1.555556p-2f
15 
16 /* Fast erff approximation based on series expansion near x rounded to
17    nearest multiple of 1/128.
18    Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
19 
20    erf(x) ~ erf(r)
21      + scale * d * [
22        + 1
23        - r d
24        + 1/3 (2 r^2 - 1) d^2
25        - 1/6 (r (2 r^2 - 3) ) d^3
26        + 1/30 (4 r^4 - 12 r^2 + 3) d^4
27      ]
28 
29    This single precision implementation uses only the following terms:
30 
31    erf(x) ~ erf(r) + scale * d * [1 - r * d - 1/3 * d^2]
32 
33    Values of erf(r) and scale are read from lookup tables.
34    For |x| > 3.9375, erf(|x|) rounds to 1.0f.
35 
36    Maximum error: 1.93 ULP
37    erff(0x1.c373e6p-9) got 0x1.fd686cp-9
38 		      want 0x1.fd6868p-9.  */
39 float
40 erff (float x)
41 {
42   /* Get absolute value and sign.  */
43   uint32_t ix = asuint (x);
44   uint32_t ia = ix & 0x7fffffff;
45   uint32_t sign = ix & ~0x7fffffff;
46 
47   /* |x| < 0x1p-62. Triggers exceptions.  */
48   if (unlikely (ia < 0x20800000))
49     return fmaf (TwoOverSqrtPiMinusOne, x, x);
50 
51   if (ia < 0x407b8000) /* |x| <  4 - 8 / 128 = 3.9375.  */
52     {
53       /* Lookup erf(r) and scale(r) in tables, e.g. set erf(r) to 0 and scale
54 	 to 2/sqrt(pi), when x reduced to r = 0.  */
55       float a = asfloat (ia);
56       float z = a + Shift;
57       uint32_t i = asuint (z) - asuint (Shift);
58       float r = z - Shift;
59       float erfr = __erff_data.tab[i].erf;
60       float scale = __erff_data.tab[i].scale;
61 
62       /* erf(x) ~ erf(r) + scale * d * (1 - r * d - 1/3 * d^2).  */
63       float d = a - r;
64       float d2 = d * d;
65       float y = -fmaf (OneThird, d, r);
66       y = fmaf (fmaf (y, d2, d), scale, erfr);
67       return asfloat (asuint (y) | sign);
68     }
69 
70   /* Special cases : erff(nan)=nan, erff(+inf)=+1 and erff(-inf)=-1.  */
71   if (unlikely (ia >= 0x7f800000))
72     return (1.0f - (float) (sign >> 30)) + 1.0f / x;
73 
74   /* Boring domain (|x| >= 4.0).  */
75   return asfloat (sign | asuint (1.0f));
76 }
77 
78 PL_SIG (S, F, 1, erf, -4.0, 4.0)
79 PL_TEST_ULP (erff, 1.43)
80 PL_TEST_SYM_INTERVAL (erff, 0, 3.9375, 40000)
81 PL_TEST_SYM_INTERVAL (erff, 3.9375, inf, 40000)
82 PL_TEST_SYM_INTERVAL (erff, 0, inf, 40000)
83