1 /*
2 * Single-precision erf(x) function.
3 *
4 * Copyright (c) 2023-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "math_config.h"
9 #include "test_sig.h"
10 #include "test_defs.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
arm_math_erff(float x)40 arm_math_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 = __v_erff_data.tab[i].erf;
60 float scale = __v_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 TEST_ULP (arm_math_erff, 1.43)
79 TEST_SYM_INTERVAL (arm_math_erff, 0, 3.9375, 40000)
80 TEST_SYM_INTERVAL (arm_math_erff, 3.9375, inf, 40000)
81 TEST_SYM_INTERVAL (arm_math_erff, 0, inf, 40000)
82