math: Consolidate CORE-MATH double-double routines

For lgamma and tgamma the muldd, mulddd, and polydd are renamed
to muldd2, mulddd2, and polydd2 respectively.

Checked on aarch64-linux-gnu and x86_64-linux-gnu.

Reviewed-by: DJ Delorie <dj@redhat.com>

diff --git a/sysdeps/ieee754/dbl-64/ddcoremath.h b/sysdeps/ieee754/dbl-64/ddcoremath.h
new file mode 100644 (file)
index 0000000..7c50b04
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/ddcoremath.h
@@ -0,0 +1,213 @@
+/* Double-double common routines used in correctly rounded implementations.
+
+Copyright (c) 2023-2025 Alexei Sibidanov.
+
+The original version of this file was copied from the CORE-MATH
+project (file src/binary64/acosh/acosh.c, revision 6d87ca23).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE. */
+
+#ifndef _DD_COREMATH_H
+#define _DD_COREMATH_H
+
+/* References:
+   [1] Tight and rigourous error bounds for basic building blocks of
+       double-word arithmetic, by Mioara Joldeş, Jean-Michel Muller,
+       and Valentina Popescu, ACM Transactions on Mathematical Software,
+       44(2), 2017.
+*/
+
+static inline double
+fasttwosum (double x, double y, double *e)
+{
+  double s = x + y, z = s - x;
+  *e = y - z;
+  return s;
+}
+
+static inline double
+fasttwosub (double x, double y, double *e)
+{
+  double s = x - y, z = x - s;
+  *e = z - y;
+  return s;
+}
+
+static inline double
+twosum (double x, double y, double *e)
+{
+  if (__glibc_likely (fabs (x) > fabs (y)))
+    return fasttwosum (x, y, e);
+  else
+    return fasttwosum (y, x, e);
+}
+
+static inline double
+fastsum (double xh, double xl, double yh, double yl, double *e)
+{
+  double sl, sh = fasttwosum (xh, yh, &sl);
+  *e = (xl + yl) + sl;
+  return sh;
+}
+
+static inline double
+sumdd (double xh, double xl, double yh, double yl, double *e)
+{
+  double sl, sh;
+  if (__glibc_likely (fabs (xh) > fabs (yh)))
+    sh = fasttwosum (xh, yh, &sl);
+  else
+    sh = fasttwosum (yh, xh, &sl);
+  *e = (xl + yl) + sl;
+  return sh;
+}
+
+static inline double
+adddd (double xh, double xl, double ch, double cl, double *l)
+{
+  double s = xh + ch, d = s - xh;
+  *l = ((ch - d) + (xh + (d - s))) + (xl + cl);
+  return s;
+}
+
+/* This function implements Algorithm 10 (DWTimesDW1) from [1]
+   Its relative error (for round-to-nearest ties-to-even) is bounded by 5u^2
+   (Theorem 2.6 of [2]), where u = 2^-53 for double precision,
+   assuming xh = RN(xh + xl), which implies |xl| <= 1/2 ulp(xh),
+   and similarly for ch, cl. */
+static inline double
+muldd (double xh, double xl, double ch, double cl, double *l)
+{
+  double ahlh = ch * xl, alhh = cl * xh, ahhh = ch * xh,
+        ahhl = fma (ch, xh, -ahhh);
+  ahhl += alhh + ahlh;
+  ch = ahhh + ahhl;
+  *l = (ahhh - ch) + ahhl;
+  return ch;
+}
+
+static inline double
+muldd2 (double xh, double xl, double ch, double cl, double *l)
+{
+  double ahhh = ch * xh;
+  *l = (ch * xl + cl * xh) + fma (ch, xh, -ahhh);
+  return ahhh;
+}
+
+static inline double
+muldd3 (double xh, double xl, double yh, double yl, double *l)
+{
+  double ch = xh * yh, cl1 = fma (xh, yh, -ch);
+  double tl0 = xl * yl;
+  double tl1 = tl0 + xh * yl;
+  double cl2 = tl1 + xl * yh;
+  double cl3 = cl1 + cl2;
+  return fasttwosum (ch, cl3, l);
+}
+
+static inline double
+mulddd (double xh, double xl, double ch, double *l)
+{
+  double ahlh = ch * xl, ahhh = ch * xh, ahhl = fma (ch, xh, -ahhh);
+  ahhl += ahlh;
+  ch = ahhh + ahhl;
+  *l = (ahhh - ch) + ahhl;
+  return ch;
+}
+
+static inline double
+mulddd2 (double x, double ch, double cl, double *l)
+{
+  double ahhh = ch * x;
+  *l = cl * x + fma (ch, x, -ahhh);
+  return ahhh;
+}
+
+static inline double
+polydd (double xh, double xl, int n, const double c[][2], double *l)
+{
+  int i = n - 1;
+  double ch = c[i][0] + *l, cl = ((c[i][0] - ch) + *l) + c[i][1];
+  while (--i >= 0)
+    {
+      ch = muldd (xh, xl, ch, cl, &cl);
+      double th = ch + c[i][0], tl = (c[i][0] - th) + ch;
+      ch = th;
+      cl += tl + c[i][1];
+    }
+  *l = cl;
+  return ch;
+}
+
+static inline double
+polydd2 (double xh, double xl, int n, const double c[][2], double *l)
+{
+  int i = n - 1;
+  double cl, ch = fasttwosum (c[i][0], *l, &cl);
+  cl += c[i][1];
+  while (--i >= 0)
+    {
+      ch = muldd2 (xh, xl, ch, cl, &cl);
+      ch = fastsum (c[i][0], c[i][1], ch, cl, &cl);
+    }
+  *l = cl;
+  return ch;
+}
+
+static inline double
+polyddd (double x, int n, const double c[][2], double *l)
+{
+  int i = n - 1;
+  double cl, ch = fasttwosum (c[i][0], *l, &cl);
+  cl += c[i][1];
+  while (--i >= 0)
+    {
+      ch = mulddd2 (x, ch, cl, &cl);
+      ch = sumdd (c[i][0], c[i][1], ch, cl, &cl);
+    }
+  *l = cl;
+  return ch;
+}
+
+static inline double
+polydddfst (double x, int n, const double c[][2], double *l)
+{
+  int i = n - 1;
+  double cl, ch = fasttwosum (c[i][0], *l, &cl);
+  cl += c[i][1];
+  while (--i >= 0)
+    {
+      ch = mulddd2 (x, ch, cl, &cl);
+      ch = fastsum (c[i][0], c[i][1], ch, cl, &cl);
+    }
+  *l = cl;
+  return ch;
+}
+
+static inline double
+polyd (double x, int n, const double c[][2])
+{
+  int i = n - 1;
+  double ch = c[i][0];
+  while (--i >= 0)
+    ch = c[i][0] + x * ch;
+  return ch;
+}
+
+#endif

