aarch64/fpu: Add vector variants of pow

Plus a small amount of moving includes around in order to be able to
remove duplicate definition of asuint64.

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

21 files changed, 2236 insertions, 12 deletions

diff --git a/sysdeps/aarch64/fpu/Makefile b/sysdeps/aarch64/fpu/Makefile
index 990d1135b9..234a6c457c 100644
--- a/sysdeps/aarch64/fpu/Makefile
+++ b/sysdeps/aarch64/fpu/Makefile
@@ -19,6 +19,7 @@ libmvec-supported-funcs = acos \
                           log10 \
                           log1p \
                           log2 \
+                          pow \
                           sin \
                           sinh \
                           tan \
@@ -44,7 +45,10 @@ libmvec-support = $(addsuffix f_advsimd,$(float-advsimd-funcs)) \
                   sv_erff_data \
                   v_exp_tail_data \
                   erfc_data \
-                  erfcf_data
+                  erfcf_data \
+                  v_pow_exp_data \
+                  v_pow_log_data \
+                  v_powf_data
 endif
 
 sve-cflags = -march=armv8-a+sve
diff --git a/sysdeps/aarch64/fpu/Versions b/sysdeps/aarch64/fpu/Versions
index 36a9e4df1e..cc15ce2d1e 100644
--- a/sysdeps/aarch64/fpu/Versions
+++ b/sysdeps/aarch64/fpu/Versions
@@ -119,6 +119,11 @@ libmvec {
     _ZGVnN2vv_hypot;
     _ZGVsMxvv_hypotf;
     _ZGVsMxvv_hypot;
+    _ZGVnN4vv_powf;
+    _ZGVnN2vv_powf;
+    _ZGVnN2vv_pow;
+    _ZGVsMxvv_powf;
+    _ZGVsMxvv_pow;
     _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 54858efd8a..097d403ffe 100644
--- a/sysdeps/aarch64/fpu/advsimd_f32_protos.h
+++ b/sysdeps/aarch64/fpu/advsimd_f32_protos.h
@@ -37,6 +37,7 @@ libmvec_hidden_proto (V_NAME_F1(log10));
 libmvec_hidden_proto (V_NAME_F1(log1p));
 libmvec_hidden_proto (V_NAME_F1(log2));
 libmvec_hidden_proto (V_NAME_F1(log));
+libmvec_hidden_proto (V_NAME_F2(pow));
 libmvec_hidden_proto (V_NAME_F1(sin));
 libmvec_hidden_proto (V_NAME_F1(sinh));
 libmvec_hidden_proto (V_NAME_F1(tan));
diff --git a/sysdeps/aarch64/fpu/atan2_advsimd.c b/sysdeps/aarch64/fpu/atan2_advsimd.c
index 2fd6164134..b1e7a9b8fc 100644
--- a/sysdeps/aarch64/fpu/atan2_advsimd.c
+++ b/sysdeps/aarch64/fpu/atan2_advsimd.c
@@ -17,6 +17,7 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
+#include "math_config.h"
 #include "v_math.h"
 #include "poly_advsimd_f64.h"
 
diff --git a/sysdeps/aarch64/fpu/atan2_sve.c b/sysdeps/aarch64/fpu/atan2_sve.c
index 04fa71fa37..ed9f683436 100644
--- a/sysdeps/aarch64/fpu/atan2_sve.c
+++ b/sysdeps/aarch64/fpu/atan2_sve.c
@@ -17,6 +17,7 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
+#include "math_config.h"
 #include "sv_math.h"
 #include "poly_sve_f64.h"
 
diff --git a/sysdeps/aarch64/fpu/bits/math-vector.h b/sysdeps/aarch64/fpu/bits/math-vector.h
index b1c024fe13..7484150131 100644
--- a/sysdeps/aarch64/fpu/bits/math-vector.h
+++ b/sysdeps/aarch64/fpu/bits/math-vector.h
@@ -113,6 +113,10 @@
 # define __DECL_SIMD_log2 __DECL_SIMD_aarch64
 # undef __DECL_SIMD_log2f
 # define __DECL_SIMD_log2f __DECL_SIMD_aarch64
+# undef __DECL_SIMD_pow
+# define __DECL_SIMD_pow __DECL_SIMD_aarch64
+# undef __DECL_SIMD_powf
+# define __DECL_SIMD_powf __DECL_SIMD_aarch64
 # undef __DECL_SIMD_sin
 # define __DECL_SIMD_sin __DECL_SIMD_aarch64
 # undef __DECL_SIMD_sinf
@@ -176,6 +180,7 @@ __vpcs __f32x4_t _ZGVnN4v_logf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_log10f (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_log1pf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_log2f (__f32x4_t);
+__vpcs __f32x4_t _ZGVnN4vv_powf (__f32x4_t, __f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_sinf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_sinhf (__f32x4_t);
 __vpcs __f32x4_t _ZGVnN4v_tanf (__f32x4_t);
@@ -202,6 +207,7 @@ __vpcs __f64x2_t _ZGVnN2v_log (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_log10 (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_log1p (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_log2 (__f64x2_t);
+__vpcs __f64x2_t _ZGVnN2vv_pow (__f64x2_t, __f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_sin (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_sinh (__f64x2_t);
 __vpcs __f64x2_t _ZGVnN2v_tan (__f64x2_t);
@@ -233,6 +239,7 @@ __sv_f32_t _ZGVsMxv_logf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_log10f (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_log1pf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_log2f (__sv_f32_t, __sv_bool_t);
+__sv_f32_t _ZGVsMxvv_powf (__sv_f32_t, __sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_sinf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_sinhf (__sv_f32_t, __sv_bool_t);
 __sv_f32_t _ZGVsMxv_tanf (__sv_f32_t, __sv_bool_t);
@@ -259,6 +266,7 @@ __sv_f64_t _ZGVsMxv_log (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_log10 (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_log1p (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_log2 (__sv_f64_t, __sv_bool_t);
+__sv_f64_t _ZGVsMxvv_pow (__sv_f64_t, __sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_sin (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_sinh (__sv_f64_t, __sv_bool_t);
 __sv_f64_t _ZGVsMxv_tan (__sv_f64_t, __sv_bool_t);
diff --git a/sysdeps/aarch64/fpu/finite_pow.h b/sysdeps/aarch64/fpu/finite_pow.h
new file mode 100644
index 0000000000..84c93d4048
--- /dev/null
+++ b/sysdeps/aarch64/fpu/finite_pow.h
@@ -0,0 +1,373 @@
+/* Double-precision x^y 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 "math_config.h"
+
+/* Scalar version of pow used for fallbacks in vector implementations.  */
+
+/* Data is defined in v_pow_log_data.c.  */
+#define N_LOG (1 << V_POW_LOG_TABLE_BITS)
+#define Off 0x3fe6955500000000
+#define As __v_pow_log_data.poly
+
+/* Data is defined in v_pow_exp_data.c.  */
+#define N_EXP (1 << V_POW_EXP_TABLE_BITS)
+#define SignBias (0x800 << V_POW_EXP_TABLE_BITS)
+#define SmallExp 0x3c9 /* top12(0x1p-54).  */
+#define BigExp 0x408   /* top12(512.0).  */
+#define ThresExp 0x03f /* BigExp - SmallExp.  */
+#define InvLn2N __v_pow_exp_data.n_over_ln2
+#define Ln2HiN __v_pow_exp_data.ln2_over_n_hi
+#define Ln2LoN __v_pow_exp_data.ln2_over_n_lo
+#define SBits __v_pow_exp_data.sbits
+#define Cs __v_pow_exp_data.poly
+
+/* Constants associated with pow.  */
+#define SmallPowX 0x001 /* top12(0x1p-126).  */
+#define BigPowX 0x7ff	/* top12(INFINITY).  */
+#define ThresPowX 0x7fe /* BigPowX - SmallPowX.  */
+#define SmallPowY 0x3be /* top12(0x1.e7b6p-65).  */
+#define BigPowY 0x43e	/* top12(0x1.749p62).  */
+#define ThresPowY 0x080 /* BigPowY - SmallPowY.  */
+
+/* Top 12 bits of a double (sign and exponent bits).  */
+static inline uint32_t
+top12 (double x)
+{
+  return asuint64 (x) >> 52;
+}
+
+/* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
+   additional 15 bits precision.  IX is the bit representation of x, but
+   normalized in the subnormal range using the sign bit for the exponent.  */
+static inline double
+log_inline (uint64_t ix, double *tail)
+{
+  /* 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.  */
+  uint64_t tmp = ix - Off;
+  int i = (tmp >> (52 - V_POW_LOG_TABLE_BITS)) & (N_LOG - 1);
+  int k = (int64_t) tmp >> 52; /* arithmetic shift.  */
+  uint64_t iz = ix - (tmp & 0xfffULL << 52);
+  double z = asdouble (iz);
+  double kd = (double) k;
+
+  /* log(x) = k*Ln2 + log(c) + log1p(z/c-1).  */
+  double invc = __v_pow_log_data.invc[i];
+  double logc = __v_pow_log_data.logc[i];
+  double logctail = __v_pow_log_data.logctail[i];
+
+  /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
+     |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible.  */
+  double r = fma (z, invc, -1.0);
+
+  /* k*Ln2 + log(c) + r.  */
+  double t1 = kd * __v_pow_log_data.ln2_hi + logc;
+  double t2 = t1 + r;
+  double lo1 = kd * __v_pow_log_data.ln2_lo + logctail;
+  double lo2 = t1 - t2 + r;
+
+  /* Evaluation is optimized assuming superscalar pipelined execution.  */
+  double ar = As[0] * r;
+  double ar2 = r * ar;
+  double ar3 = r * ar2;
+  /* k*Ln2 + log(c) + r + A[0]*r*r.  */
+  double hi = t2 + ar2;
+  double lo3 = fma (ar, r, -ar2);
+  double lo4 = t2 - hi + ar2;
+  /* p = log1p(r) - r - A[0]*r*r.  */
+  double p = (ar3
+	      * (As[1] + r * As[2]
+		 + ar2 * (As[3] + r * As[4] + ar2 * (As[5] + r * As[6]))));
+  double lo = lo1 + lo2 + lo3 + lo4 + p;
+  double y = hi + lo;
+  *tail = hi - y + lo;
+  return y;
+}
+
+/* Handle cases that may overflow or underflow when computing the result that
+   is scale*(1+TMP) without intermediate rounding.  The bit representation of
+   scale is in SBITS, however it has a computed exponent that may have
+   overflown into the sign bit so that needs to be adjusted before using it as
+   a double.  (int32_t)KI is the k used in the argument reduction and exponent
+   adjustment of scale, positive k here means the result may overflow and
+   negative k means the result may underflow.  */
+static inline double
+special_case (double tmp, uint64_t sbits, uint64_t ki)
+{
+  double scale, y;
+
+  if ((ki & 0x80000000) == 0)
+    {
+      /* k > 0, the exponent of scale might have overflowed by <= 460.  */
+      sbits -= 1009ull << 52;
+      scale = asdouble (sbits);
+      y = 0x1p1009 * (scale + scale * tmp);
+      return y;
+    }
+  /* k < 0, need special care in the subnormal range.  */
+  sbits += 1022ull << 52;
+  /* Note: sbits is signed scale.  */
+  scale = asdouble (sbits);
+  y = scale + scale * tmp;
+#if WANT_SIMD_EXCEPT
+  if (fabs (y) < 1.0)
+    {
+      /* Round y to the right precision before scaling it into the subnormal
+	 range to avoid double rounding that can cause 0.5+E/2 ulp error where
+	 E is the worst-case ulp error outside the subnormal range.  So this
+	 is only useful if the goal is better than 1 ulp worst-case error.  */
+      double hi, lo, one = 1.0;
+      if (y < 0.0)
+	one = -1.0;
+      lo = scale - y + scale * tmp;
+      hi = one + y;
+      lo = one - hi + y + lo;
+      y = (hi + lo) - one;
+      /* Fix the sign of 0.  */
+      if (y == 0.0)
+	y = asdouble (sbits & 0x8000000000000000);
+      /* The underflow exception needs to be signaled explicitly.  */
+      force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022);
+    }
+#endif
+  y = 0x1p-1022 * y;
+  return y;
+}
+
+/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+   The sign_bias argument is SignBias or 0 and sets the sign to -1 or 1.  */
+static inline double
+exp_inline (double x, double xtail, uint32_t sign_bias)
+{
+  uint32_t abstop = top12 (x) & 0x7ff;
+  if (__glibc_unlikely (abstop - SmallExp >= ThresExp))
+    {
+      if (abstop - SmallExp >= 0x80000000)
+	{
+	  /* Avoid spurious underflow for tiny x.  */
+	  /* Note: 0 is common input.  */
+	  return sign_bias ? -1.0 : 1.0;
+	}
+      if (abstop >= top12 (1024.0))
+	{
+	  /* Note: inf and nan are already handled.  */
+	  /* Skip errno handling.  */
+#if WANT_SIMD_EXCEPT
+	  return asuint64 (x) >> 63 ? __math_uflow (sign_bias)
+				    : __math_oflow (sign_bias);
+#else
+	  double res_uoflow = asuint64 (x) >> 63 ? 0.0 : INFINITY;
+	  return sign_bias ? -res_uoflow : res_uoflow;
+#endif
+	}
+      /* Large x is special cased below.  */
+      abstop = 0;
+    }
+
+  /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
+  /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
+  double z = InvLn2N * x;
+  double kd = round (z);
+  uint64_t ki = lround (z);
+  double r = x - kd * Ln2HiN - kd * Ln2LoN;
+  /* The code assumes 2^-200 < |xtail| < 2^-8/N.  */
+  r += xtail;
+  /* 2^(k/N) ~= scale.  */
+  uint64_t idx = ki & (N_EXP - 1);
+  uint64_t top = (ki + sign_bias) << (52 - V_POW_EXP_TABLE_BITS);
+  /* This is only a valid scale when -1023*N < k < 1024*N.  */
+  uint64_t sbits = SBits[idx] + top;
+  /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1).  */
+  /* Evaluation is optimized assuming superscalar pipelined execution.  */
+  double r2 = r * r;
+  double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]);
+  if (__glibc_unlikely (abstop == 0))
+    return special_case (tmp, sbits, ki);
+  double scale = asdouble (sbits);
+  /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+     is no spurious underflow here even without fma.  */
+  return scale + scale * tmp;
+}
+
+/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+   A version of exp_inline that is not inlined and for which sign_bias is
+   equal to 0.  */
+static double NOINLINE
+exp_nosignbias (double x, double xtail)
+{
+  uint32_t abstop = top12 (x) & 0x7ff;
+  if (__glibc_unlikely (abstop - SmallExp >= ThresExp))
+    {
+      /* Avoid spurious underflow for tiny x.  */
+      if (abstop - SmallExp >= 0x80000000)
+	return 1.0;
+      /* Note: inf and nan are already handled.  */
+      if (abstop >= top12 (1024.0))
+#if WANT_SIMD_EXCEPT
+	return asuint64 (x) >> 63 ? __math_uflow (0) : __math_oflow (0);
+#else
+	return asuint64 (x) >> 63 ? 0.0 : INFINITY;
+#endif
+      /* Large x is special cased below.  */
+      abstop = 0;
+    }
+
+  /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
+  /* x = ln2/N*k + r, with k integer and r in [-ln2/2N, ln2/2N].  */
+  double z = InvLn2N * x;
+  double kd = round (z);
+  uint64_t ki = lround (z);
+  double r = x - kd * Ln2HiN - kd * Ln2LoN;
+  /* The code assumes 2^-200 < |xtail| < 2^-8/N.  */
+  r += xtail;
+  /* 2^(k/N) ~= scale.  */
+  uint64_t idx = ki & (N_EXP - 1);
+  uint64_t top = ki << (52 - V_POW_EXP_TABLE_BITS);
+  /* This is only a valid scale when -1023*N < k < 1024*N.  */
+  uint64_t sbits = SBits[idx] + top;
+  /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
+  double r2 = r * r;
+  double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]);
+  if (__glibc_unlikely (abstop == 0))
+    return special_case (tmp, sbits, ki);
+  double scale = asdouble (sbits);
+  /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+     is no spurious underflow here even without fma.  */
+  return scale + scale * tmp;
+}
+
+/* Returns 0 if not int, 1 if odd int, 2 if even int.  The argument is
+   the bit representation of a non-zero finite floating-point value.  */
+static inline int
+checkint (uint64_t iy)
+{
+  int e = iy >> 52 & 0x7ff;
+  if (e < 0x3ff)
+    return 0;
+  if (e > 0x3ff + 52)
+    return 2;
+  if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
+    return 0;
+  if (iy & (1ULL << (0x3ff + 52 - e)))
+    return 1;
+  return 2;
+}
+
+/* Returns 1 if input is the bit representation of 0, infinity or nan.  */
+static inline int
+zeroinfnan (uint64_t i)
+{
+  return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1;
+}
+
+static double NOINLINE
+pow_scalar_special_case (double x, double y)
+{
+  uint32_t sign_bias = 0;
+  uint64_t ix, iy;
+  uint32_t topx, topy;
+
+  ix = asuint64 (x);
+  iy = asuint64 (y);
+  topx = top12 (x);
+  topy = top12 (y);
+  if (__glibc_unlikely (topx - SmallPowX >= ThresPowX
+		|| (topy & 0x7ff) - SmallPowY >= ThresPowY))
+    {
+      /* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then pow(x,y) = inf/0
+	 and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1.  */
+      /* Special cases: (x < 0x1p-126 or inf or nan) or
+	 (|y| < 0x1p-65 or |y| >= 0x1p63 or nan).  */
+      if (__glibc_unlikely (zeroinfnan (iy)))
+	{
+	  if (2 * iy == 0)
+	    return issignaling_inline (x) ? x + y : 1.0;
+	  if (ix == asuint64 (1.0))
+	    return issignaling_inline (y) ? x + y : 1.0;
+	  if (2 * ix > 2 * asuint64 (INFINITY)
+	      || 2 * iy > 2 * asuint64 (INFINITY))
+	    return x + y;
+	  if (2 * ix == 2 * asuint64 (1.0))
+	    return 1.0;
+	  if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63))
+	    return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf.  */
+	  return y * y;
+	}
+      if (__glibc_unlikely (zeroinfnan (ix)))
+	{
+	  double x2 = x * x;
+	  if (ix >> 63 && checkint (iy) == 1)
+	    {
+	      x2 = -x2;
+	      sign_bias = 1;
+	    }
+#if WANT_SIMD_EXCEPT
+	  if (2 * ix == 0 && iy >> 63)
+	    return __math_divzero (sign_bias);
+#endif
+	  return iy >> 63 ? 1 / x2 : x2;
+	}
+      /* Here x and y are non-zero finite.  */
+      if (ix >> 63)
+	{
+	  /* Finite x < 0.  */
+	  int yint = checkint (iy);
+	  if (yint == 0)
+#if WANT_SIMD_EXCEPT
+	    return __math_invalid (x);
+#else
+	    return __builtin_nan ("");
+#endif
+	  if (yint == 1)
+	    sign_bias = SignBias;
+	  ix &= 0x7fffffffffffffff;
+	  topx &= 0x7ff;
+	}
+      if ((topy & 0x7ff) - SmallPowY >= ThresPowY)
+	{
+	  /* Note: sign_bias == 0 here because y is not odd.  */
+	  if (ix == asuint64 (1.0))
+	    return 1.0;
+	  /* |y| < 2^-65, x^y ~= 1 + y*log(x).  */
+	  if ((topy & 0x7ff) < SmallPowY)
+	    return 1.0;
+#if WANT_SIMD_EXCEPT
+	  return (ix > asuint64 (1.0)) == (topy < 0x800) ? __math_oflow (0)
+							 : __math_uflow (0);
+#else
+	  return (ix > asuint64 (1.0)) == (topy < 0x800) ? INFINITY : 0;
+#endif
+	}
+      if (topx == 0)
+	{
+	  /* Normalize subnormal x so exponent becomes negative.  */
+	  ix = asuint64 (x * 0x1p52);
+	  ix &= 0x7fffffffffffffff;
+	  ix -= 52ULL << 52;
+	}
+    }
+
+  double lo;
+  double hi = log_inline (ix, &lo);
+  double ehi = y * hi;
+  double elo = y * l