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