diff --git a/sysdeps/ieee754/dbl-64/e_acosh.c b/sysdeps/ieee754/dbl-64/e_acosh.c
index faeef67e7e78dff94e92dc678b00aa6bdaa57773..959d347238d631e0690c2efd0e9b3c1f2af98231 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_acosh.c
+++ b/sysdeps/ieee754/dbl-64/e_acosh.c
@@ -41,64 +41,7 @@ SOFTWARE. */
 #include <libm-alias-finite.h>
 #include "math_config.h"
 #include "s_asincosh_data.h"
-
-static inline double
-fasttwosum (double x, double y, double *e)
-{
-  double s = x + y, z = s - x;
-  *e = y - z;
-  return s;
-}
-
-static inline double
-adddd (double xh, double xl, double ch, double cl, double *l)
-{
-  double s = xh + ch, d = s - xh;
-  *l = ((ch - d) + (xh + (d - s))) + (xl + cl);
-  return s;
-}
-
-/* This function implements Algorithm 10 (DWTimesDW1) from [1]
-   Its relative error (for round-to-nearest ties-to-even) is bounded by 5u^2
-   (Theorem 2.6 of [2]), where u = 2^-53 for double precision,
-   assuming xh = RN(xh + xl), which implies |xl| <= 1/2 ulp(xh),
-   and similarly for ch, cl. */
-static inline double
-muldd (double xh, double xl, double ch, double cl, double *l)
-{
-  double ahlh = ch * xl, alhh = cl * xh, ahhh = ch * xh,
-        ahhl = fma (ch, xh, -ahhh);
-  ahhl += alhh + ahlh;
-  ch = ahhh + ahhl;
-  *l = (ahhh - ch) + ahhl;
-  return ch;
-}
-
-static inline double
-mulddd (double xh, double xl, double ch, double *l)
-{
-  double ahlh = ch * xl, ahhh = ch * xh, ahhl = fma (ch, xh, -ahhh);
-  ahhl += ahlh;
-  ch = ahhh + ahhl;
-  *l = (ahhh - ch) + ahhl;
-  return ch;
-}
-
-static inline double
-polydd (double xh, double xl, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double ch = c[i][0] + *l, cl = ((c[i][0] - ch) + *l) + c[i][1];
-  while (--i >= 0)
-    {
-      ch = muldd (xh, xl, ch, cl, &cl);
-      double th = ch + c[i][0], tl = (c[i][0] - th) + ch;
-      ch = th;
-      cl += tl + c[i][1];
-    }
-  *l = cl;
-  return ch;
-}
+#include "ddcoremath.h"
 
 static double __attribute__ ((noinline)) as_acosh_refine (double, double);
 static double __attribute__ ((noinline))

diff --git a/sysdeps/ieee754/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c
index 471535bd3e555c57e215dc4f976db6db5329cc4b..9ebc7bab8471e828db87ec5e8c4c13fb4a8585cb 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_atanh.c
+++ b/sysdeps/ieee754/dbl-64/e_atanh.c
@@ -29,67 +29,7 @@ SOFTWARE. */
 #include <libm-alias-finite.h>
 #include "math_config.h"
 #include "s_atanh_data.h"
