aarch64: Optimise AdvSIMD log

Optimise AdvSIMD log by vectorising the special case.
For subnormal input values, use the same scaling technique as
described in the single precision equivalent.
Then check for inf, nan and x<=0.

diff --git a/sysdeps/aarch64/fpu/log_advsimd.c b/sysdeps/aarch64/fpu/log_advsimd.c
index 92c2c1d39482cb2aa822fc1b0d1a5da7c01eee10..bc293e3890629de3799e68b7fa77a36efa15d187 100644 (file)
--- a/sysdeps/aarch64/fpu/log_advsimd.c
+++ b/sysdeps/aarch64/fpu/log_advsimd.c
@@ -22,9 +22,11 @@
 static const struct data
 {
   uint64x2_t off, sign_exp_mask, offset_lower_bound;
-  uint32x4_t special_bound;
+  uint64x2_t special_bound;
+  uint32x4_t special_bound_u32;
   float64x2_t c0, c2;
   double c1, c3, ln2, c4;
+  float64x2_t pinf, minf, nan;
 } data = {
   /* Rel error: 0x1.6272e588p-56 in [ -0x1.fc1p-9 0x1.009p-8 ].  */
   .c0 = V2 (-0x1.ffffffffffff7p-2),
@@ -39,7 +41,12 @@ static const struct data
      optimised register use subnormals are detected after offset has been
      subtracted, so lower bound - offset (which wraps around).  */
   .offset_lower_bound = V2 (0x0010000000000000 - 0x3fe6900900000000),
-  .special_bound = V4 (0x7fe00000), /* asuint64(inf) -  asuint64(0x1p-126).  */
+  .special_bound = V2 (0x7ffe000000000000),
+  .special_bound_u32
+  = V4 (0x7fe00000), /* asuint64(inf) -  asuint64(0x1p-1022).  */
+  .pinf = V2 (INFINITY),
+  .minf = V2 (-INFINITY),
+  .nan = V2 (NAN),
 };
 
 #define N (1 << V_LOG_TABLE_BITS)
@@ -65,28 +72,9 @@ lookup (uint64x2_t i)
   return e;
 }
 
-static float64x2_t VPCS_ATTR NOINLINE
-special_case (float64x2_t hi, uint64x2_t u_off, float64x2_t y, float64x2_t r2,
-             uint32x2_t special, const struct data *d)
+static inline float64x2_t VPCS_ATTR
+inline_log (uint64x2_t u, uint64x2_t u_off, const struct data *d)
 {
-  float64x2_t x = vreinterpretq_f64_u64 (vaddq_u64 (u_off, d->off));
-  return v_call_f64 (log, x, vfmaq_f64 (hi, y, r2), vmovl_u32 (special));
-}
-
-/* Double-precision vector log routine.
-   The maximum observed error is 2.17 ULP:
-   _ZGVnN2v_log(0x1.a6129884398a3p+0) got 0x1.ffffff1cca043p-2
-                                    want 0x1.ffffff1cca045p-2.  */
-float64x2_t VPCS_ATTR V_NAME_D1 (log) (float64x2_t x)
-{
-  const struct data *d = ptr_barrier (&data);
-
-  /* To avoid having to mov x out of the way, keep u after offset has been
-     applied, and recover x by adding the offset back in the special-case
-     handler.  */
-  uint64x2_t u = vreinterpretq_u64_f64 (x);
-  uint64x2_t u_off = vsubq_u64 (u, d->off);
-
   /* x = 2^k z; where z is in range [Off,2*Off) and exact.
      The range is split into N subintervals.
      The ith subinterval contains z and c is near its center.  */
@@ -96,9 +84,6 @@ float64x2_t VPCS_ATTR V_NAME_D1 (log) (float64x2_t x)
 
   struct entry e = lookup (u_off);
 
-  uint32x2_t special = vcge_u32 (vsubhn_u64 (u_off, d->offset_lower_bound),
-                                vget_low_u32 (d->special_bound));
-
   /* log(x) = log1p(z/c-1) + log(c) + k*Ln2.  */
   float64x2_t r = vfmaq_f64 (v_f64 (-1.0), z, e.invc);
   float64x2_t kd = vcvtq_f64_s64 (k);
@@ -115,7 +100,58 @@ float64x2_t VPCS_ATTR V_NAME_D1 (log) (float64x2_t x)
   y = vfmaq_laneq_f64 (y, r2, ln2_and_c4, 1);
   y = vfmaq_f64 (p, r2, y);
 
-  if (__glibc_unlikely (v_any_u32h (special)))
-    return special_case (hi, u_off, y, r2, special, d);
   return vfmaq_f64 (hi, y, r2);
 }
+
+static inline float64x2_t VPCS_ATTR
+special_case (uint64x2_t u_off, const struct data *d)
+{
+  float64x2_t x = vreinterpretq_f64_u64 (vaddq_u64 (u_off, d->off));
+  /* If x is special, compute 2log(sqrt(x)), else compute log(x).
+     x might be subnormal, and sqrting it makes it larger.
+     And the above two expressions are equivalent.  */
+  uint64x2_t special
+      = vcgeq_u64 (vsubq_u64 (u_off, d->offset_lower_bound), d->special_bound);
+  float64x2_t x_sqrt = vbslq_f64 (special, vsqrtq_f64 (x), x);
+
+  u_off = vsubq_u64 (vreinterpretq_u64_f64 (x_sqrt), d->off);
+
+  /* Don't pass u into this, it isn't using x_sqrt.  */
+  float64x2_t y = inline_log (vreinterpretq_u64_f64 (x_sqrt), u_off, d);
+
+  y = vbslq_f64 (special, vmulq_f64 (y, v_f64 (2.0f)), y);
+
+  /* Is true for +/- inf, +/- nan as well as all negative numbers.  */
+  uint64x2_t is_infnan
+      = vcgeq_u64 (vreinterpretq_u64_f64 (x), vreinterpretq_u64_f64 (d->pinf));
+  uint64x2_t infnan_or_zero = vorrq_u64 (is_infnan, vceqzq_f64 (x));
+
+  y = vbslq_f64 (infnan_or_zero, d->nan, y);
+  uint64x2_t ret_pinf = vceqq_f64 (x, d->pinf);
+  uint64x2_t ret_minf = vceqzq_f64 (x);
+  y = vbslq_f64 (ret_pinf, d->pinf, y);
+  y = vbslq_f64 (ret_minf, d->minf, y);
+  return y;
+}
+
+/* Double-precision vector log routine.
+   The maximum observed error is 1.67 + 0.5 ULP:
+   _ZGVnN2v_log(0x1.a6129884398a3p+0) got 0x1.ffffff1cca043p-2
+                                    want 0x1.ffffff1cca045p-2.  */
+float64x2_t VPCS_ATTR V_NAME_D1 (log) (float64x2_t x)
+{
+  const struct data *d = ptr_barrier (&data);
+
+  /* To avoid having to mov x out of the way, keep u after offset has been
+     applied, and recover x by adding the offset back in the special-case
+     handler.  */
+  uint64x2_t u = vreinterpretq_u64_f64 (x);
+  uint64x2_t u_off = vsubq_u64 (u, d->off);
+
+  uint32x2_t special_u32 = vcge_u32 (vsubhn_u64 (u_off, d->offset_lower_bound),
+                                    vget_low_u32 (d->special_bound_u32));
+
+  if (__glibc_unlikely (v_any_u32h (special_u32)))
+    return special_case (u_off, d);
+  return inline_log (u, u_off, d);
+}