1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision vector erfc(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 "sv_math.h"
9*f3087befSAndrew Turner #include "test_sig.h"
10*f3087befSAndrew Turner #include "test_defs.h"
11*f3087befSAndrew Turner
12*f3087befSAndrew Turner static const struct data
13*f3087befSAndrew Turner {
14*f3087befSAndrew Turner uint64_t off_idx, off_arr;
15*f3087befSAndrew Turner double max, shift;
16*f3087befSAndrew Turner double p20, p40, p41, p42;
17*f3087befSAndrew Turner double p51, p52;
18*f3087befSAndrew Turner double q5, r5;
19*f3087befSAndrew Turner double q6, r6;
20*f3087befSAndrew Turner double q7, r7;
21*f3087befSAndrew Turner double q8, r8;
22*f3087befSAndrew Turner double q9, r9;
23*f3087befSAndrew Turner uint64_t table_scale;
24*f3087befSAndrew Turner } data = {
25*f3087befSAndrew Turner /* Set an offset so the range of the index used for lookup is 3487, and it
26*f3087befSAndrew Turner can be clamped using a saturated add on an offset index.
27*f3087befSAndrew Turner Index offset is 0xffffffffffffffff - asuint64(shift) - 3487. */
28*f3087befSAndrew Turner .off_idx = 0xbd3ffffffffff260,
29*f3087befSAndrew Turner .off_arr = 0xfffffffffffff260, /* 0xffffffffffffffff - 3487. */
30*f3087befSAndrew Turner .max = 0x1.b3ep+4, /* 3487/128. */
31*f3087befSAndrew Turner .shift = 0x1p45,
32*f3087befSAndrew Turner .table_scale = 0x37f0000000000000, /* asuint64(0x1p-128). */
33*f3087befSAndrew Turner .p20 = 0x1.5555555555555p-2, /* 1/3, used to compute 2/3 and 1/6. */
34*f3087befSAndrew Turner .p40 = -0x1.999999999999ap-4, /* 1/10. */
35*f3087befSAndrew Turner .p41 = -0x1.999999999999ap-2, /* 2/5. */
36*f3087befSAndrew Turner .p42 = 0x1.1111111111111p-3, /* 2/15. */
37*f3087befSAndrew Turner .p51 = -0x1.c71c71c71c71cp-3, /* 2/9. */
38*f3087befSAndrew Turner .p52 = 0x1.6c16c16c16c17p-5, /* 2/45. */
39*f3087befSAndrew Turner /* Qi = (i+1) / i, for i = 5, ..., 9. */
40*f3087befSAndrew Turner .q5 = 0x1.3333333333333p0,
41*f3087befSAndrew Turner .q6 = 0x1.2aaaaaaaaaaabp0,
42*f3087befSAndrew Turner .q7 = 0x1.2492492492492p0,
43*f3087befSAndrew Turner .q8 = 0x1.2p0,
44*f3087befSAndrew Turner .q9 = 0x1.1c71c71c71c72p0,
45*f3087befSAndrew Turner /* Ri = -2 * i / ((i+1)*(i+2)), for i = 5, ..., 9. */
46*f3087befSAndrew Turner .r5 = -0x1.e79e79e79e79ep-3,
47*f3087befSAndrew Turner .r6 = -0x1.b6db6db6db6dbp-3,
48*f3087befSAndrew Turner .r7 = -0x1.8e38e38e38e39p-3,
49*f3087befSAndrew Turner .r8 = -0x1.6c16c16c16c17p-3,
50*f3087befSAndrew Turner .r9 = -0x1.4f2094f2094f2p-3,
51*f3087befSAndrew Turner };
52*f3087befSAndrew Turner
53*f3087befSAndrew Turner /* Optimized double-precision vector erfc(x).
54*f3087befSAndrew Turner Approximation based on series expansion near x rounded to
55*f3087befSAndrew Turner nearest multiple of 1/128.
56*f3087befSAndrew Turner Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
57*f3087befSAndrew Turner
58*f3087befSAndrew Turner erfc(x) ~ erfc(r) - scale * d * poly(r, d), with
59*f3087befSAndrew Turner
60*f3087befSAndrew Turner poly(r, d) = 1 - r d + (2/3 r^2 - 1/3) d^2 - r (1/3 r^2 - 1/2) d^3
61*f3087befSAndrew Turner + (2/15 r^4 - 2/5 r^2 + 1/10) d^4
62*f3087befSAndrew Turner - r * (2/45 r^4 - 2/9 r^2 + 1/6) d^5
63*f3087befSAndrew Turner + p6(r) d^6 + ... + p10(r) d^10
64*f3087befSAndrew Turner
65*f3087befSAndrew Turner Polynomials p6(r) to p10(r) are computed using recurrence relation
66*f3087befSAndrew Turner
67*f3087befSAndrew Turner 2(i+1)p_i + 2r(i+2)p_{i+1} + (i+2)(i+3)p_{i+2} = 0,
68*f3087befSAndrew Turner with p0 = 1, and p1(r) = -r.
69*f3087befSAndrew Turner
70*f3087befSAndrew Turner Values of erfc(r) and scale are read from lookup tables. Stored values
71*f3087befSAndrew Turner are scaled to avoid hitting the subnormal range.
72*f3087befSAndrew Turner
73*f3087befSAndrew Turner Note that for x < 0, erfc(x) = 2.0 - erfc(-x).
74*f3087befSAndrew Turner
75*f3087befSAndrew Turner Maximum measured error: 1.71 ULP
76*f3087befSAndrew Turner _ZGVsMxv_erfc(0x1.46cfe976733p+4) got 0x1.e15fcbea3e7afp-608
77*f3087befSAndrew Turner want 0x1.e15fcbea3e7adp-608. */
SV_NAME_D1(erfc)78*f3087befSAndrew Turner svfloat64_t SV_NAME_D1 (erfc) (svfloat64_t x, const svbool_t pg)
79*f3087befSAndrew Turner {
80*f3087befSAndrew Turner const struct data *dat = ptr_barrier (&data);
81*f3087befSAndrew Turner
82*f3087befSAndrew Turner svfloat64_t a = svabs_x (pg, x);
83*f3087befSAndrew Turner
84*f3087befSAndrew Turner /* Clamp input at |x| <= 3487/128. */
85*f3087befSAndrew Turner a = svmin_x (pg, a, dat->max);
86*f3087befSAndrew Turner
87*f3087befSAndrew Turner /* Reduce x to the nearest multiple of 1/128. */
88*f3087befSAndrew Turner svfloat64_t shift = sv_f64 (dat->shift);
89*f3087befSAndrew Turner svfloat64_t z = svadd_x (pg, a, shift);
90*f3087befSAndrew Turner
91*f3087befSAndrew Turner /* Saturate index for the NaN case. */
92*f3087befSAndrew Turner svuint64_t i = svqadd (svreinterpret_u64 (z), dat->off_idx);
93*f3087befSAndrew Turner
94*f3087befSAndrew Turner /* Lookup erfc(r) and 2/sqrt(pi)*exp(-r^2) in tables. */
95*f3087befSAndrew Turner i = svadd_x (pg, i, i);
96*f3087befSAndrew Turner const float64_t *p = &__v_erfc_data.tab[0].erfc - 2 * dat->off_arr;
97*f3087befSAndrew Turner svfloat64_t erfcr = svld1_gather_index (pg, p, i);
98*f3087befSAndrew Turner svfloat64_t scale = svld1_gather_index (pg, p + 1, i);
99*f3087befSAndrew Turner
100*f3087befSAndrew Turner /* erfc(x) ~ erfc(r) - scale * d * poly(r, d). */
101*f3087befSAndrew Turner svfloat64_t r = svsub_x (pg, z, shift);
102*f3087befSAndrew Turner svfloat64_t d = svsub_x (pg, a, r);
103*f3087befSAndrew Turner svfloat64_t d2 = svmul_x (pg, d, d);
104*f3087befSAndrew Turner svfloat64_t r2 = svmul_x (pg, r, r);
105*f3087befSAndrew Turner
106*f3087befSAndrew Turner /* poly (d, r) = 1 + p1(r) * d + p2(r) * d^2 + ... + p9(r) * d^9. */
107*f3087befSAndrew Turner svfloat64_t p1 = r;
108*f3087befSAndrew Turner svfloat64_t third = sv_f64 (dat->p20);
109*f3087befSAndrew Turner svfloat64_t twothird = svmul_x (pg, third, 2.0);
110*f3087befSAndrew Turner svfloat64_t sixth = svmul_x (pg, third, 0.5);
111*f3087befSAndrew Turner svfloat64_t p2 = svmls_x (pg, third, r2, twothird);
112*f3087befSAndrew Turner svfloat64_t p3 = svmad_x (pg, r2, third, -0.5);
113*f3087befSAndrew Turner p3 = svmul_x (pg, r, p3);
114*f3087befSAndrew Turner svfloat64_t p4 = svmla_x (pg, sv_f64 (dat->p41), r2, dat->p42);
115*f3087befSAndrew Turner p4 = svmls_x (pg, sv_f64 (dat->p40), r2, p4);
116*f3087befSAndrew Turner svfloat64_t p5 = svmla_x (pg, sv_f64 (dat->p51), r2, dat->p52);
117*f3087befSAndrew Turner p5 = svmla_x (pg, sixth, r2, p5);
118*f3087befSAndrew Turner p5 = svmul_x (pg, r, p5);
119*f3087befSAndrew Turner /* Compute p_i using recurrence relation:
120*f3087befSAndrew Turner p_{i+2} = (p_i + r * Q_{i+1} * p_{i+1}) * R_{i+1}. */
121*f3087befSAndrew Turner svfloat64_t qr5 = svld1rq (svptrue_b64 (), &dat->q5);
122*f3087befSAndrew Turner svfloat64_t qr6 = svld1rq (svptrue_b64 (), &dat->q6);
123*f3087befSAndrew Turner svfloat64_t qr7 = svld1rq (svptrue_b64 (), &dat->q7);
124*f3087befSAndrew Turner svfloat64_t qr8 = svld1rq (svptrue_b64 (), &dat->q8);
125*f3087befSAndrew Turner svfloat64_t qr9 = svld1rq (svptrue_b64 (), &dat->q9);
126*f3087befSAndrew Turner svfloat64_t p6 = svmla_x (pg, p4, p5, svmul_lane (r, qr5, 0));
127*f3087befSAndrew Turner p6 = svmul_lane (p6, qr5, 1);
128*f3087befSAndrew Turner svfloat64_t p7 = svmla_x (pg, p5, p6, svmul_lane (r, qr6, 0));
129*f3087befSAndrew Turner p7 = svmul_lane (p7, qr6, 1);
130*f3087befSAndrew Turner svfloat64_t p8 = svmla_x (pg, p6, p7, svmul_lane (r, qr7, 0));
131*f3087befSAndrew Turner p8 = svmul_lane (p8, qr7, 1);
132*f3087befSAndrew Turner svfloat64_t p9 = svmla_x (pg, p7, p8, svmul_lane (r, qr8, 0));
133*f3087befSAndrew Turner p9 = svmul_lane (p9, qr8, 1);
134*f3087befSAndrew Turner svfloat64_t p10 = svmla_x (pg, p8, p9, svmul_lane (r, qr9, 0));
135*f3087befSAndrew Turner p10 = svmul_lane (p10, qr9, 1);
136*f3087befSAndrew Turner /* Compute polynomial in d using pairwise Horner scheme. */
137*f3087befSAndrew Turner svfloat64_t p90 = svmla_x (pg, p9, d, p10);
138*f3087befSAndrew Turner svfloat64_t p78 = svmla_x (pg, p7, d, p8);
139*f3087befSAndrew Turner svfloat64_t p56 = svmla_x (pg, p5, d, p6);
140*f3087befSAndrew Turner svfloat64_t p34 = svmla_x (pg, p3, d, p4);
141*f3087befSAndrew Turner svfloat64_t p12 = svmla_x (pg, p1, d, p2);
142*f3087befSAndrew Turner svfloat64_t y = svmla_x (pg, p78, d2, p90);
143*f3087befSAndrew Turner y = svmla_x (pg, p56, d2, y);
144*f3087befSAndrew Turner y = svmla_x (pg, p34, d2, y);
145*f3087befSAndrew Turner y = svmla_x (pg, p12, d2, y);
146*f3087befSAndrew Turner
147*f3087befSAndrew Turner y = svmls_x (pg, erfcr, scale, svmls_x (pg, d, d2, y));
148*f3087befSAndrew Turner
149*f3087befSAndrew Turner /* Offset equals 2.0 if sign, else 0.0. */
150*f3087befSAndrew Turner svuint64_t sign = svand_x (pg, svreinterpret_u64 (x), 0x8000000000000000);
151*f3087befSAndrew Turner svfloat64_t off = svreinterpret_f64 (svlsr_x (pg, sign, 1));
152*f3087befSAndrew Turner /* Handle sign and scale back in a single fma. */
153*f3087befSAndrew Turner svfloat64_t fac = svreinterpret_f64 (svorr_x (pg, sign, dat->table_scale));
154*f3087befSAndrew Turner
155*f3087befSAndrew Turner return svmla_x (pg, off, fac, y);
156*f3087befSAndrew Turner }
157*f3087befSAndrew Turner
158*f3087befSAndrew Turner TEST_SIG (SV, D, 1, erfc, -6.0, 28.0)
159*f3087befSAndrew Turner TEST_ULP (SV_NAME_D1 (erfc), 1.21)
160*f3087befSAndrew Turner TEST_DISABLE_FENV (SV_NAME_D1 (erfc))
161*f3087befSAndrew Turner TEST_SYM_INTERVAL (SV_NAME_D1 (erfc), 0.0, 0x1p-26, 40000)
162*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_D1 (erfc), 0x1p-26, 28.0, 40000)
163*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_D1 (erfc), -0x1p-26, -6.0, 40000)
164*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_D1 (erfc), 28.0, inf, 40000)
165*f3087befSAndrew Turner TEST_INTERVAL (SV_NAME_D1 (erfc), 6.0, -inf, 40000)
166*f3087befSAndrew Turner CLOSE_SVE_ATTR
167