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 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