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