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