-
-static inline double
-fasttwosum (double x, double y, double *e)
-{
-  double s = x + y, z = s - x;
-  *e = y - z;
-  return s;
-}
-
-static inline double
-fasttwosub (double x, double y, double *e)
-{
-  double s = x - y, z = x - s;
-  *e = z - y;
-  return s;
-}
-
-static inline double
-adddd (double xh, double xl, double ch, double cl, double *l)
-{
-  double s = xh + ch, d = s - xh;
-  *l = ((ch - d) + (xh + (d - s))) + (xl + cl);
-  return s;
-}
-
-static inline double
-muldd (double xh, double xl, double ch, double cl, double *l)
-{
-  double ahlh = ch * xl, alhh = cl * xh, ahhh = ch * xh,
-        ahhl = fma (ch, xh, -ahhh);
-  ahhl += alhh + ahlh;
-  ch = ahhh + ahhl;
-  *l = (ahhh - ch) + ahhl;
-  return ch;
-}
-
-static inline double
-mulddd (double xh, double xl, double ch, double *l)
-{
-  double ahlh = ch * xl, ahhh = ch * xh, ahhl = fma (ch, xh, -ahhh);
-  ahhl += ahlh;
-  ch = ahhh + ahhl;
-  *l = (ahhh - ch) + ahhl;
-  return ch;
-}
-
-static inline double
-polydd (double xh, double xl, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double ch = c[i][0] + *l, cl = ((c[i][0] - ch) + *l) + c[i][1];
-  while (--i >= 0)
-    {
-      ch = muldd (xh, xl, ch, cl, &cl);
-      double th = ch + c[i][0], tl = (c[i][0] - th) + ch;
-      ch = th;
-      cl += tl + c[i][1];
-    }
-  *l = cl;
-  return ch;
-}
+#include "ddcoremath.h"
 
 static double __attribute__ ((noinline)) as_atanh_refine (double, double,
                                                          double, double);

diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c
index fb241897ab60b50f2e787c10425937c0934d0e58..997c63a7dfe66a9a3d0ffb88c0559a5c6298e472 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_gamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c
@@ -29,110 +29,7 @@ SOFTWARE.
 #include <float.h>
 #include <libm-alias-finite.h>
 #include "math_config.h"
-
-static inline double
-fasttwosum (double x, double y, double *e)
-{
-  double s = x + y, z = s - x;
-  *e = y - z;
-  return s;
-}
-
-static inline double
-fastsum (double xh, double xl, double yh, double yl, double *e)
-{
-  double sl, sh = fasttwosum (xh, yh, &sl);
-  *e = (xl + yl) + sl;
-  return sh;
-}
-
-static inline double
-sumdd (double xh, double xl, double yh, double yl, double *e)
-{
-  double sl, sh;
-  if (__glibc_likely (fabs (xh) > fabs (yh)))
-    sh = fasttwosum (xh, yh, &sl);
-  else
-    sh = fasttwosum (yh, xh, &sl);
-  *e = (xl + yl) + sl;
-  return sh;
-}
-
-static inline double
-twosum (double x, double y, double *e)
-{
-  if (__glibc_likely (fabs (x) > fabs (y)))
-    return fasttwosum (x, y, e);
-  else
-    return fasttwosum (y, x, e);
-}
-
-static inline double
-muldd (double xh, double xl, double ch, double cl, double *l)
-{
-  double ahhh = ch * xh;
-  *l = (ch * xl + cl * xh) + fma (ch, xh, -ahhh);
-  return ahhh;
-}
-
-static inline double
-muldd3 (double xh, double xl, double yh, double yl, double *l)
-{
-  double ch = xh * yh, cl1 = fma (xh, yh, -ch);
-  double tl0 = xl * yl;
-  double tl1 = tl0 + xh * yl;
-  double cl2 = tl1 + xl * yh;
-  double cl3 = cl1 + cl2;
-  return fasttwosum (ch, cl3, l);
-}
-
-static inline double
-mulddd (double x, double ch, double cl, double *l)
-{
-  double ahhh = ch * x;
-  *l = cl * x + fma (ch, x, -ahhh);
-  return ahhh;
-}
-
-static inline double
-polydd (double xh, double xl, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double cl, ch = fasttwosum (c[i][0], *l, &cl);
-  cl += c[i][1];
-  while (--i >= 0)
-    {
-      ch = muldd (xh, xl, ch, cl, &cl);
-      ch = fastsum (c[i][0], c[i][1], ch, cl, &cl);
-    }
-  *l = cl;
-  return ch;
-}
-
-static inline double
-polyddd (double x, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double cl, ch = fasttwosum (c[i][0], *l, &cl);
-  cl += c[i][1];
-  while (--i >= 0)
-    {
-      ch = mulddd (x, ch, cl, &cl);
-      ch = sumdd (c[i][0], c[i][1], ch, cl, &cl);
-    }
-  *l = cl;
-  return ch;
-}
-
-static inline double
-polyd (double x, int n, const double c[][2])
-{
-  int i = n - 1;
-  double ch = c[i][0];
-  while (--i >= 0)
-    ch = c[i][0] + x * ch;
-  return ch;
-}
+#include "ddcoremath.h"
 
 static double __attribute__ ((noinline)) as_logd (double, double *);
 static double __attribute__ ((noinline)) as_expd (double, double *, int *);
@@ -173,7 +70,7 @@ poly3 (double d, unsigned imax, const double ch[], unsigned jmax,
     {
       t0 = sprod (d, t0, t1, t2, &t1, &t2);
       s0 += t0 * ch[i];
-      double fl, fh = mulddd (ch[i], t1, t2, &fl);
+      double fl, fh = mulddd2 (ch[i], t1, t2, &fl);
       s1 = sumdd (s1, s2, fh, fl, &s2);
     }
   double fl = 0, fh = polyddd (d, jmax, cl, &fl);
@@ -884,7 +781,7 @@ __ieee754_gamma_r (double x, int *signgamp)
        return __math_invalid (0);
       double t0h = 1, t0l = 0, x0 = 1;
       for (int i = 1; i < k; i++, x0 += 1.0)
-       t0h = mulddd (x0, t0h, t0l, &t0l);
+       t0h = mulddd2 (x0, t0h, t0l, &t0l);
       return t0h + t0l;
     }
 
