1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision 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 "math_config.h"
9*f3087befSAndrew Turner #include "test_sig.h"
10*f3087befSAndrew Turner #include "test_defs.h"
11*f3087befSAndrew Turner
12*f3087befSAndrew Turner #define TwoOverSqrtPiMinusOne 0x1.06eba8214db69p-3
13*f3087befSAndrew Turner #define Shift 0x1p45
14*f3087befSAndrew Turner
15*f3087befSAndrew Turner /* Polynomial coefficients. */
16*f3087befSAndrew Turner #define OneThird 0x1.5555555555555p-2
17*f3087befSAndrew Turner #define TwoThird 0x1.5555555555555p-1
18*f3087befSAndrew Turner
19*f3087befSAndrew Turner #define TwoOverFifteen 0x1.1111111111111p-3
20*f3087befSAndrew Turner #define TwoOverFive 0x1.999999999999ap-2
21*f3087befSAndrew Turner #define Tenth 0x1.999999999999ap-4
22*f3087befSAndrew Turner
23*f3087befSAndrew Turner #define TwoOverNine 0x1.c71c71c71c71cp-3
24*f3087befSAndrew Turner #define TwoOverFortyFive 0x1.6c16c16c16c17p-5
25*f3087befSAndrew Turner #define Sixth 0x1.555555555555p-3
26*f3087befSAndrew Turner
27*f3087befSAndrew Turner /* Fast erf approximation based on series expansion near x rounded to
28*f3087befSAndrew Turner nearest multiple of 1/128.
29*f3087befSAndrew Turner Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
30*f3087befSAndrew Turner
31*f3087befSAndrew Turner erf(x) ~ erf(r)
32*f3087befSAndrew Turner + scale * d * [
33*f3087befSAndrew Turner + 1
34*f3087befSAndrew Turner - r d
35*f3087befSAndrew Turner + 1/3 (2 r^2 - 1) d^2
36*f3087befSAndrew Turner - 1/6 (r (2 r^2 - 3)) d^3
37*f3087befSAndrew Turner + 1/30 (4 r^4 - 12 r^2 + 3) d^4
38*f3087befSAndrew Turner - 1/90 (4 r^4 - 20 r^2 + 15) d^5
39*f3087befSAndrew Turner ]
40*f3087befSAndrew Turner
41*f3087befSAndrew Turner Maximum measure error: 2.29 ULP
42*f3087befSAndrew Turner erf(-0x1.00003c924e5d1p-8) got -0x1.20dd59132ebadp-8
43*f3087befSAndrew Turner want -0x1.20dd59132ebafp-8. */
44*f3087befSAndrew Turner double
arm_math_erf(double x)45*f3087befSAndrew Turner arm_math_erf (double x)
46*f3087befSAndrew Turner {
47*f3087befSAndrew Turner /* Get absolute value and sign. */
48*f3087befSAndrew Turner uint64_t ix = asuint64 (x);
49*f3087befSAndrew Turner uint64_t ia = ix & 0x7fffffffffffffff;
50*f3087befSAndrew Turner uint64_t sign = ix & ~0x7fffffffffffffff;
51*f3087befSAndrew Turner
52*f3087befSAndrew Turner /* |x| < 0x1p-508. Triggers exceptions. */
53*f3087befSAndrew Turner if (unlikely (ia < 0x2030000000000000))
54*f3087befSAndrew Turner return fma (TwoOverSqrtPiMinusOne, x, x);
55*f3087befSAndrew Turner
56*f3087befSAndrew Turner if (ia < 0x4017f80000000000) /* |x| < 6 - 1 / 128 = 5.9921875. */
57*f3087befSAndrew Turner {
58*f3087befSAndrew Turner /* Set r to multiple of 1/128 nearest to |x|. */
59*f3087befSAndrew Turner double a = asdouble (ia);
60*f3087befSAndrew Turner double z = a + Shift;
61*f3087befSAndrew Turner uint64_t i = asuint64 (z) - asuint64 (Shift);
62*f3087befSAndrew Turner double r = z - Shift;
63*f3087befSAndrew Turner /* Lookup erf(r) and scale(r) in table.
64*f3087befSAndrew Turner Set erf(r) to 0 and scale to 2/sqrt(pi) for |x| <= 0x1.cp-9. */
65*f3087befSAndrew Turner double erfr = __v_erf_data.tab[i].erf;
66*f3087befSAndrew Turner double scale = __v_erf_data.tab[i].scale;
67*f3087befSAndrew Turner
68*f3087befSAndrew Turner /* erf(x) ~ erf(r) + scale * d * poly (d, r). */
69*f3087befSAndrew Turner double d = a - r;
70*f3087befSAndrew Turner double r2 = r * r;
71*f3087befSAndrew Turner double d2 = d * d;
72*f3087befSAndrew Turner
73*f3087befSAndrew Turner /* poly (d, r) = 1 + p1(r) * d + p2(r) * d^2 + ... + p5(r) * d^5. */
74*f3087befSAndrew Turner double p1 = -r;
75*f3087befSAndrew Turner double p2 = fma (TwoThird, r2, -OneThird);
76*f3087befSAndrew Turner double p3 = -r * fma (OneThird, r2, -0.5);
77*f3087befSAndrew Turner double p4 = fma (fma (TwoOverFifteen, r2, -TwoOverFive), r2, Tenth);
78*f3087befSAndrew Turner double p5
79*f3087befSAndrew Turner = -r * fma (fma (TwoOverFortyFive, r2, -TwoOverNine), r2, Sixth);
80*f3087befSAndrew Turner
81*f3087befSAndrew Turner double p34 = fma (p4, d, p3);
82*f3087befSAndrew Turner double p12 = fma (p2, d, p1);
83*f3087befSAndrew Turner double y = fma (p5, d2, p34);
84*f3087befSAndrew Turner y = fma (y, d2, p12);
85*f3087befSAndrew Turner
86*f3087befSAndrew Turner y = fma (fma (y, d2, d), scale, erfr);
87*f3087befSAndrew Turner return asdouble (asuint64 (y) | sign);
88*f3087befSAndrew Turner }
89*f3087befSAndrew Turner
90*f3087befSAndrew Turner /* Special cases : erf(nan)=nan, erf(+inf)=+1 and erf(-inf)=-1. */
91*f3087befSAndrew Turner if (unlikely (ia >= 0x7ff0000000000000))
92*f3087befSAndrew Turner return (1.0 - (double) (sign >> 62)) + 1.0 / x;
93*f3087befSAndrew Turner
94*f3087befSAndrew Turner /* Boring domain (|x| >= 6.0). */
95*f3087befSAndrew Turner return asdouble (sign | asuint64 (1.0));
96*f3087befSAndrew Turner }
97*f3087befSAndrew Turner
98*f3087befSAndrew Turner TEST_ULP (arm_math_erf, 1.79)
99*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_erf, 0, 5.9921875, 40000)
100*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_erf, 5.9921875, inf, 40000)
101*f3087befSAndrew Turner TEST_SYM_INTERVAL (arm_math_erf, 0, inf, 40000)
102