aarch64/fpu: Add vector variants of erfc

Reviewed-by: Szabolcs Nagy <szabolcs.nagy@arm.com>

15 files changed, 4884 insertions, 1 deletions

diff --git a/sysdeps/aarch64/fpu/Makefile b/sysdeps/aarch64/fpu/Makefile
index e5f418ae42..e8af35099d 100644
--- a/sysdeps/aarch64/fpu/Makefile
+++ b/sysdeps/aarch64/fpu/Makefile
@@ -8,6 +8,7 @@ libmvec-supported-funcs = acos \
                           cos \
                           cosh \
                           erf \
+                          erfc \
                           exp \
                           exp10 \
                           exp2 \
@@ -39,7 +40,9 @@ libmvec-support = $(addsuffix f_advsimd,$(float-advsimd-funcs)) \
                   erff_data \
                   sv_erf_data \
                   sv_erff_data \
-                  v_exp_tail_data
+                  v_exp_tail_data \
+                  erfc_data \
+                  erfcf_data
 endif
 
 sve-cflags = -march=armv8-a+sve
diff --git a/sysdeps/aarch64/fpu/Versions b/sysdeps/aarch64/fpu/Versions
index 4dbf3d3244..3cb1b82bd2 100644
--- a/sysdeps/aarch64/fpu/Versions
+++ b/sysdeps/aarch64/fpu/Versions
@@ -104,6 +104,11 @@ libmvec {
     _ZGVnN4v_erff;
     _ZGVsMxv_erf;
     _ZGVsMxv_erff;
+    _ZGVnN2v_erfc;
+    _ZGVnN2v_erfcf;
+    _ZGVnN4v_erfcf;
+    _ZGVsMxv_erfc;
+    _ZGVsMxv_erfcf;
     _ZGVnN2v_sinh;
     _ZGVnN2v_sinhf;
     _ZGVnN4v_sinhf;
diff --git a/sysdeps/aarch64/fpu/advsimd_f32_protos.h b/sysdeps/aarch64/fpu/advsimd_f32_protos.h
index 4ff191c324..383c436972 100644
--- a/sysdeps/aarch64/fpu/advsimd_f32_protos.h
+++ b/sysdeps/aarch64/fpu/advsimd_f32_protos.h
@@ -26,6 +26,7 @@ libmvec_hidden_proto (V_NAME_F1(atanh));
 libmvec_hidden_proto (V_NAME_F1(cos));
 libmvec_hidden_proto (V_NAME_F1(cosh));
 libmvec_hidden_proto (V_NAME_F1(erf));
+libmvec_hidden_proto (V_NAME_F1(erfc));
 libmvec_hidden_proto (V_NAME_F1(exp10));
 libmvec_hidden_proto (V_NAME_F1(exp2));
 libmvec_hidden_proto (V_NAME_F1(exp));
diff --git a/sysdeps/aarch64/fpu/bits/math-vector.h b/sysdeps/aarch64/fpu/bits/math-vector.h
index 585e022082..e29b2d1c09 100644
--- a/sysdeps/aarch64/fpu/bits/math-vector.h
+++ b/sysdeps/aarch64/fpu/bits/math-vector.h
@@ -69,6 +69,10 @@
 # define __DECL_SIMD_erf __DECL_SIMD_aarch64
 # undef __DECL_SIMD_erff
 # define __DECL_SIMD_erff __DECL_SIMD_aarch64
+# undef __DECL_SIMD_erfc
+# define __DECL_SIMD_erfc __DECL_SIMD_aarch64
+# undef __DECL_SIMD_erfcf
+# define __DECL_SIMD_erfcf __DECL_SIMD_aarch64
 # undef __DECL_SIMD_exp
 # define __DECL_SIMD_exp __DECL_SIMD_aarch64
 # undef __DECL_SIMD_expf
@@ -153,6 +157,7 @@ __vpcs __f32x4_t _ZGVnN4v_atanhf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_cosf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_coshf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_erff (__f32x4_t);
+__vpcs __f32x4_t _ZGVnN4v_erfcf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_expf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_exp10f (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_exp2f (__f32x4_t);
@@ -176,6 +181,7 @@ __vpcs __f64x2_t _ZGVnN2v_atanh (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_cos (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_cosh (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_erf (__f64x2_t);
+__vpcs __f64x2_t _ZGVnN2v_erfc (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_exp (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_exp10 (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_exp2 (__f64x2_t);
@@ -204,6 +210,7 @@ __sv_f32_t _ZGVsMxv_atanhf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_cosf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_coshf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_erff (__sv_f32_t, __sv_bool_t);
+__sv_f32_t _ZGVsMxv_erfcf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_expf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_exp10f (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_exp2f (__sv_f32_t, __sv_bool_t);
@@ -227,6 +234,7 @@ __sv_f64_t _ZGVsMxv_atanh (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_cos (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_cosh (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_erf (__sv_f64_t, __sv_bool_t);
+__sv_f64_t _ZGVsMxv_erfc (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_exp (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_exp10 (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_exp2 (__sv_f64_t, __sv_bool_t);
diff --git a/sysdeps/aarch64/fpu/erfc_advsimd.c b/sysdeps/aarch64/fpu/erfc_advsimd.c
new file mode 100644
index 0000000000..548f21a3d6
--- /dev/null
+++ b/sysdeps/aarch64/fpu/erfc_advsimd.c
@@ -0,0 +1,201 @@
+/* Double-precision vector (Advanced SIMD) erfc function
+
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include "v_math.h"
+#include "vecmath_config.h"
+
+static const struct data
+{
+  uint64x2_t offset, table_scale;
+  float64x2_t max, shift;
+  float64x2_t p20, p40, p41, p42;
+  float64x2_t p51, p52;
+  float64x2_t qr5, qr6, qr7, qr8, qr9;
+#if WANT_SIMD_EXCEPT
+  float64x2_t uflow_bound;
+#endif
+} data = {
+  /* Set an offset so the range of the index used for lookup is 3487, and it
+     can be clamped using a saturated add on an offset index.
+     Index offset is 0xffffffffffffffff - asuint64(shift) - 3487.  */
+  .offset = V2 (0xbd3ffffffffff260),
+  .table_scale = V2 (0x37f0000000000000 << 1), /* asuint64 (2^-128) << 1.  */
+  .max = V2 (0x1.b3ep+4),		       /* 3487/128.  */
+  .shift = V2 (0x1p45),
+  .p20 = V2 (0x1.5555555555555p-2),  /* 1/3, used to compute 2/3 and 1/6.  */
+  .p40 = V2 (-0x1.999999999999ap-4), /* 1/10.  */
+  .p41 = V2 (-0x1.999999999999ap-2), /* 2/5.  */
+  .p42 = V2 (0x1.1111111111111p-3),  /* 2/15.  */
+  .p51 = V2 (-0x1.c71c71c71c71cp-3), /* 2/9.  */
+  .p52 = V2 (0x1.6c16c16c16c17p-5),  /* 2/45.  */
+  /* Qi = (i+1) / i, Ri = -2 * i / ((i+1)*(i+2)), for i = 5, ..., 9.  */
+  .qr5 = { 0x1.3333333333333p0, -0x1.e79e79e79e79ep-3 },
+  .qr6 = { 0x1.2aaaaaaaaaaabp0, -0x1.b6db6db6db6dbp-3 },
+  .qr7 = { 0x1.2492492492492p0, -0x1.8e38e38e38e39p-3 },
+  .qr8 = { 0x1.2p0, -0x1.6c16c16c16c17p-3 },
+  .qr9 = { 0x1.1c71c71c71c72p0, -0x1.4f2094f2094f2p-3 },
+#if WANT_SIMD_EXCEPT
+  .uflow_bound = V2 (0x1.a8b12fc6e4892p+4),
+#endif
+};
+
+#define TinyBound 0x4000000000000000 /* 0x1p-511 << 1.  */
+#define Off 0xfffffffffffff260	     /* 0xffffffffffffffff - 3487.  */
+
+struct entry
+{
+  float64x2_t erfc;
+  float64x2_t scale;
+};
+
+static inline struct entry
+lookup (uint64x2_t i)
+{
+  struct entry e;
+  float64x2_t e1 = vld1q_f64 ((float64_t *) (__erfc_data.tab - Off + i[0])),
+	      e2 = vld1q_f64 ((float64_t *) (__erfc_data.tab - Off + i[1]));
+  e.erfc = vuzp1q_f64 (e1, e2);
+  e.scale = vuzp2q_f64 (e1, e2);
+  return e;
+}
+
+#if WANT_SIMD_EXCEPT
+static float64x2_t VPCS_ATTR NOINLINE
+special_case (float64x2_t x, float64x2_t y, uint64x2_t cmp)
+{
+  return v_call_f64 (erfc, x, y, cmp);
+}
+#endif
+
+/* Optimized double-precision vector erfc(x).
+   Approximation based on series expansion near x rounded to
+   nearest multiple of 1/128.
+
+   Let d = x - r, and scale = 2 / sqrt(pi) * exp(-r^2). For x near r,
+
+   erfc(x) ~ erfc(r) - scale * d * poly(r, d), with
+
+   poly(r, d) = 1 - r d + (2/3 r^2 - 1/3) d^2 - r (1/3 r^2 - 1/2) d^3
+		+ (2/15 r^4 - 2/5 r^2 + 1/10) d^4
+		- r * (2/45 r^4 - 2/9 r^2 + 1/6) d^5
+		+ p6(r) d^6 + ... + p10(r) d^10
+
+   Polynomials p6(r) to p10(r) are computed using recurrence relation
+
+   2(i+1)p_i + 2r(i+2)p_{i+1} + (i+2)(i+3)p_{i+2} = 0,
+   with p0 = 1, and p1(r) = -r.
+
+   Values of erfc(r) and scale are read from lookup tables. Stored values
+   are scaled to avoid hitting the subnormal range.
+
+   Note that for x < 0, erfc(x) = 2.0 - erfc(-x).
+
+   Maximum measured error: 1.71 ULP
+   V_NAME_D1 (erfc)(0x1.46cfe976733p+4) got 0x1.e15fcbea3e7afp-608
+				       want 0x1.e15fcbea3e7adp-608.  */
+VPCS_ATTR
+float64x2_t V_NAME_D1 (erfc) (float64x2_t x)
+{
+  const struct data *dat = ptr_barrier (&data);
+
+#if WANT_SIMD_EXCEPT
+  /* |x| < 2^-511. Avoid fabs by left-shifting by 1.  */
+  uint64x2_t ix = vreinterpretq_u64_f64 (x);
+  uint64x2_t cmp = vcltq_u64 (vaddq_u64 (ix, ix), v_u64 (TinyBound));
+  /* x >= ~26.54 (into subnormal case and uflow case). Comparison is done in
+     integer domain to avoid raising exceptions in presence of nans.  */
+  uint64x2_t uflow = vcgeq_s64 (vreinterpretq_s64_f64 (x),
+				vreinterpretq_s64_f64 (dat->uflow_bound));
+  cmp = vorrq_u64 (cmp, uflow);
+  float64x2_t xm = x;
+  /* If any lanes are special, mask them with 0 and retain a copy of x to allow
+     special case handler to fix special lanes later. This is only necessary if
+     fenv exceptions are to be triggered correctly.  */
+  if (__glibc_unlikely (v_any_u64 (cmp)))
+    x = v_zerofy_f64 (x, cmp);
+#endif
+
+  float64x2_t a = vabsq_f64 (x);
+  a = vminq_f64 (a, dat->max);
+
+  /* Lookup erfc(r) and scale(r) in tables, e.g. set erfc(r) to 0 and scale to
+     2/sqrt(pi), when x reduced to r = 0.  */
+  float64x2_t shift = dat->shift;
+  float64x2_t z = vaddq_f64 (a, shift);
+
+  /* Clamp index to a range of 3487. A naive approach would use a subtract and
+     min. Instead we offset the table address and the index, then use a
+     saturating add.  */
+  uint64x2_t i = vqaddq_u64 (vreinterpretq_u64_f64 (z), dat->offset);
+
+  struct entry e = lookup (i);
+
+  /* erfc(x) ~ erfc(r) - scale * d * poly(r, d).  */
+  float64x2_t r = vsubq_f64 (z, shift);
+  float64x2_t d = vsubq_f64 (a, r);
+  float64x2_t d2 = vmulq_f64 (d, d);
+  float64x2_t r2 = vmulq_f64 (r, r);
+
+  float64x2_t p1 = r;
+  float64x2_t p2 = vfmsq_f64 (dat->p20, r2, vaddq_f64 (dat->p20, dat->p20));
+  float64x2_t p3 = vmulq_f64 (r, vfmaq_f64 (v_f64 (-0.5), r2, dat->p20));
+  float64x2_t p4 = vfmaq_f64 (dat->p41, r2, dat->p42);
+  p4 = vfmsq_f64 (dat->p40, r2, p4);
+  float64x2_t p5 = vfmaq_f64 (dat->p51, r2, dat->p52);
+  p5 = vmulq_f64 (r, vfmaq_f64 (vmulq_f64 (v_f64 (0.5), dat->p20), r2, p5));
+  /* Compute p_i using recurrence relation:
+     p_{i+2} = (p_i + r * Q_{i+1} * p_{i+1}) * R_{i+1}.  */
+  float64x2_t p6 = vfmaq_f64 (p4, p5, vmulq_laneq_f64 (r, dat->qr5, 0));
+  p6 = vmulq_laneq_f64 (p6, dat->qr5, 1);
+  float64x2_t p7 = vfmaq_f64 (p5, p6, vmulq_laneq_f64 (r, dat->qr6, 0));
+  p7 = vmulq_laneq_f64 (p7, dat->qr6, 1);
+  float64x2_t p8 = vfmaq_f64 (p6, p7, vmulq_laneq_f64 (r, dat->qr7, 0));
+  p8 = vmulq_laneq_f64 (p8, dat->qr7, 1);
+  float64x2_t p9 = vfmaq_f64 (p7, p8, vmulq_laneq_f64 (r, dat->qr8, 0));
+  p9 = vmulq_laneq_f64 (p9, dat->qr8, 1);
+  float64x2_t p10 = vfmaq_f64 (p8, p9, vmulq_laneq_f64 (r, dat->qr9, 0));
+  p10 = vmulq_laneq_f64 (p10, dat->qr9, 1);
+  /* Compute polynomial in d using pairwise Horner scheme.  */
+  float64x2_t p90 = vfmaq_f64 (p9, d, p10);
+  float64x2_t p78 = vfmaq_f64 (p7, d, p8);
+  float64x2_t p56 = vfmaq_f64 (p5, d, p6);
+  float64x2_t p34 = vfmaq_f64 (p3, d, p4);
+  float64x2_t p12 = vfmaq_f64 (p1, d, p2);
+  float64x2_t y = vfmaq_f64 (p78, d2, p90);
+  y = vfmaq_f64 (p56, d2, y);
+  y = vfmaq_f64 (p34, d2, y);
+  y = vfmaq_f64 (p12, d2, y);
+
+  y = vfmsq_f64 (e.erfc, e.scale, vfmsq_f64 (d, d2, y));
+
+  /* Offset equals 2.0 if sign, else 0.0.  */
+  uint64x2_t sign = vshrq_n_u64 (vreinterpretq_u64_f64 (x), 63);
+  float64x2_t off = vreinterpretq_f64_u64 (vshlq_n_u64 (sign, 62));
+  /* Copy sign and scale back in a single fma. Since the bit patterns do not
+     overlap, then logical or and addition are equivalent here.  */
+  float64x2_t fac = vreinterpretq_f64_u64 (
+      vsraq_n_u64 (vshlq_n_u64 (sign, 63), dat->table_scale, 1));
+
+#if WANT_SIMD_EXCEPT
+  if (__glibc_unlikely (v_any_u64 (cmp)))
+    return special_case (xm, vfmaq_f64 (off, fac, y), cmp);
+#endif
+
+  return vfmaq_f64 (off, fac, y);
+}
diff --git a/sysdeps/aarch64/fpu/erfc_data.c b/sysdeps/aarch64/fpu/erfc_data.c
new file mode 100644
index 0000000000..76a94e4681
--- /dev/null
+++ b/sysdeps/aarch64/fpu/erfc_data.c
@@ -0,0 +1,3519 @@
+/* Table for Advanced SIMD erfc
+
+   Copyright (C) 2024 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
+
+#include "vecmath_config.h"
+
+/* Lookup table used in erfc.
+   For each possible rounded input r (multiples of 1/128), between
+   r = 0.0 and r = ~27.0 (3488 values):
+   - the first entry __erfc_data.tab.erfc contains the values of erfc(r),
+   - the second entry __erfc_data.tab.scale contains the values of
+   2/sqrt(pi)*exp(-r^2). Both values may go into subnormal range, therefore
+   they are scaled by a large enough value 2^128 (fits in 8bit).  */
+const struct erfc_data __erfc_data = {
+  .tab = { { 0x1p128, 0x1.20dd750429b6dp128 },
+	   { 0x1.fb7c9030853b3p127, 0x1.20d8f1975c85dp128 },
+	   { 0x1.f6f9447be0743p127, 0x1.20cb67bd452c7p128 },
+	   { 0x1.f27640f9853d9p127, 0x1.20b4d8bac36c1p128 },
+	   { 0x1.edf3a9ba22dadp127, 0x1.209546ad13ccfp128 },
+	   { 0x1.e971a2c4436aep127, 0x1.206cb4897b148p128 },
+	   { 0x1.e4f05010eca8cp127, 0x1.203b261cd0053p128 },
+	   { 0x1.e06fd58842c7ep127, 0x1.2000a00ae3804p128 },
+	   { 0x1.dbf056fe2df35p127, 0x1.1fbd27cdc72d3p128 },
+	   { 0x1.d771f82f02f4ep127, 0x1.1f70c3b4f2cc8p128 },
+	   { 0x1.d2f4dcbc2f894p127, 0x1.1f1b7ae44867fp128 },
+	   { 0x1.ce792828eae5cp127, 0x1.1ebd5552f795bp128 },
+	   { 0x1.c9fefdd6eaf19p127, 0x1.1e565bca400d4p128 },
+	   { 0x1.c58681031eb6ap127, 0x1.1de697e413d29p128 },
+	   { 0x1.c10fd4c26e896p127, 0x1.1d6e14099944ap128 },
+	   { 0x1.bc9b1bfe82687p127, 0x1.1cecdb718d61cp128 },
+	   { 0x1.b82879728f11ep127, 0x1.1c62fa1e869b6p128 },
+	   { 0x1.b3b80fa82a4bbp127, 0x1.1bd07cdd189acp128 },
+	   { 0x1.af4a00f426daap127, 0x1.1b357141d95d5p128 },
+	   { 0x1.aade6f7378a0ep127, 0x1.1a91e5a748165p128 },
+	   { 0x1.a6757d08215d8p127, 0x1.19e5e92b964abp128 },
+	   { 0x1.a20f4b5626818p127, 0x1.19318bae53a04p128 },
+	   { 0x1.9dabfbc090901p127, 0x1.1874ddcdfce24p128 },
+	   { 0x1.994baf66747adp127, 0x1.17aff0e56ec1p128 },
+	   { 0x1.94ee8720076b6p127, 0x1.16e2d7093cd8cp128 },
+	   { 0x1.9094a37bbd66ep127, 0x1.160da304ed92fp128 },
+	   { 0x1.8c3e24bb73372p127, 0x1.153068581b781p128 },
+	   { 0x1.87eb2ad1a4032p127, 0x1.144b3b337c90cp128 },
+	   { 0x1.839bd55eaafc8p127, 0x1.135e3075d076bp128 },
+	   { 0x1.7f5043ae11862p127, 0x1.12695da8b5bdep128 },
+	   { 0x1.7b0894b3ea35cp127, 0x1.116cd8fd67618p128 },
+	   { 0x1.76c4e70a390e7p127, 0x1.1068b94962e5ep128 },
+	   { 0x1.728558ee694fcp127, 0x1.0f5d1602f7e41p128 },
+	   { 0x1.6e4a083ed132fp127, 0x1.0e4a073dc1b91p128 },
+	   { 0x1.6a13127843ec1p127, 0x1.0d2fa5a70c168p128 },
+	   { 0x1.65e094b3b2413p127, 0x1.0c0e0a8223359p128 },
+	   { 0x1.61b2aba3da093p127, 0x1.0ae54fa490723p128 },
+	   { 0x1.5d89739304dcfp127, 0x1.09b58f724416bp128 },
+	   { 0x1.59650860d6469p127, 0x1.087ee4d9ad247p128 },
+	   { 0x1.5545858029b39p127, 0x1.07416b4fbfe7cp128 },
+	   { 0x1.512b05f5006e1p127, 0x1.05fd3ecbec298p128 },
+	   { 0x1.4d15a4527fdc7p127, 0x1.04b27bc403d3p128 },
+	   { 0x1.49057ab900447p127, 0x1.03613f2812dafp128 },
+	   { 0x1.44faa2d42c4ap127, 0x1.0209a65e29545p128 },
+	   { 0x1.40f535d93160ep127, 0x1.00abcf3e187a9p128 },
+	   { 0x1.3cf54c850162p127, 0x1.fe8fb01a47307p127 },
+	   { 0x1.38faff1aa574ap127, 0x1.fbbbbef34b4b2p127 },
+	   { 0x1.35066561a275dp127, 0x1.f8dc092d58ff8p127 },
+	   { 0x1.311796a46f064p127, 0x1.f5f0cdaf15313p127 },
+	   { 0x1.2d2ea9aefb636p127, 0x1.f2fa4c16c0019p127 },
+	   { 0x1.294bb4cd4b2bdp127, 0x1.eff8c4b1375dbp127 },
+	   { 0x1.256ecdca212ccp127, 0x1.ecec7870ebca8p127 },
+	   { 0x1.219809edbd524p127, 0x1.e9d5a8e4c934ep127 },
+	   { 0x1.1dc77dfcacd02p127, 0x1.e6b4982f158b9p127 },
+	   { 0x1.19fd3e36ac96ap127, 0x1.e38988fc46e72p127 },
+	   { 0x1.16395e559e218p127, 0x1.e054be79d3042p127 },
+	   { 0x1.127bf18c8eadcp127, 0x1.dd167c4cf9d2ap127 },
+	   { 0x1.0ec50a86d0dd4p127, 0x1.d9cf06898cdafp127 },
+	   { 0x1.0b14bb6728cd8p127, 0x1.d67ea1a8b5368p127 },
+	   { 0x1.076b15c70aa28p127, 0x1.d325927fb9d89p127 },
+	   { 0x1.03c82ab5eb831p127, 0x1.cfc41e36c7df9p127 },
+	   { 0x1.002c0ab8a5018p127, 0x1.cc5a8a3fbea4p127 },
+	   { 0x1.f92d8b91d5cc7p126, 0x1.c8e91c4d01368p127 },
+	   { 0x1.f210d6a9a6a31p126, 0x1.c5701a484ef9dp127 },
+	   { 0x1.eb02147ce245cp126, 0x1.c1efca49a5011p127 },
+	   { 0x1.e40161b701275p126, 0x1.be68728e29d5ep127 },
+	   { 0x1.dd0ed9ea4bdd6p126, 0x1.bada596f25436p127 },
+	   { 0x1.d62a978f7c957p126, 0x1.b745c55905bf8p127 },
+	   { 0x1.cf54b4058455fp126, 0x1.b3aafcc27502ep127 },
+	   { 0x1.c88d479173ccep126, 0x1.b00a46237d5bep127 },
+	   { 0x1.c1d4695e87644p126, 0x1.ac63e7ecc1411p127 },
+	   { 0x1.bb2a2f7e5652p126, 0x1.a8b8287ec6a09p127 },
+	   { 0x1.b48eaee924501p126, 0x1.a5074e215762p127 },
+	   { 0x1.ae01fb7e55a66p126, 0x1.a1519efaf889ep127 },
+	   { 0x1.a78428050527ep126, 0x1.9d97610879642p127 },
+	   { 0x1.a115462cbbc17p126, 0x1.99d8da149c13fp127 },
+	   { 0x1.9ab5668e4930ap126, 0x1.96164fafd8de3p127 },
+	   { 0x1.946498acbd766p126, 0x1.925007283d7aap127 },
+	   { 0x1.8e22eaf68291ep126, 0x1.8e86458169af8p127 },
+	   { 0x1.87f06ac6960c4p126, 0x1.8ab94f6caa71dp127 },
+	   { 0x1.81cd2465e1d96p126, 0x1.86e9694134b9ep127 },
+	   { 0x1.7bb9230cb40b4p126, 0x1.8316d6f48133dp127 },
+	   { 0x1.75b470e454d35p126, 0x1.7f41dc12c9e89p127 },
+	   { 0x1.6fbf1708ba47cp126, 0x1.7b6abbb7aaf19p127 },
+	   { 0x1.69d91d8a595dap126, 0x1.7791b886e7403p127 },
+	   { 0x1.64028b7013867p126, 0x1.73b714a552763p127 },
+	   { 0x1.5e3b66b9405a9p126, 0x1.6fdb11b1e0c34p127 },
+	   { 0x1.5883b45fd2b63p126, 0x1.6bfdf0beddaf5p127 },
+	   { 0x1.52db785a98acap126, 0x1.681ff24b4ab04p127 },
+	   { 0x1.4d42b59f95afap126, 0x1.6441563c665d4p127 },
+	   { 0x1.47b96e267647ap126, 0x1.60625bd75d07bp127 },
+	   { 0x1.423fa2eb1cb59p126, 0x1.5c8341bb23767p127 },
+	   { 0x1.3cd553f045d45p126, 0x1.58a445da7c74cp127 },
+	   { 0x1.377a8042458d1p126, 0x1.54c5a57629dbp127 },
+	   { 0x1.322f25f9da2fdp126, 0x1.50e79d1749ac9p127 },
+	   { 0x1.2cf3423f15fdfp126, 0x1.4d0a6889dfd9fp127 },
+	   { 0x1.27c6d14c5e341p126, 0x1.492e42d78d2c5p127 },
+	   { 0x1.22a9ce717edcbp126, 0x1.4553664273d24p127 },
+	   { 0x1.1d9c3416d2b4bp126, 0x1.417a0c4049fdp127 },
+	   { 0x1.189dfbc07e69p126, 0x1.3da26d759aef5p127 },
+	   { 0x1.13af1e11be721p126, 0x1.39ccc1b136d5ap127 },
+	   { 0x1.0ecf92d046d22p126, 0x1.35f93fe7d1b3dp127 },
+	   { 0x1.09ff50e7b3f93p126, 0x1.32281e2fd1a92p127 },
+	   { 0x1.053e4e6d0c10bp126, 0x1.2e5991bd4cbfcp127 },
+	   { 0x1.008c80a24ff1p126, 0x1.2a8dcede3673bp127 },
+	   { 0x1.f7d3b7f436013p125, 0x1.26c508f6bd0ffp127 },
+	   { 0x1.eeaca836a27ccp125, 0x1.22ff727dd6f7bp127 },
+	   { 0x1.e5a3b7c9b56dap125, 0x1.1f3d3cf9ffe5ap127 },
+	   { 0x1.dcb8cae2d747fp125, 0x1.1b7e98fe26217p127 },
+	   { 0x1.d3ebc436b0f26p125, 0x1.17c3b626c7a12p127 },
+	   { 0x1.cb3c8500ea349p125, 0x1.140cc3173f007p127 },
+	   { 0x1.c2aaed0bfcfeep125, 0x1.1059ed7740313p127 },
+	   { 0x1.ba36dab91c0e9p125, 0x1.0cab61f084b93p127 },
+	   { 0x1.b1e02b082b72p125, 0x1.09014c2ca74dap127 },
+	   { 0x1.a9a6b99fc973bp125, 0x1.055bd6d32e8d7p127 },
+	   { 0x1.a18a60d56673ep125, 0x1.01bb2b87c6968p127 },
+	   { 0x1.998af9b56a3aep125, 0x1.fc3ee5d1524bp126 },
+	   { 0x1.91a85c0b65519p125, 0x1.f511a91a67d2ap126 },
+	   { 0x1.89e25e6a4cef9p125, 0x1.edeeee0959518p126 },
+	   { 0x1.8238d634c0127p125, 0x1.e6d6ffaa65a25p126 },
+	   { 0x1.7aab97a554544p125, 0x1.dfca26f5bbf88p126 },
+	   { 0x1.733a75d6e91b8p125, 0x1.d8c8aace11e63p126 },
+	   { 0x1.6be542ccffc2fp125, 0x1.d1d2cfff91594p126 },
+	   { 0x1.64abcf7c175b4p125, 0x1.cae8d93f1d7b7p126 },
+	   { 0x1.5d8debd20aacep125, 0x1.c40b0729ed548p126 },
+	   { 0x1.568b66be6f268p125, 0x1.bd3998457afdbp126 },
+	   { 0x1.4fa40e3af3674p125, 0x1.b674c8ffc6283p126 },
+	   { 0x1.48d7af53bc19fp125, 0x1.afbcd3afe8ab6p126 },
+	   { 0x1.4226162fbddd5p125, 0x1.a911f096fbc26p126 },
+	   { 0x1.3b8f0e1912f7p125, 0x1.a27455e14c93cp126 },
+	   { 0x1.351261854b991p125, 0x1.9be437a7de946p126 },
+	   { 0x1.2eafda1db784ap125, 0x1.9561c7f23a47bp126 },
+	   { 0x1.286740c7a7dabp125, 0x1.8eed36b886d93p126 },
+	   { 0x1.22385daca7f47p125, 0x1.8886b1e5ecfd1p126 },
+	   { 0x1.1c22f842ac1f2p125, 0x1.822e655b417e7p126 },
+	   { 0x1.1626d7543522p125, 0x1.7be47af1f5d89p126 },
+	   { 0x1.1043c1086777dp125, 0x1.75a91a7f4d2edp126 },
+	   { 0x1.0a797aeb152f2p125, 0x1.6f7c69