@@ -904,7 +801,7 @@ __ieee754_gamma_r (double x, int *signgamp)
       double ix = floor (x), dx = x - ix;
       int ip = ix;
       double sl, sh = as_sinpid (dx, &sl);
-      lh = muldd (sh, sl, lh, ll, &ll);
+      lh = muldd2 (sh, sl, lh, ll, &ll);
       const double pih = 0x1.921fb54442d18p+1, pil = 0x1.1a62633145c07p-53;
       double rcp = 1 / lh, rh = rcp * pih,
             rl = rcp * (pil - ll * rh - fma (rh, lh, -pih));
@@ -1015,11 +912,11 @@ __ieee754_gamma_r (double x, int *signgamp)
          else
            xph = fasttwosum (1, xph, &l);
          xpl += l;
-         wh = muldd (xph, xpl, wh, wl, &wl);
+         wh = muldd2 (xph, xpl, wh, wl, &wl);
        }
     }
   double rh = 1.0 / wh, rl = (fma (rh, -wh, 1.0) - wl * rh) * rh;
-  fh = muldd (rh, rl, fh, fl, &fl);
+  fh = muldd2 (rh, rl, fh, fl, &fl);
   double eps = fh * 0x1.2e3b40a0e9b4fp-70;
   double ub = fh + (fl + eps), lb = fh + (fl - eps);
   if (ub != lb)
@@ -1242,13 +1139,13 @@ as_sinpid (double x, double *l)
   double Q = d2 * (s[1] + d2 * (s[2] + d2 * s[3]));
   double ql, qh = fasttwosum (s[0], Q, &ql);
   ql += s0;
-  ch = muldd (qh, ql, ch, cl, &cl);
+  ch = muldd2 (qh, ql, ch, cl, &cl);
   double tl, th = fasttwosum (c[0], P, &tl);
   tl += c0;
-  th = mulddd (d, th, tl, &tl);
-  double pl, ph = muldd (th, tl, sh, sl, &pl);
+  th = mulddd2 (d, th, tl, &tl);
+  double pl, ph = muldd2 (th, tl, sh, sl, &pl);
   ch = fastsum (ch, cl, ph, pl, &cl);
-  ch = mulddd (d, ch, cl, &cl);
+  ch = mulddd2 (d, ch, cl, &cl);
   sh = fastsum (sh, sl, ch, cl, l);
   return sh;
 }
@@ -1325,12 +1222,12 @@ as_expd (double x, double *l, int *e)
 {
   const double ln2h = 0x1.71547652b82fep+10, ln2l = 0x1.777d0ffda0d24p-46;
   double xh = x, xl = *l;
-  xh = muldd (xh, xl, ln2h, ln2l, &xl);
+  xh = muldd2 (xh, xl, ln2h, ln2l, &xl);
   double ix = roundeven_finite (xh);
   xh = fasttwosum (xh - ix, xl, &xl);
   int k = ix, i0 = (k >> 5) & 31, i1 = k & 31;
   *e = k >> 10;
-  double rl, rh = muldd (E0[i0][1], E0[i0][0], E1[i1][1], E1[i1][0], &rl);
+  double rl, rh = muldd2 (E0[i0][1], E0[i0][0], E1[i1][1], E1[i1][0], &rl);
   static const double c[][2]
       = { { 0x1.62e42fefa39efp-11, 0x1.abc9e3bf9d4d1p-66 },
          { 0x1.ebfbdff82c58ep-23, 0x1.ec07243b4e585p-77 },
@@ -1340,11 +1237,11 @@ as_expd (double x, double *l, int *e)
   const int m = 1;
   double fh, fl, el;
   fl = xh * polyd (xh, 5 - m, c + m);
-  fh = polydd (xh, xl, m, c, &fl);
-  fh = muldd (xh, xl, fh, fl, &fl);
+  fh = polydd2 (xh, xl, m, c, &fl);
+  fh = muldd2 (xh, xl, fh, fl, &fl);
   fh = fasttwosum (1, fh, &el);
   fl += el;
-  rh = muldd (rh, rl, fh, fl, &rl);
+  rh = muldd2 (rh, rl, fh, fl, &rl);
   *l = rl;
   return rh;
 }
@@ -1355,8 +1252,8 @@ as_lgamma_asym (double xh, double *xl)
   double zh = 1.0 / xh, dz = *xl * zh, zl = (fma (zh, -xh, 1.0) - dz) * zh;
   double ll, lh = as_logd (xh, &ll);
   ll += dz;
-  lh = muldd (xh - 0.5, *xl, lh - 1, ll, &ll);
-  double z2l, z2h = muldd (zh, zl, zh, zl, &z2l);
+  lh = muldd2 (xh - 0.5, *xl, lh - 1, ll, &ll);
+  double z2l, z2h = muldd2 (zh, zl, zh, zl, &z2l);
   double fh, fl;
   double x2 = z2h * z2h;
   if (xh > 11.5)
@@ -1377,7 +1274,7 @@ as_lgamma_asym (double xh, double *xl)
       double q2 = q[2][0] + z2h * q[3][0];
       double q4 = q[4][0] + z2h * q[5][0];
       fl = z2h * (q0 + x2 * (q2 + x2 * q4));
-      fh = polydd (z2h, z2l, k, b, &fl);
+      fh = polydd2 (z2h, z2l, k, b, &fl);
     }
   else
     {
@@ -1408,8 +1305,8 @@ as_lgamma_asym (double xh, double *xl)
       q0 += x2 * q2;
       q0 += x4 * q4;
       fl = z2h * q0;
-      fh = polydd (z2h, z2l, k, b, &fl);
+      fh = polydd2 (z2h, z2l, k, b, &fl);
     }
