1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision inverse error function (SVE variant).
3*f3087befSAndrew Turner *
4*f3087befSAndrew Turner * Copyright (c) 2024, Arm Limited.
5*f3087befSAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*f3087befSAndrew Turner */
7*f3087befSAndrew Turner #include "sv_math.h"
8*f3087befSAndrew Turner #include "test_defs.h"
9*f3087befSAndrew Turner #include "math_config.h"
10*f3087befSAndrew Turner #include "test_sig.h"
11*f3087befSAndrew Turner #include "sv_poly_f64.h"
12*f3087befSAndrew Turner #define SV_LOG_INLINE_POLY_ORDER 4
13*f3087befSAndrew Turner #include "sv_log_inline.h"
14*f3087befSAndrew Turner
15*f3087befSAndrew Turner const static struct data
16*f3087befSAndrew Turner {
17*f3087befSAndrew Turner /* We use P_N and Q_N to refer to arrays of coefficients, where P_N is the
18*f3087befSAndrew Turner coeffs of the numerator in table N of Blair et al, and Q_N is the coeffs
19*f3087befSAndrew Turner of the denominator. P is interleaved P_17 and P_37, similar for Q. */
20*f3087befSAndrew Turner double P[7][2], Q[7][2];
21*f3087befSAndrew Turner double P_57[9], Q_57[9], tailshift, P37_0;
22*f3087befSAndrew Turner struct sv_log_inline_data log_tbl;
23*f3087befSAndrew Turner } data = {
24*f3087befSAndrew Turner .P37_0 = -0x1.f3596123109edp-7,
25*f3087befSAndrew Turner .tailshift = -0.87890625,
26*f3087befSAndrew Turner .P = { { 0x1.007ce8f01b2e8p+4, 0x1.60b8fe375999ep-2 },
27*f3087befSAndrew Turner { -0x1.6b23cc5c6c6d7p+6, -0x1.779bb9bef7c0fp+1 },
28*f3087befSAndrew Turner { 0x1.74e5f6ceb3548p+7, 0x1.786ea384470a2p+3 },
29*f3087befSAndrew Turner { -0x1.5200bb15cc6bbp+7, -0x1.6a7c1453c85d3p+4 },
30*f3087befSAndrew Turner { 0x1.05d193233a849p+6, 0x1.31f0fc5613142p+4 },
31*f3087befSAndrew Turner { -0x1.148c5474ee5e1p+3, -0x1.5ea6c007d4dbbp+2 },
32*f3087befSAndrew Turner { 0x1.689181bbafd0cp-3, 0x1.e66f265ce9e5p-3 } },
33*f3087befSAndrew Turner .Q = { { 0x1.d8fb0f913bd7bp+3, -0x1.636b2dcf4edbep-7 },
34*f3087befSAndrew Turner { -0x1.6d7f25a3f1c24p+6, 0x1.0b5411e2acf29p-2 },
35*f3087befSAndrew Turner { 0x1.a450d8e7f4cbbp+7, -0x1.3413109467a0bp+1 },
36*f3087befSAndrew Turner { -0x1.bc3480485857p+7, 0x1.563e8136c554ap+3 },
37*f3087befSAndrew Turner { 0x1.ae6b0c504ee02p+6, -0x1.7b77aab1dcafbp+4 },
38*f3087befSAndrew Turner { -0x1.499dfec1a7f5fp+4, 0x1.8a3e174e05ddcp+4 },
39*f3087befSAndrew Turner { 0x1p+0, -0x1.4075c56404eecp+3 } },
40*f3087befSAndrew Turner .P_57 = { 0x1.b874f9516f7f1p-14, 0x1.5921f2916c1c4p-7, 0x1.145ae7d5b8fa4p-2,
41*f3087befSAndrew Turner 0x1.29d6dcc3b2fb7p+1, 0x1.cabe2209a7985p+2, 0x1.11859f0745c4p+3,
42*f3087befSAndrew Turner 0x1.b7ec7bc6a2ce5p+2, 0x1.d0419e0bb42aep+1, 0x1.c5aa03eef7258p-1 },
43*f3087befSAndrew Turner .Q_57 = { 0x1.b8747e12691f1p-14, 0x1.59240d8ed1e0ap-7, 0x1.14aef2b181e2p-2,
44*f3087befSAndrew Turner 0x1.2cd181bcea52p+1, 0x1.e6e63e0b7aa4cp+2, 0x1.65cf8da94aa3ap+3,
45*f3087befSAndrew Turner 0x1.7e5c787b10a36p+3, 0x1.0626d68b6cea3p+3, 0x1.065c5f193abf6p+2 },
46*f3087befSAndrew Turner .log_tbl = SV_LOG_CONSTANTS
47*f3087befSAndrew Turner };
48*f3087befSAndrew Turner
49*f3087befSAndrew Turner static inline svfloat64_t
special(svbool_t pg,svfloat64_t x,const struct data * d)50*f3087befSAndrew Turner special (svbool_t pg, svfloat64_t x, const struct data *d)
51*f3087befSAndrew Turner {
52*f3087befSAndrew Turner /* Note erfinv(inf) should return NaN, and erfinv(1) should return Inf.
53*f3087befSAndrew Turner By using log here, instead of log1p, we return finite values for both
54*f3087befSAndrew Turner these inputs, and values outside [-1, 1]. This is non-compliant, but is an
55*f3087befSAndrew Turner acceptable optimisation at Ofast. To get correct behaviour for all finite
56*f3087befSAndrew Turner values use the log1p_inline helper on -abs(x) - note that erfinv(inf)
57*f3087befSAndrew Turner will still be finite. */
58*f3087befSAndrew Turner svfloat64_t ax = svabs_x (pg, x);
59*f3087befSAndrew Turner svfloat64_t t
60*f3087befSAndrew Turner = svneg_x (pg, sv_log_inline (pg, svsubr_x (pg, ax, 1), &d->log_tbl));
61*f3087befSAndrew Turner t = svdivr_x (pg, svsqrt_x (pg, t), 1);
62*f3087befSAndrew Turner svuint64_t sign
63*f3087befSAndrew Turner = sveor_x (pg, svreinterpret_u64 (ax), svreinterpret_u64 (x));
64*f3087befSAndrew Turner svfloat64_t ts
65*f3087befSAndrew Turner = svreinterpret_f64 (svorr_x (pg, sign, svreinterpret_u64 (t)));
66*f3087befSAndrew Turner
67*f3087befSAndrew Turner svfloat64_t q = svadd_x (pg, t, d->Q_57[8]);
68*f3087befSAndrew Turner for (int i = 7; i >= 0; i--)
69*f3087befSAndrew Turner q = svmad_x (pg, q, t, d->Q_57[i]);
70*f3087befSAndrew Turner
71*f3087befSAndrew Turner return svdiv_x (pg, sv_horner_8_f64_x (pg, t, d->P_57), svmul_x (pg, ts, q));
72*f3087befSAndrew Turner }
73*f3087befSAndrew Turner
74*f3087befSAndrew Turner static inline svfloat64_t
lookup(const double * c,svuint64_t idx)75*f3087befSAndrew Turner lookup (const double *c, svuint64_t idx)
76*f3087befSAndrew Turner {
77*f3087befSAndrew Turner svfloat64_t x = svld1rq_f64 (svptrue_b64 (), c);
78*f3087befSAndrew Turner return svtbl (x, idx);
79*f3087befSAndrew Turner }
80*f3087befSAndrew Turner
81*f3087befSAndrew Turner static inline svfloat64_t
notails(svbool_t pg,svfloat64_t x,const struct data * d)82*f3087befSAndrew Turner notails (svbool_t pg, svfloat64_t x, const struct data *d)
83*f3087befSAndrew Turner {
84*f3087befSAndrew Turner svfloat64_t t = svmad_x (pg, x, x, -0.5625);
85*f3087befSAndrew Turner svfloat64_t p = svmla_x (pg, sv_f64 (d->P[5][0]), t, d->P[6][0]);
86*f3087befSAndrew Turner svfloat64_t q = svadd_x (pg, t, d->Q[5][0]);
87*f3087befSAndrew Turner for (int i = 4; i >= 0; i--)
88*f3087befSAndrew Turner {
89*f3087befSAndrew Turner p = svmad_x (pg, t, p, d->P[i][0]);
90*f3087befSAndrew Turner q = svmad_x (pg, t, q, d->Q[i][0]);
91*f3087befSAndrew Turner }
92*f3087befSAndrew Turner p = svmul_x (pg, p, x);
93*f3087befSAndrew Turner return svdiv_x (pg, p, q);
94*f3087befSAndrew Turner }
95*f3087befSAndrew Turner
96*f3087befSAndrew Turner /* Vector implementation of Blair et al's rational approximation to inverse
97*f3087befSAndrew Turner error function in double precision. Largest observed error is 24.75 ULP:
98*f3087befSAndrew Turner _ZGVsMxv_erfinv(0x1.fc861d81c2ba8p-1) got 0x1.ea05472686625p+0
99*f3087befSAndrew Turner want 0x1.ea0547268660cp+0. */
SV_NAME_D1(erfinv)100*f3087befSAndrew Turner svfloat64_t SV_NAME_D1 (erfinv) (svfloat64_t x, svbool_t pg)
101*f3087befSAndrew Turner {
102*f3087befSAndrew Turner const struct data *d = ptr_barrier (&data);
103*f3087befSAndrew Turner /* Calculate inverse error using algorithm described in
104*f3087befSAndrew Turner J. M. Blair, C. A. Edwards, and J. H. Johnson,
105*f3087befSAndrew Turner "Rational Chebyshev approximations for the inverse of the error function",
106*f3087befSAndrew Turner Math. Comp. 30, pp. 827--830 (1976).
107*f3087befSAndrew Turner https://doi.org/10.1090/S0025-5718-1976-0421040-7.
108*f3087befSAndrew Turner
109*f3087befSAndrew Turner Algorithm has 3 intervals:
110*f3087befSAndrew Turner - 'Normal' region [-0.75, 0.75]
111*f3087befSAndrew Turner - Tail region [0.75, 0.9375] U [-0.9375, -0.75]
112*f3087befSAndrew Turner - Extreme tail [-1, -0.9375] U [0.9375, 1]
113*f3087befSAndrew Turner Normal and tail are both rational approximation of similar order on
114*f3087befSAndrew Turner shifted input - these are typically performed in parallel using gather
115*f3087befSAndrew Turner loads to obtain correct coefficients depending on interval. */
116*f3087befSAndrew Turner
117*f3087befSAndrew Turner svbool_t no_tail = svacle (pg, x, 0.75);
118*f3087befSAndrew Turner if (unlikely (!svptest_any (pg, svnot_z (pg, no_tail))))
119*f3087befSAndrew Turner return notails (pg, x, d);
120*f3087befSAndrew Turner
121*f3087befSAndrew Turner svbool_t is_tail = svnot_z (pg, no_tail);
122*f3087befSAndrew Turner svbool_t extreme_tail = svacgt (pg, x, 0.9375);
123*f3087befSAndrew Turner svuint64_t idx = svdup_n_u64_z (is_tail, 1);
124*f3087befSAndrew Turner
125*f3087befSAndrew Turner svfloat64_t t = svsel_f64 (is_tail, sv_f64 (d->tailshift), sv_f64 (-0.5625));
126*f3087befSAndrew Turner t = svmla_x (pg, t, x, x);
127*f3087befSAndrew Turner
128*f3087befSAndrew Turner svfloat64_t p = lookup (&d->P[6][0], idx);
129*f3087befSAndrew Turner svfloat64_t q
130*f3087befSAndrew Turner = svmla_x (pg, lookup (&d->Q[6][0], idx), svdup_n_f64_z (is_tail, 1), t);
131*f3087befSAndrew Turner for (int i = 5; i >= 0; i--)
132*f3087befSAndrew Turner {
133*f3087befSAndrew Turner p = svmla_x (pg, lookup (&d->P[i][0], idx), p, t);
134*f3087befSAndrew Turner q = svmla_x (pg, lookup (&d->Q[i][0], idx), q, t);
135*f3087befSAndrew Turner }
136*f3087befSAndrew Turner p = svmad_m (is_tail, p, t, d->P37_0);
137*f3087befSAndrew Turner p = svmul_x (pg, p, x);
138*f3087befSAndrew Turner
139*f3087befSAndrew Turner if (likely (svptest_any (pg, extreme_tail)))
140*f3087befSAndrew Turner return svsel (extreme_tail, special (pg, x, d), svdiv_x (pg, p, q));
141*f3087befSAndrew Turner return svdiv_x (pg, p, q);
142*f3087befSAndrew Turner }
143*f3087befSAndrew Turner
144*f3087befSAndrew Turner #if USE_MPFR
145*f3087befSAndrew Turner # warning Not generating tests for _ZGVsMxv_erfinv, as MPFR has no suitable reference
146*f3087befSAndrew Turner #else
147*f3087befSAndrew Turner TEST_SIG (SV, D, 1, erfinv, -0.99, 0.99)
148*f3087befSAndrew Turner TEST_ULP (SV_NAME_D1 (erfinv), 24.5)
149*f3087befSAndrew Turner TEST_DISABLE_FENV (SV_NAME_D1 (erfinv))
150*f3087befSAndrew Turner /* Test with control lane in each interval. */
151*f3087befSAndrew Turner TEST_SYM_INTERVAL (SV_NAME_F1 (erfinv), 0, 1, 100000)
152*f3087befSAndrew Turner TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.5)
153*f3087befSAndrew Turner TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.8)
154*f3087befSAndrew Turner TEST_CONTROL_VALUE (SV_NAME_F1 (erfinv), 0.95)
155*f3087befSAndrew Turner #endif
156*f3087befSAndrew Turner CLOSE_SVE_ATTR
157