1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Single-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_sig.h"
9*f3087befSAndrew Turner #include "test_defs.h"
10*f3087befSAndrew Turner #include "sv_poly_f32.h"
11*f3087befSAndrew Turner #include "sv_logf_inline.h"
12*f3087befSAndrew Turner
13*f3087befSAndrew Turner const static struct data
14*f3087befSAndrew Turner {
15*f3087befSAndrew Turner /* We use P_N and Q_N to refer to arrays of coefficients, where P_N
16*f3087befSAndrew Turner is the coeffs of the numerator in table N of Blair et al, and
17*f3087befSAndrew Turner Q_N is the coeffs of the denominator. Coefficients stored in
18*f3087befSAndrew Turner interleaved format to support lookup scheme. */
19*f3087befSAndrew Turner float P10_2, P29_3, Q10_2, Q29_2;
20*f3087befSAndrew Turner float P10_0, P29_1, P10_1, P29_2;
21*f3087befSAndrew Turner float Q10_0, Q29_0, Q10_1, Q29_1;
22*f3087befSAndrew Turner float P29_0, P_50[6], Q_50[2], tailshift;
23*f3087befSAndrew Turner struct sv_logf_data logf_tbl;
24*f3087befSAndrew Turner } data = { .P10_0 = -0x1.a31268p+3,
25*f3087befSAndrew Turner .P10_1 = 0x1.ac9048p+4,
26*f3087befSAndrew Turner .P10_2 = -0x1.293ff6p+3,
27*f3087befSAndrew Turner .P29_0 = -0x1.fc0252p-4,
28*f3087befSAndrew Turner .P29_1 = 0x1.119d44p+0,
29*f3087befSAndrew Turner .P29_2 = -0x1.f59ee2p+0,
30*f3087befSAndrew Turner .P29_3 = 0x1.b13626p-2,
31*f3087befSAndrew Turner .Q10_0 = -0x1.8265eep+3,
32*f3087befSAndrew Turner .Q10_1 = 0x1.ef5eaep+4,
33*f3087befSAndrew Turner .Q10_2 = -0x1.12665p+4,
34*f3087befSAndrew Turner .Q29_0 = -0x1.69952p-4,
35*f3087befSAndrew Turner .Q29_1 = 0x1.c7b7d2p-1,
36*f3087befSAndrew Turner .Q29_2 = -0x1.167d7p+1,
37*f3087befSAndrew Turner .P_50 = { 0x1.3d8948p-3, 0x1.61f9eap+0, 0x1.61c6bcp-1,
38*f3087befSAndrew Turner -0x1.20c9f2p+0, 0x1.5c704cp-1, -0x1.50c6bep-3 },
39*f3087befSAndrew Turner .Q_50 = { 0x1.3d7dacp-3, 0x1.629e5p+0 },
40*f3087befSAndrew Turner .tailshift = -0.87890625,
41*f3087befSAndrew Turner .logf_tbl = SV_LOGF_CONSTANTS };
42*f3087befSAndrew Turner
43*f3087befSAndrew Turner static inline svfloat32_t
special(svbool_t pg,svfloat32_t x,const struct data * d)44*f3087befSAndrew Turner special (svbool_t pg, svfloat32_t x, const struct data *d)
45*f3087befSAndrew Turner {
46*f3087befSAndrew Turner svfloat32_t ax = svabs_x (pg, x);
47*f3087befSAndrew Turner svfloat32_t t = svdivr_x (
48*f3087befSAndrew Turner pg,
49*f3087befSAndrew Turner svsqrt_x (pg, svneg_x (pg, sv_logf_inline (pg, svsubr_x (pg, ax, 1),
50*f3087befSAndrew Turner &d->logf_tbl))),
51*f3087befSAndrew Turner 1);
52*f3087befSAndrew Turner svuint32_t sign
53*f3087befSAndrew Turner = sveor_x (pg, svreinterpret_u32 (ax), svreinterpret_u32 (x));
54*f3087befSAndrew Turner svfloat32_t ts
55*f3087befSAndrew Turner = svreinterpret_f32 (svorr_x (pg, sign, svreinterpret_u32 (t)));
56*f3087befSAndrew Turner svfloat32_t q
57*f3087befSAndrew Turner = svmla_x (pg, sv_f32 (d->Q_50[0]), svadd_x (pg, t, d->Q_50[1]), t);
58*f3087befSAndrew Turner return svdiv_x (pg, sv_horner_5_f32_x (pg, t, d->P_50), svmul_x (pg, ts, q));
59*f3087befSAndrew Turner }
60*f3087befSAndrew Turner
61*f3087befSAndrew Turner static inline svfloat32_t
notails(svbool_t pg,svfloat32_t x,const struct data * d)62*f3087befSAndrew Turner notails (svbool_t pg, svfloat32_t x, const struct data *d)
63*f3087befSAndrew Turner {
64*f3087befSAndrew Turner /* Shortcut when no input is in a tail region - no need to gather shift or
65*f3087befSAndrew Turner coefficients. */
66*f3087befSAndrew Turner svfloat32_t t = svmad_x (pg, x, x, -0.5625);
67*f3087befSAndrew Turner svfloat32_t q = svadd_x (pg, t, d->Q10_2);
68*f3087befSAndrew Turner q = svmad_x (pg, t, q, d->Q10_1);
69*f3087befSAndrew Turner q = svmad_x (pg, t, q, d->Q10_0);
70*f3087befSAndrew Turner
71*f3087befSAndrew Turner svfloat32_t p = svmla_x (pg, sv_f32 (d->P10_1), t, d->P10_2);
72*f3087befSAndrew Turner p = svmad_x (pg, p, t, d->P10_0);
73*f3087befSAndrew Turner
74*f3087befSAndrew Turner return svdiv_x (pg, svmul_x (pg, x, p), q);
75*f3087befSAndrew Turner }
76*f3087befSAndrew Turner
77*f3087befSAndrew Turner /* Vector implementation of Blair et al's rational approximation to inverse
78*f3087befSAndrew Turner error function in single-precision. Worst-case error is 4.71 ULP, in the
79*f3087befSAndrew Turner tail region:
80*f3087befSAndrew Turner _ZGVsMxv_erfinvf(0x1.f84e9ap-1) got 0x1.b8326ap+0
81*f3087befSAndrew Turner want 0x1.b83274p+0. */
SV_NAME_F1(erfinv)82*f3087befSAndrew Turner svfloat32_t SV_NAME_F1 (erfinv) (svfloat32_t x, svbool_t pg)
83*f3087befSAndrew Turner {
84*f3087befSAndrew Turner const struct data *d = ptr_barrier (&data);
85*f3087befSAndrew Turner
86*f3087befSAndrew Turner /* Calculate inverse error using algorithm described in
87*f3087befSAndrew Turner J. M. Blair, C. A. Edwards, and J. H. Johnson,
88*f3087befSAndrew Turner "Rational Chebyshev approximations for the inverse of the error function",
89*f3087befSAndrew Turner Math. Comp. 30, pp. 827--830 (1976).
90*f3087befSAndrew Turner https://doi.org/10.1090/S0025-5718-1976-0421040-7. */
91*f3087befSAndrew Turner
92*f3087befSAndrew Turner /* Algorithm has 3 intervals:
93*f3087befSAndrew Turner - 'Normal' region [-0.75, 0.75]
94*f3087befSAndrew Turner - Tail region [0.75, 0.9375] U [-0.9375, -0.75]
95*f3087befSAndrew Turner - Extreme tail [-1, -0.9375] U [0.9375, 1]
96*f3087befSAndrew Turner Normal and tail are both rational approximation of similar order on
97*f3087befSAndrew Turner shifted input - these are typically performed in parallel using gather
98*f3087befSAndrew Turner loads to obtain correct coefficients depending on interval. */
99*f3087befSAndrew Turner svbool_t is_tail = svacge (pg, x, 0.75);
100*f3087befSAndrew Turner svbool_t extreme_tail = svacge (pg, x, 0.9375);
101*f3087befSAndrew Turner
102*f3087befSAndrew Turner if (likely (!svptest_any (pg, is_tail)))
103*f3087befSAndrew Turner return notails (pg, x, d);
104*f3087befSAndrew Turner
105*f3087befSAndrew Turner /* Select requisite shift depending on interval: polynomial is evaluated on
106*f3087befSAndrew Turner x * x - shift.
107*f3087befSAndrew Turner Normal shift = 0.5625
108*f3087befSAndrew Turner Tail shift = 0.87890625. */
109*f3087befSAndrew Turner svfloat32_t t = svmla_x (
110*f3087befSAndrew Turner pg, svsel (is_tail, sv_f32 (d->tailshift), sv_f32 (-0.5625)), x, x);
111*f3087befSAndrew Turner
112*f3087befSAndrew Turner svuint32_t idx = svdup_u32_z (is_tail, 1);
113*f3087befSAndrew Turner svuint32_t idxhi = svadd_x (pg, idx, 2);
114*f3087befSAndrew Turner
115*f3087befSAndrew Turner /* Load coeffs in quadwords and select them according to interval. */
116*f3087befSAndrew Turner svfloat32_t pqhi = svld1rq (svptrue_b32 (), &d->P10_2);
117*f3087befSAndrew Turner svfloat32_t plo = svld1rq (svptrue_b32 (), &d->P10_0);
118*f3087befSAndrew Turner svfloat32_t qlo = svld1rq (svptrue_b32 (), &d->Q10_0);
119*f3087befSAndrew Turner
120*f3087befSAndrew Turner svfloat32_t p2 = svtbl (pqhi, idx);
121*f3087befSAndrew Turner svfloat32_t p1 = svtbl (plo, idxhi);
122*f3087befSAndrew Turner svfloat32_t p0 = svtbl (plo, idx);
123*f3087befSAndrew Turner svfloat32_t q0 = svtbl (qlo, idx);
124*f3087befSAndrew Turner svfloat32_t q1 = svtbl (qlo, idxhi);
125*f3087befSAndrew Turner svfloat32_t q2 = svtbl (pqhi, idxhi);
126*f3087befSAndrew Turner
127*f3087befSAndrew Turner svfloat32_t p = svmla_x (pg, p1, p2, t);
128*f3087befSAndrew Turner p = svmla_x (pg, p0, p, t);
129*f3087befSAndrew Turner /* Tail polynomial has higher order - merge with normal lanes. */
130*f3087befSAndrew Turner p = svmad_m (is_tail, p, t, d->P29_0);
131*f3087befSAndrew Turner svfloat32_t y = svmul_x (pg, x, p);
132*f3087befSAndrew Turner
133*f3087befSAndrew Turner /* Least significant term of both Q polynomials is 1, so no need to generate
134*f3087befSAndrew Turner it. */
135*f3087befSAndrew Turner svfloat32_t q = svadd_x (pg, t, q2);
136*f3087befSAndrew Turner q = svmla_x (pg, q1, q, t);
137*f3087befSAndrew Turner q = svmla_x (pg, q0, q, t);
138*f3087befSAndrew Turner
139*f3087befSAndrew Turner if (unlikely (svptest_any (pg, extreme_tail)))
140*f3087befSAndrew Turner return svsel (extreme_tail, special (extreme_tail, x, d),
141*f3087befSAndrew Turner svdiv_x (pg, y, q));
142*f3087befSAndrew Turner return svdiv_x (pg, y, q);
143*f3087befSAndrew Turner }
144*f3087befSAndrew Turner
145*f3087befSAndrew Turner #if USE_MPFR
146*f3087befSAndrew Turner # warning Not generating tests for _ZGVsMxv_erfinvf, as MPFR has no suitable reference
147*f3087befSAndrew Turner #else
148*f3087befSAndrew Turner TEST_SIG (SV, F, 1, erfinv, -0.99, 0.99)
149*f3087befSAndrew Turner TEST_ULP (SV_NAME_F1 (erfinv), 4.09)
150*f3087befSAndrew Turner TEST_DISABLE_FENV (SV_NAME_F1 (erfinv))
151*f3087befSAndrew Turner TEST_SYM_INTERVAL (SV_NAME_F1 (erfinv), 0, 1, 40000)
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