-  fh = muldd (zh, zl, fh, fl, &fl);
+  fh = muldd2 (zh, zl, fh, fl, &fl);
   return fastsum (lh, ll, fh, fl, xl);
 }

diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
index c630e0913fda867a97aff8affbd5d91276a7fdad..aab7b4d1443b9146f51b73bcf167bd4f16360ea7 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
@@ -31,101 +31,7 @@ SOFTWARE.
 #include <math.h>
 #include <libm-alias-finite.h>
 #include "math_config.h"
-
-static inline double
-fasttwosum (double x, double y, double *e)
-{
-  double s = x + y, z = s - x;
-  *e = y - z;
-  return s;
-}
-
-static inline double
-twosum (double x, double y, double *e)
-{
-  if (__glibc_likely (fabs (x) > fabs (y)))
-    return fasttwosum (x, y, e);
-  else
-    return fasttwosum (y, x, e);
-}
-
-static inline double
-fastsum (double xh, double xl, double yh, double yl, double *e)
-{
-  double sl, sh = fasttwosum (xh, yh, &sl);
-  *e = (xl + yl) + sl;
-  return sh;
-}
-
-static inline double
-sumdd (double xh, double xl, double yh, double yl, double *e)
-{
-  double sl, sh;
-  char o = fabs (xh) > fabs (yh);
-  if (__glibc_likely (o))
-    sh = fasttwosum (xh, yh, &sl);
-  else
-    sh = fasttwosum (yh, xh, &sl);
-  sl += xl + yl;
-  *e = sl;
-  return sh;
-}
-
-static inline double
-muldd (double xh, double xl, double ch, double cl, double *l)
-{
-  double ahhh = ch * xh;
-  *l = (ch * xl + cl * xh) + fma (ch, xh, -ahhh);
-  return ahhh;
-}
-
-static inline double
-mulddd (double x, double ch, double cl, double *l)
-{
-  double ahhh = ch * x;
-  *l = cl * x + fma (ch, x, -ahhh);
-  return ahhh;
-}
-
-static inline double
-polydd (double xh, double xl, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double cl, ch = fasttwosum (c[i][0], *l, &cl);
-  cl += c[i][1];
-  while (--i >= 0)
-    {
-      ch = muldd (xh, xl, ch, cl, &cl);
-      ch = fastsum (c[i][0], c[i][1], ch, cl, &cl);
-    }
-  *l = cl;
-  return ch;
-}
-
-static inline double
-polydddfst (double x, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double cl, ch = fasttwosum (c[i][0], *l, &cl);
-  cl += c[i][1];
-  while (--i >= 0)
-    {
-      ch = mulddd (x, ch, cl, &cl);
-      ch = fastsum (c[i][0], c[i][1], ch, cl, &cl);
-    }
-  *l = cl;
-  return ch;
-}
-
-static inline double
-polyd (double x, int n, const double c[][2])
-{
-  int i = n - 1;
-  double ch = c[i][0];
-  while (--i >= 0)
-    ch = c[i][0] + x * ch;
-  return ch;
-}
+#include "ddcoremath.h"
 
 static double __attribute__ ((noinline)) as_logd (double, double *);
 static double __attribute__ ((noinline)) as_logd_accurate (double, double *,
@@ -268,7 +174,7 @@ as_lgamma_accurate (double x)
   else if (x < 0x1p-2)
     {
       fh = polydddfst (sx, array_length (c0), c0, &fl);
-      fh = mulddd (sx, fh, fl, &fl);
+      fh = mulddd2 (sx, fh, fl, &fl);
       double lll, ll, lh = as_logd_accurate (x, &ll, &lll);
       fh = sumdd (fh, fl, -ll, -lll, &fl);
       fh = twosum (-lh, fh, &fll);
@@ -328,7 +234,7 @@ as_lgamma_accurate (double x)
       else if (fabs (x - 1) < 0x1p-2)
        {
          fh = polydddfst (x - 1, array_length (c0), c0, &fl);
-         fh = mulddd (x - 1, fh, fl, &fl);
+         fh = mulddd2 (x - 1, fh, fl, &fl);
          if (sx < 0)
            {
              double lll, ll, lh = as_logd_accurate (x, &ll, &lll);
@@ -352,7 +258,7 @@ as_lgamma_accurate (double x)
        {
          double lll, ll, lh = as_logd_accurate (x - 1, &ll, &lll);
          fh = polydddfst (x - 2, array_length (c0), c0, &fl);
-         fh = mulddd (x - 2, fh, fl, &fl);
+         fh = mulddd2 (x - 2, fh, fl, &fl);
          fl = sumdd (fl, 0, ll, lll, &fll);
          fh = twosum (fh, lh, &lh);
          fl = sumdd (fl, fll, lh, 0, &fll);
@@ -372,7 +278,7 @@ as_lgamma_accurate (double x)
          l1h = fasttwosum (l1h, l2h, &l2h);
          l1l = sumdd (l1l, l1ll, l2h, 0, &l1ll);
          fh = polydddfst (x - 3, array_length (c0), c0, &fl);
-         fh = mulddd (x - 3, fh, fl, &fl);
+         fh = mulddd2 (x - 3, fh, fl, &fl);
          fl = sumdd (fl, 0, l1l, l1ll, &fll);
          fh = twosum (fh, l1h, &l1h);
          fl = sumdd (fl, fll, l1h, 0, &fll);
@@ -458,8 +364,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.168p-4 && sx > -3 && sx < SX_BND)
        {
@@ -500,8 +406,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.2d4p-4 && sx > -3.5 && sx < -3)
        {
@@ -535,8 +441,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 7;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.efp-5 && sx > -4 && sx < -3.5)
        {
@@ -570,8 +476,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 7;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1p-3 && sx > -4.5 && sx < -4)
        {
@@ -609,8 +515,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 7;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.08p-4 && sx > -5 && sx < -4.5)
        {
@@ -643,8 +549,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 7;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.1p-4 && sx > -5.5 && sx < -5)
        {
@@ -677,8 +583,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.13p-3 && sx > -6 && sx < -5.5)
        {
@@ -716,8 +622,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.15p-3 && sx > -6.5 && sx < -6)
        {
@@ -755,8 +661,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.34p-3 && sx > -7.0 && sx < -6.5)
        {
@@ -795,8 +701,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.145p-3 && sx > -7.5 && sx < -7)
        {
@@ -834,8 +740,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.76bp-3 && sx > -8 && sx < -7.5)
        {
@@ -876,8 +782,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.76cp-3 && sx > -0x1.10f5c28f5c28fp+3
               && sx < -8)
@@ -919,8 +825,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.99p-3 && sx > -9 && sx < -8.5)
        {
@@ -962,8 +868,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.99p-3 && sx > -9.5 && sx < -9)
        {
@@ -1005,8 +911,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.76bp-3 && sx > -10 && sx < -9.5)
        {
@@ -1047,8 +953,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
       else if (fabs (fh) < 0x1.76bp-3 && sx > -10.5 && sx < -10)
        {
@@ -1089,8 +995,8 @@ as_lgamma_accurate (double x)
          double sh = zh * sc, sl = zl * sc;
          int n = array_length (c), k = 6;
          fl = sh * polyd (sh, k, c + n - k);
-         fh = polydd (sh, sl, n - k, c, &fl);
-         fh = muldd (zh, zl, fh, fl, &fl);
+         fh = polydd2 (sh, sl, n - k, c, &fl);
+         fh = muldd2 (zh, zl, fh, fl, &fl);
        }
     }
 
@@ -1344,7 +1250,7 @@ __ieee754_lgamma_r (double x, int *signgamp)
                 q4 = q[4] + z * q[5], q6 = q[6] + z * q[7];
          fl = z * ((q0 + z2 * q2) + z4 * (q4 + z2 * q6));
          fh = polydddfst (z, 4, c0, &fl);
-         fh = mulddd (x, fh, fl, &fl);
+         fh = mulddd2 (x, fh, fl, &fl);
          fh = sumdd (-lh, -ll, fh, fl, &fl);
          eps = 1.5e-22;
        }
@@ -1367,7 +1273,7 @@ __ieee754_lgamma_r (double x, int *signgamp)
          if (__glibc_unlikely (j == 4))
            { // treat the region around the root at 1
              z = -x;
-             fh = mulddd (z, fh, fl, &fl);
+             fh = mulddd2 (z, fh, fl, &fl);
            }
          eps = fabs (fh) * 8.3e-20;
          fh = sumdd (-lh, -ll, fh, fl, &fl);
@@ -1388,13 +1294,13 @@ __ieee754_lgamma_r (double x, int *signgamp)
              if (__glibc_unlikely (au >= 0x3fabaa6))
                lh = fasttwosum (lh, ll, &ll); // x>=0x1.754cp+1014
              double hlh = lh * 0.5;
-             lh = mulddd (ax, lh, ll, &ll);
+             lh = mulddd2 (ax, lh, ll, &ll);
              ll -= hlh;
            }
          else
            {
              // for other |x| use a simple product
-             lh = mulddd (ax - 0.5, lh, ll, &ll);
+             lh = mulddd2 (ax - 0.5, lh, ll, &ll);
            }
          static const double c[][2]
              = { { 0x1.acfe390c97d6ap-2, -0x1.1d9792ced423ap-58 },
@@ -1412,7 +1318,7 @@ __ieee754_lgamma_r (double x, int *signgamp)
              fl = z2h * (q0 + z4h * (q2 + z4h * q4));
              fh = fasttwosum (c[1][0], fl, &fl);
              fl += c[1][1];
-             fh = muldd (fh, fl, zh, zl, &fl);
+             fh = muldd2 (fh, fl, zh, zl, &fl);
            }
          else
            {
@@ -1436,19 +1342,19 @@ __ieee754_lgamma_r (double x, int *signgamp)
          if (__glibc_unlikely (j == 4))
            { // treat the region around the root at 1
              z = 1 - ax;
-             fh = mulddd (z, fh, fl, &fl);
+             fh = mulddd2 (z, fh, fl, &fl);
            }
          if (__glibc_unlikely (j == 10))
            { // treat the region around the root at 2
              z = ax - 2;
-             fh = mulddd (z, fh, fl, &fl);
+             fh = mulddd2 (z, fh, fl, &fl);
            }
          eps = fabs (fh) * 8.3e-20 + 1e-24;
        }
       if (t >> 63)
        { // x<0 so use reflection formula
          double sl, sh = as_sinpipid (x - floor (x), &sl);
-         sh = mulddd (-x, sh, sl, &sl);
+         sh = mulddd2 (-x, sh, sl, &sl);
          double ll, lh = as_logd (sh, &ll);
          ll += sl / sh;
          fh = -sumdd (fh, fl, lh, ll, &fl);
@@ -1741,8 +1647,8 @@ as_logd_accurate (double x, double *l, double *l_)
          { 0x1.c71e5ec7051f6p-4, 0x1.217ec3dcb2f03p-58 } };
   dxh = fasttwosum (dxh, dxl, &dxl);
   double fl = dxh * (c[6][0] + dxh * (c[7][0] + dxh * (c[8][0]))),
-        fh = polydd (dxh, dxl, 6, c, &fl);
-  fh = muldd (dxh, dxl, fh, fl, &fl);
+        fh = polydd2 (dxh, dxl, 6, c, &fl);
+  fh = muldd2 (dxh, dxl, fh, fl, &fl);
   double s2 = h1[i1][2] + h2[i2][2], s1 = h1[i1][1] + h2[i2][1],
         s0 = h1[i1][0] + h2[i2][0];
   double L0 = 0x1.62e42fefa38p-1 * ed, L1 = 0x1.ef35793c76p-45 * ed,
@@ -1845,8 +1751,8 @@ as_sinpipid (double x, double *l)
       double fl = z2 * (cl[0] + z2 * (cl[1] + z2 * (cl[2]))),
             fh = fasttwosum (c[0], fl, &fl), e;
       fl += c[1];
-      fh = muldd (z2, z2l, fh, fl, &fl);
-      fh = mulddd (z, fh, fl, &fl);
+      fh = muldd2 (z2, z2l, fh, fl, &fl);
+      fh = mulddd2 (z, fh, fl, &fl);
       fh = fasttwosum (z, fh, &e);
       fl += e;
       *l = fl;
@@ -1872,13 +1778,13 @@ as_sinpipid (double x, double *l)
 
   double ql, qh = fasttwosum (s[0], Q, &ql);
   ql += s0;
-  ch = muldd (qh, ql, ch, cl, &cl);
+  ch = muldd2 (qh, ql, ch, cl, &cl);
   double tl, th = fasttwosum (c[0], P, &tl);
   tl += c0;
-  th = mulddd (d, th, tl, &tl);
-  double pl, ph = muldd (th, tl, sh, sl, &pl);
+  th = mulddd2 (d, th, tl, &tl);
+  double pl, ph = muldd2 (th, tl, sh, sl, &pl);
   ch = fastsum (ch, cl, ph, pl, &cl);
-  ch = mulddd (d, ch, cl, &cl);
+  ch = mulddd2 (d, ch, cl, &cl);
   sh = fastsum (sh, sl, ch, cl, l);
   return sh;
 }
