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