1/*
2 * Double-precision vector erfc(x) function.
3 *
4 * Copyright (c) 2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8#include "v_math.h"
9#include "pl_sig.h"
10#include "pl_test.h"
11
12static const struct data
13{
14  uint64x2_t offset, table_scale;
15  float64x2_t max, shift;
16  float64x2_t p20, p40, p41, p42;
17  float64x2_t p51, p52;
18  float64x2_t qr5, qr6, qr7, qr8, qr9;
19#if WANT_SIMD_EXCEPT
20  float64x2_t uflow_bound;
21#endif
22} data = {
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.  */
28  .max = V2 (0x1.b3ep+4),		       /* 3487/128.  */
29  .shift = V2 (0x1p45),
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 = V2 (0x1.1111111111111p-3),  /* 2/15.  */
34  .p51 = V2 (-0x1.c71c71c71c71cp-3), /* 2/9.  */
35  .p52 = V2 (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 },
42#if WANT_SIMD_EXCEPT
43  .uflow_bound = V2 (0x1.a8b12fc6e4892p+4),
44#endif
45};
46
47#define TinyBound 0x4000000000000000 /* 0x1p-511 << 1.  */
48#define Off 0xfffffffffffff260	     /* 0xffffffffffffffff - 3487.  */
49
50struct entry
51{
52  float64x2_t erfc;
53  float64x2_t scale;
54};
55
56static inline struct entry
57lookup (uint64x2_t i)
58{
59  struct entry e;
60  float64x2_t e1 = vld1q_f64 ((float64_t *) (__erfc_data.tab - Off + i[0])),
61	      e2 = vld1q_f64 ((float64_t *) (__erfc_data.tab - Off + i[1]));
62  e.erfc = vuzp1q_f64 (e1, e2);
63  e.scale = vuzp2q_f64 (e1, e2);
64  return e;
65}
66
67#if WANT_SIMD_EXCEPT
68static float64x2_t VPCS_ATTR NOINLINE
69special_case (float64x2_t x, float64x2_t y, uint64x2_t cmp)
70{
71  return v_call_f64 (erfc, x, y, cmp);
72}
73#endif
74
75/* Optimized double-precision vector erfc(x).
76   Approximation based on series expansion near x rounded to
77   nearest multiple of 1/128.
78
79   Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
80
81   erfc(x) ~ erfc(r) - scale * d * poly(r, d), with
82
83   poly(r, d) = 1 - r d + (2/3 r^2 - 1/3) d^2 - r (1/3 r^2 - 1/2) d^3
84		+ (2/15 r^4 - 2/5 r^2 + 1/10) d^4
85		- r * (2/45 r^4 - 2/9 r^2 + 1/6) d^5
86		+ p6(r) d^6 + ... + p10(r) d^10
87
88   Polynomials p6(r) to p10(r) are computed using recurrence relation
89
90   2(i+1)p_i + 2r(i+2)p_{i+1} + (i+2)(i+3)p_{i+2} = 0,
91   with p0 = 1, and p1(r) = -r.
92
93   Values of erfc(r) and scale are read from lookup tables. Stored values
94   are scaled to avoid hitting the subnormal range.
95
96   Note that for x < 0, erfc(x) = 2.0 - erfc(-x).
97
98   Maximum measured error: 1.71 ULP
99   V_NAME_D1 (erfc)(0x1.46cfe976733p+4) got 0x1.e15fcbea3e7afp-608
100				       want 0x1.e15fcbea3e7adp-608.  */
101VPCS_ATTR
102float64x2_t V_NAME_D1 (erfc) (float64x2_t x)
103{
104  const struct data *dat = ptr_barrier (&data);
105
106#if WANT_SIMD_EXCEPT
107  /* |x| < 2^-511. Avoid fabs by left-shifting by 1.  */
108  uint64x2_t ix = vreinterpretq_u64_f64 (x);
109  uint64x2_t cmp = vcltq_u64 (vaddq_u64 (ix, ix), v_u64 (TinyBound));
110  /* x >= ~26.54 (into subnormal case and uflow case). Comparison is done in
111     integer domain to avoid raising exceptions in presence of nans.  */
112  uint64x2_t uflow = vcgeq_s64 (vreinterpretq_s64_f64 (x),
113				vreinterpretq_s64_f64 (dat->uflow_bound));
114  cmp = vorrq_u64 (cmp, uflow);
115  float64x2_t xm = x;
116  /* If any lanes are special, mask them with 0 and retain a copy of x to allow
117     special case handler to fix special lanes later. This is only necessary if
118     fenv exceptions are to be triggered correctly.  */
119  if (unlikely (v_any_u64 (cmp)))
120    x = v_zerofy_f64 (x, cmp);
121#endif
122
123  float64x2_t a = vabsq_f64 (x);
124  a = vminq_f64 (a, dat->max);
125
126  /* Lookup erfc(r) and scale(r) in tables, e.g. set erfc(r) to 0 and scale to
127     2/sqrt(pi), when x reduced to r = 0.  */
128  float64x2_t shift = dat->shift;
129  float64x2_t z = vaddq_f64 (a, shift);
130
131  /* Clamp index to a range of 3487. A naive approach would use a subtract and
132     min. Instead we offset the table address and the index, then use a
133     saturating add.  */
134  uint64x2_t i = vqaddq_u64 (vreinterpretq_u64_f64 (z), dat->offset);
135
136  struct entry e = lookup (i);
137
138  /* erfc(x) ~ erfc(r) - scale * d * poly(r, d).  */
139  float64x2_t r = vsubq_f64 (z, shift);
140  float64x2_t d = vsubq_f64 (a, r);
141  float64x2_t d2 = vmulq_f64 (d, d);
142  float64x2_t r2 = vmulq_f64 (r, r);
143
144  float64x2_t p1 = r;
145  float64x2_t p2 = vfmsq_f64 (dat->p20, r2, vaddq_f64 (dat->p20, dat->p20));
146  float64x2_t p3 = vmulq_f64 (r, vfmaq_f64 (v_f64 (-0.5), r2, dat->p20));
147  float64x2_t p4 = vfmaq_f64 (dat->p41, r2, dat->p42);
148  p4 = vfmsq_f64 (dat->p40, r2, p4);
149  float64x2_t p5 = vfmaq_f64 (dat->p51, r2, dat->p52);
150  p5 = vmulq_f64 (r, vfmaq_f64 (vmulq_f64 (v_f64 (0.5), dat->p20), r2, p5));
151  /* Compute p_i using recurrence relation:
152     p_{i+2} = (p_i + r * Q_{i+1} * p_{i+1}) * R_{i+1}.  */
153  float64x2_t p6 = vfmaq_f64 (p4, p5, vmulq_laneq_f64 (r, dat->qr5, 0));
154  p6 = vmulq_laneq_f64 (p6, dat->qr5, 1);
155  float64x2_t p7 = vfmaq_f64 (p5, p6, vmulq_laneq_f64 (r, dat->qr6, 0));
156  p7 = vmulq_laneq_f64 (p7, dat->qr6, 1);
157  float64x2_t p8 = vfmaq_f64 (p6, p7, vmulq_laneq_f64 (r, dat->qr7, 0));
158  p8 = vmulq_laneq_f64 (p8, dat->qr7, 1);
159  float64x2_t p9 = vfmaq_f64 (p7, p8, vmulq_laneq_f64 (r, dat->qr8, 0));
160  p9 = vmulq_laneq_f64 (p9, dat->qr8, 1);
161  float64x2_t p10 = vfmaq_f64 (p8, p9, vmulq_laneq_f64 (r, dat->qr9, 0));
162  p10 = vmulq_laneq_f64 (p10, dat->qr9, 1);
163  /* Compute polynomial in d using pairwise Horner scheme.  */
164  float64x2_t p90 = vfmaq_f64 (p9, d, p10);
165  float64x2_t p78 = vfmaq_f64 (p7, d, p8);
166  float64x2_t p56 = vfmaq_f64 (p5, d, p6);
167  float64x2_t p34 = vfmaq_f64 (p3, d, p4);
168  float64x2_t p12 = vfmaq_f64 (p1, d, p2);
169  float64x2_t y = vfmaq_f64 (p78, d2, p90);
170  y = vfmaq_f64 (p56, d2, y);
171  y = vfmaq_f64 (p34, d2, y);
172  y = vfmaq_f64 (p12, d2, y);
173
174  y = vfmsq_f64 (e.erfc, e.scale, vfmsq_f64 (d, d2, y));
175
176  /* Offset equals 2.0 if sign, else 0.0.  */
177  uint64x2_t sign = vshrq_n_u64 (vreinterpretq_u64_f64 (x), 63);
178  float64x2_t off = vreinterpretq_f64_u64 (vshlq_n_u64 (sign, 62));
179  /* Copy sign and scale back in a single fma. Since the bit patterns do not
180     overlap, then logical or and addition are equivalent here.  */
181  float64x2_t fac = vreinterpretq_f64_u64 (
182      vsraq_n_u64 (vshlq_n_u64 (sign, 63), dat->table_scale, 1));
183
184#if WANT_SIMD_EXCEPT
185  if (unlikely (v_any_u64 (cmp)))
186    return special_case (xm, vfmaq_f64 (off, fac, y), cmp);
187#endif
188
189  return vfmaq_f64 (off, fac, y);
190}
191
192PL_SIG (V, D, 1, erfc, -6.0, 28.0)
193PL_TEST_ULP (V_NAME_D1 (erfc), 1.21)
194PL_TEST_SYM_INTERVAL (V_NAME_D1 (erfc), 0, 0x1p-26, 40000)
195PL_TEST_INTERVAL (V_NAME_D1 (erfc), 0x1p-26, 28.0, 40000)
196PL_TEST_INTERVAL (V_NAME_D1 (erfc), -0x1p-26, -6.0, 40000)
197PL_TEST_INTERVAL (V_NAME_D1 (erfc), 28.0, inf, 40000)
198PL_TEST_INTERVAL (V_NAME_D1 (erfc), -6.0, -inf, 40000)
199