@@ -1912,15 +1818,15 @@ as_sinpipid_accurate (double x, double *l)
          { -0x1.e306ec8cf7c02p-85, 0x1.5d2601f85289ap-139 } };
 
   double d2h = d * d, d2l = fma (d, d, -d2h);
-  double Pl = 0, Ph = polydd (d2h, d2l, 5, c, &Pl);
-  double Ql = 0, Qh = polydd (d2h, d2l, 6, s, &Ql);
+  double Pl = 0, Ph = polydd2 (d2h, d2l, 5, c, &Pl);
+  double Ql = 0, Qh = polydd2 (d2h, d2l, 6, s, &Ql);
 
-  Ph = mulddd (d, Ph, Pl, &Pl);
-  Ph = muldd (sh, sl, Ph, Pl, &Pl);
-  Qh = muldd (ch, cl, Qh, Ql, &Ql);
+  Ph = mulddd2 (d, Ph, Pl, &Pl);
+  Ph = muldd2 (sh, sl, Ph, Pl, &Pl);
+  Qh = muldd2 (ch, cl, Qh, Ql, &Ql);
 
   ch = fastsum (Qh, Ql, Ph, Pl, &cl);
-  ch = mulddd (d, ch, cl, &cl);
+  ch = mulddd2 (d, ch, cl, &cl);
   sh = fastsum (sh, sl, ch, cl, l);
   return sh;
 }
