2  * Double-precision vector erfc(x) function.
4  * Copyright (c) 2023-2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
14   uint64x2_t offset, table_scale;  member
23   /* Set an offset so the range of the index used for lookup is 3487, and it
24      can be clamped using a saturated add on an offset index.
25      Index offset is 0xffffffffffffffff - asuint64(shift) - 3487.  */
26   .offset = V2 (0xbd3ffffffffff260),
27   .table_scale = V2 (0x37f0000000000000 << 1), /* asuint64 (2^-128) << 1.  */
30   .p20 = V2 (0x1.5555555555555p-2),  /* 1/3, used to compute 2/3 and 1/6.  */
31   .p40 = V2 (-0x1.999999999999ap-4), /* 1/10.  */
32   .p41 = V2 (-0x1.999999999999ap-2), /* 2/5.  */
33   .p42 = 0x1.1111111111111p-3,	     /* 2/15.  */
34   .p51 = V2 (-0x1.c71c71c71c71cp-3), /* 2/9.  */
35   .p52 = 0x1.6c16c16c16c17p-5,	     /* 2/45.  */
36   /* Qi = (i+1) / i, Ri = -2 * i / ((i+1)*(i+2)), for i = 5, ..., 9.  */
37   .qr5 = { 0x1.3333333333333p0, -0x1.e79e79e79e79ep-3 },
38   .qr6 = { 0x1.2aaaaaaaaaaabp0, -0x1.b6db6db6db6dbp-3 },
39   .qr7 = { 0x1.2492492492492p0, -0x1.8e38e38e38e39p-3 },
40   .qr8 = { 0x1.2p0, -0x1.6c16c16c16c17p-3 },
41   .qr9 = { 0x1.1c71c71c71c72p0, -0x1.4f2094f2094f2p-3 },
47 #define TinyBound 0x4000000000000000 /* 0x1p-511 << 1.  */
48 #define Off 0xfffffffffffff260	     /* 0xffffffffffffffff - 3487.  */
61       = vld1q_f64 (&__v_erfc_data.tab[vgetq_lane_u64 (i, 0) - Off].erfc);  in lookup()
63       = vld1q_f64 (&__v_erfc_data.tab[vgetq_lane_u64 (i, 1) - Off].erfc);  in lookup()
71 special_case (float64x2_t x, float64x2_t y, uint64x2_t cmp)  in special_case()  argument
73   return v_call_f64 (erfc, x, y, cmp);  in special_case()
77 /* Optimized double-precision vector erfc(x).
81    Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
83    erfc(x) ~ erfc(r) - scale * d * poly(r, d), with
85    poly(r, d) = 1 - r d + (2/3 r^2 - 1/3) d^2 - r (1/3 r^2 - 1/2) d^3
86 		+ (2/15 r^4 - 2/5 r^2 + 1/10) d^4
87 		- r * (2/45 r^4 - 2/9 r^2 + 1/6) d^5
93    with p0 = 1, and p1(r) = -r.
98    Note that for x < 0, erfc(x) = 2.0 - erfc(-x).
101    V_NAME_D1 (erfc)(0x1.46cfe976733p+4) got 0x1.e15fcbea3e7afp-608
102 				       want 0x1.e15fcbea3e7adp-608.  */
109   /* |x| < 2^-511. Avoid fabs by left-shifting by 1.  */  in V_NAME_D1()
115 				vreinterpretq_s64_f64 (dat->uflow_bound));  in V_NAME_D1()
126   a = vminq_f64 (a, dat->max);  in V_NAME_D1()
130   float64x2_t shift = dat->shift;  in V_NAME_D1()
134      min. Instead we offset the table address and the index, then use a  in V_NAME_D1()
136   uint64x2_t i = vqaddq_u64 (vreinterpretq_u64_f64 (z), dat->offset);  in V_NAME_D1()
140   /* erfc(x) ~ erfc(r) - scale * d * poly(r, d).  */  in V_NAME_D1()
147   float64x2_t p2 = vfmsq_f64 (dat->p20, r2, vaddq_f64 (dat->p20, dat->p20));  in V_NAME_D1()
148   float64x2_t p3 = vmulq_f64 (r, vfmaq_f64 (v_f64 (-0.5), r2, dat->p20));  in V_NAME_D1()
149   float64x2_t p42_p52 = vld1q_f64 (&dat->p42);  in V_NAME_D1()
150   float64x2_t p4 = vfmaq_laneq_f64 (dat->p41, r2, p42_p52, 0);  in V_NAME_D1()
151   p4 = vfmsq_f64 (dat->p40, r2, p4);  in V_NAME_D1()
152   float64x2_t p5 = vfmaq_laneq_f64 (dat->p51, r2, p42_p52, 1);  in V_NAME_D1()
153   p5 = vmulq_f64 (r, vfmaq_f64 (vmulq_f64 (v_f64 (0.5), dat->p20), r2, p5));  in V_NAME_D1()
156   float64x2_t qr5 = vld1q_f64 (dat->qr5), qr6 = vld1q_f64 (dat->qr6),  in V_NAME_D1()
157 	      qr7 = vld1q_f64 (dat->qr7), qr8 = vld1q_f64 (dat->qr8),  in V_NAME_D1()
158 	      qr9 = vld1q_f64 (dat->qr9);  in V_NAME_D1()
175   float64x2_t y = vfmaq_f64 (p78, d2, p90);  in V_NAME_D1()  local
176   y = vfmaq_f64 (p56, d2, y);  in V_NAME_D1()
177   y = vfmaq_f64 (p34, d2, y);  in V_NAME_D1()
178   y = vfmaq_f64 (p12, d2, y);  in V_NAME_D1()
180   y = vfmsq_f64 (e.erfc, e.scale, vfmsq_f64 (d, d2, y));  in V_NAME_D1()
182   /* Offset equals 2.0 if sign, else 0.0.  */  in V_NAME_D1()
188       vsraq_n_u64 (vshlq_n_u64 (sign, 63), dat->table_scale, 1));  in V_NAME_D1()
192     return special_case (xm, vfmaq_f64 (off, fac, y), cmp);  in V_NAME_D1()
195   return vfmaq_f64 (off, fac, y);  in V_NAME_D1()
198 TEST_SIG (V, D, 1, erfc, -6.0, 28.0)
201 TEST_SYM_INTERVAL (V_NAME_D1 (erfc), 0, 0x1p-26, 40000)
202 TEST_INTERVAL (V_NAME_D1 (erfc), 0x1p-26, 28.0, 40000)
203 TEST_INTERVAL (V_NAME_D1 (erfc), -0x1p-26, -6.0, 40000)
205 TEST_INTERVAL (V_NAME_D1 (erfc), -6.0, -inf, 40000)