1 /*
2 * Extended precision inverse error function.
3 *
4 * Copyright (c) 2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7 #define _GNU_SOURCE
8 #include <math.h>
9 #include <stdbool.h>
10 #include <float.h>
11
12 #include "math_config.h"
13 #include "poly_scalar_f64.h"
14
15 #define SQRT_PIl 0x1.c5bf891b4ef6aa79c3b0520d5db9p0l
16 #define HF_SQRT_PIl 0x1.c5bf891b4ef6aa79c3b0520d5db9p-1l
17
18 const static struct
19 {
20 /* We use P_N and Q_N to refer to arrays of coefficients, where P_N is the
21 coeffs of the numerator in table N of Blair et al, and Q_N is the coeffs
22 of the denominator. */
23 double P_17[7], Q_17[7], P_37[8], Q_37[8], P_57[9], Q_57[10];
24 } data = {
25 .P_17 = { 0x1.007ce8f01b2e8p+4, -0x1.6b23cc5c6c6d7p+6, 0x1.74e5f6ceb3548p+7,
26 -0x1.5200bb15cc6bbp+7, 0x1.05d193233a849p+6, -0x1.148c5474ee5e1p+3,
27 0x1.689181bbafd0cp-3 },
28 .Q_17 = { 0x1.d8fb0f913bd7bp+3, -0x1.6d7f25a3f1c24p+6, 0x1.a450d8e7f4cbbp+7,
29 -0x1.bc3480485857p+7, 0x1.ae6b0c504ee02p+6, -0x1.499dfec1a7f5fp+4,
30 0x1p+0 },
31 .P_37 = { -0x1.f3596123109edp-7, 0x1.60b8fe375999ep-2, -0x1.779bb9bef7c0fp+1,
32 0x1.786ea384470a2p+3, -0x1.6a7c1453c85d3p+4, 0x1.31f0fc5613142p+4,
33 -0x1.5ea6c007d4dbbp+2, 0x1.e66f265ce9e5p-3 },
34 .Q_37 = { -0x1.636b2dcf4edbep-7, 0x1.0b5411e2acf29p-2, -0x1.3413109467a0bp+1,
35 0x1.563e8136c554ap+3, -0x1.7b77aab1dcafbp+4, 0x1.8a3e174e05ddcp+4,
36 -0x1.4075c56404eecp+3, 0x1p+0 },
37 .P_57 = { 0x1.b874f9516f7f1p-14, 0x1.5921f2916c1c4p-7, 0x1.145ae7d5b8fa4p-2,
38 0x1.29d6dcc3b2fb7p+1, 0x1.cabe2209a7985p+2, 0x1.11859f0745c4p+3,
39 0x1.b7ec7bc6a2ce5p+2, 0x1.d0419e0bb42aep+1, 0x1.c5aa03eef7258p-1 },
40 .Q_57 = { 0x1.b8747e12691f1p-14, 0x1.59240d8ed1e0ap-7, 0x1.14aef2b181e2p-2,
41 0x1.2cd181bcea52p+1, 0x1.e6e63e0b7aa4cp+2, 0x1.65cf8da94aa3ap+3,
42 0x1.7e5c787b10a36p+3, 0x1.0626d68b6cea3p+3, 0x1.065c5f193abf6p+2,
43 0x1p+0 }
44 };
45
46 /* Inverse error function approximation, based on rational approximation as
47 described in
48 J. M. Blair, C. A. Edwards, and J. H. Johnson,
49 "Rational Chebyshev approximations for the inverse of the error function",
50 Math. Comp. 30, pp. 827--830 (1976).
51 https://doi.org/10.1090/S0025-5718-1976-0421040-7. */
52 static inline double
__erfinv(double x)53 __erfinv (double x)
54 {
55 if (x == 1.0)
56 return __math_oflow (0);
57 if (x == -1.0)
58 return __math_oflow (1);
59
60 double a = fabs (x);
61 if (a > 1)
62 return __math_invalid (x);
63
64 if (a <= 0.75)
65 {
66 double t = x * x - 0.5625;
67 return x * horner_6_f64 (t, data.P_17) / horner_6_f64 (t, data.Q_17);
68 }
69
70 if (a <= 0.9375)
71 {
72 double t = x * x - 0.87890625;
73 return x * horner_7_f64 (t, data.P_37) / horner_7_f64 (t, data.Q_37);
74 }
75
76 double t = 1.0 / (sqrtl (-log1pl (-a)));
77 return horner_8_f64 (t, data.P_57)
78 / (copysign (t, x) * horner_9_f64 (t, data.Q_57));
79 }
80
81 /* Extended-precision variant, which uses the above (or asymptotic estimate) as
82 starting point for Newton refinement. This implementation is a port to C of
83 the version in the SpecialFunctions.jl Julia package, with relaxed stopping
84 criteria for the Newton refinement. */
85 long double
erfinvl(long double x)86 erfinvl (long double x)
87 {
88 if (x == 0)
89 return 0;
90
91 double yf = __erfinv (x);
92 long double y;
93 if (isfinite (yf))
94 y = yf;
95 else
96 {
97 /* Double overflowed, use asymptotic estimate instead. */
98 y = copysignl (sqrtl (-logl (1.0l - fabsl (x)) * SQRT_PIl), x);
99 if (!isfinite (y))
100 return y;
101 }
102
103 double eps = fabs (yf - nextafter (yf, 0));
104 while (true)
105 {
106 long double dy = HF_SQRT_PIl * (erfl (y) - x) * exp (y * y);
107 y -= dy;
108 /* Stopping criterion is different to Julia implementation, but is enough
109 to ensure result is accurate when rounded to double-precision. */
110 if (fabsl (dy) < eps)
111 break;
112 }
113 return y;
114 }
115