@@ -1939,7 +1845,7 @@ as_lgamma_asym_accurate (double xh, double *xl, double *e)
             zl = (fma (zh, -xh, 1.0) - dz) * zh;
       if (*xl != 0)
        {
-         double dl2, dl1 = mulddd (*xl, zh, fma (zh, -xh, 1.0) * zh, &dl2);
+         double dl2, dl1 = mulddd2 (*xl, zh, fma (zh, -xh, 1.0) * zh, &dl2);
          dl2 -= dl1 * dl1 / 2;
          l1 = sumdd (l1, l2, dl1, dl2, &l2);
        }
@@ -1965,7 +1871,7 @@ as_lgamma_asym_accurate (double xh, double *xl, double *e)
 
       l1x = sumdd (l1x, l2x, l0xl, l1xl, &l2x);
       l1x = sumdd (l1x, l2x, l0 * wl, l1 * wl, &l2x);
-      double z2l, z2h = muldd (zh, zl, zh, zl, &z2l);
+      double z2l, z2h = muldd2 (zh, zl, zh, zl, &z2l);
       double fh, fl;
       if (xh >= 48)
        {
@@ -1980,7 +1886,7 @@ as_lgamma_asym_accurate (double xh, double *xl, double *e)
                  { 0x1.a0a6926f4992p-8, -0x1.1f355cbf82229p-63 } };
          l1x = sumdd (l1x, l2x, c[0][0], c[0][1], &l2x);
          fl = 0;
-         fh = polydd (z2h, z2l, 7, c + 1, &fl);
+         fh = polydd2 (z2h, z2l, 7, c + 1, &fl);
        }
       else if (xh >= 14.5)
        {
@@ -1999,7 +1905,7 @@ as_lgamma_asym_accurate (double xh, double *xl, double *e)
                  { 0x1.2ea102098f818p+3, 0x1.609db97f1bc89p-51 } };
          l1x = sumdd (l1x, l2x, c[0][0], c[0][1], &l2x);
          fl = 0;
-         fh = polydd (z2h, z2l, 11, c + 1, &fl);
+         fh = polydd2 (z2h, z2l, 11, c + 1, &fl);
        }
       else
        {
@@ -2034,9 +1940,9 @@ as_lgamma_asym_accurate (double xh, double *xl, double *e)
                  { 0x1.59bad61bd81d5p+47, 0x1.0e3f2ea42a0ep-7 } };
          l1x = sumdd (l1x, l2x, c[0][0], c[0][1], &l2x);
          fl = 0;
-         fh = polydd (z2h, z2l, array_length (c) - 1, c + 1, &fl);
+         fh = polydd2 (z2h, z2l, array_length (c) - 1, c + 1, &fl);
        }
-      fh = muldd (zh, zl, fh, fl, &fl);
+      fh = muldd2 (zh, zl, fh, fl, &fl);
       l1x = sumdd (l1x, l2x, fh, fl, &l2x);
       l0x = fasttwosum (l0x, l1x, &l1x);
       l1x = fasttwosum (l1x, l2x, &l2x);

diff --git a/sysdeps/ieee754/dbl-64/s_asinh.c b/sysdeps/ieee754/dbl-64/s_asinh.c
index e19ae4e0a01a90223e69e36859d73226b6b9d359..fca50433c625229978bda4ababf1663824b51713 100644 (file)
--- a/sysdeps/ieee754/dbl-64/s_asinh.c
+++ b/sysdeps/ieee754/dbl-64/s_asinh.c
@@ -30,59 +30,7 @@ SOFTWARE.
 #include <libm-alias-double.h>
 #include "math_config.h"
 #include "s_asincosh_data.h"
-
-static inline double
-fasttwosum (double x, double y, double *e)
-{
-  double s = x + y, z = s - x;
-  *e = y - z;
-  return s;
-}
-
-static inline double
-adddd (double xh, double xl, double ch, double cl, double *l)
-{
-  double s = xh + ch, d = s - xh;
-  *l = ((ch - d) + (xh + (d - s))) + (xl + cl);
-  return s;
-}
-
-static inline double
-muldd (double xh, double xl, double ch, double cl, double *l)
-{
-  double ahlh = ch * xl, alhh = cl * xh, ahhh = ch * xh,
-        ahhl = fma (ch, xh, -ahhh);
-  ahhl += alhh + ahlh;
-  ch = ahhh + ahhl;
-  *l = (ahhh - ch) + ahhl;
-  return ch;
-}
-
-static inline double
-mulddd (double xh, double xl, double ch, double *l)
-{
-  double ahlh = ch * xl, ahhh = ch * xh, ahhl = fma (ch, xh, -ahhh);
-  ahhl += ahlh;
-  ch = ahhh + ahhl;
-  *l = (ahhh - ch) + ahhl;
-  return ch;
-}
-
-static inline double
-polydd (double xh, double xl, int n, const double c[][2], double *l)
-{
-  int i = n - 1;
-  double ch = c[i][0] + *l, cl = ((c[i][0] - ch) + *l) + c[i][1];
-  while (--i >= 0)
-    {
-      ch = muldd (xh, xl, ch, cl, &cl);
-      double th = ch + c[i][0], tl = (c[i][0] - th) + ch;
-      ch = th;
-      cl += tl + c[i][1];
-    }
-  *l = cl;
-  return ch;
-}
+#include "ddcoremath.h"
 
 static double __attribute__ ((noinline)) as_asinh_refine (double, double,
                                                          double, double);