Merge changes between r22241 and r22552 from /fsf/trunk.

git-svn-id: svn://svn.eglibc.org/trunk@22553 7b3dc134-2b1b-0410-93df-9e9f96275f8d

13 files changed, 578 insertions, 334 deletions

diff --git a/libc/sysdeps/ieee754/bits/nan.h b/libc/sysdeps/ieee754/bits/nan.h
index d3ab38ba7..935271a7c 100644
--- a/libc/sysdeps/ieee754/bits/nan.h
+++ b/libc/sysdeps/ieee754/bits/nan.h
@@ -46,7 +46,7 @@
 # endif
 
 static union { unsigned char __c[4]; float __d; } __nan_union
-    __attribute_used__ = { __nan_bytes };
+  __attribute__ ((__unused__)) = { __nan_bytes };
 # define NAN	(__nan_union.__d)
 
 #endif	/* GCC.  */
diff --git a/libc/sysdeps/ieee754/dbl-64/mpa.c b/libc/sysdeps/ieee754/dbl-64/mpa.c
index ede8ed198..8fc2626f7 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpa.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpa.c
@@ -43,6 +43,7 @@
 #include "endian.h"
 #include "mpa.h"
 #include <sys/param.h>
+#include <alloca.h>
 
 #ifndef SECTION
 # define SECTION
@@ -59,8 +60,9 @@ const mp_no mptwo = {1, {1.0, 2.0}};
 static int
 mcr (const mp_no *x, const mp_no *y, int p)
 {
-  int i;
-  for (i = 1; i <= p; i++)
+  long i;
+  long p2 = p;
+  for (i = 1; i <= p2; i++)
     {
       if (X[i] == Y[i])
 	continue;
@@ -76,7 +78,7 @@ mcr (const mp_no *x, const mp_no *y, int p)
 int
 __acr (const mp_no *x, const mp_no *y, int p)
 {
-  int i;
+  long i;
 
   if (X[0] == ZERO)
     {
@@ -107,8 +109,10 @@ __acr (const mp_no *x, const mp_no *y, int p)
 void
 __cpy (const mp_no *x, mp_no *y, int p)
 {
+  long i;
+
   EY = EX;
-  for (int i = 0; i <= p; i++)
+  for (i = 0; i <= p; i++)
     Y[i] = X[i];
 }
 #endif
@@ -119,8 +123,8 @@ __cpy (const mp_no *x, mp_no *y, int p)
 static void
 norm (const mp_no *x, double *y, int p)
 {
-#define R  RADIXI
-  int i;
+#define R RADIXI
+  long i;
   double a, c, u, v, z[5];
   if (p < 5)
     {
@@ -194,17 +198,18 @@ norm (const mp_no *x, double *y, int p)
 static void
 denorm (const mp_no *x, double *y, int p)
 {
-  int i, k;
+  long i, k;
+  long p2 = p;
   double c, u, z[5];
 
-#define R  RADIXI
+#define R RADIXI
   if (EX < -44 || (EX == -44 && X[1] < TWO5))
     {
       *y = ZERO;
       return;
     }
 
-  if (p == 1)
+  if (p2 == 1)
     {
       if (EX == -42)
 	{
@@ -228,7 +233,7 @@ denorm (const mp_no *x, double *y, int p)
 	  k = 1;
 	}
     }
-  else if (p == 2)
+  else if (p2 == 2)
     {
       if (EX == -42)
 	{
@@ -281,7 +286,7 @@ denorm (const mp_no *x, double *y, int p)
 
   if (u == z[3])
     {
-      for (i = k + 1; i <= p; i++)
+      for (i = k + 1; i <= p2; i++)
 	{
 	  if (X[i] == ZERO)
 	    continue;
@@ -323,7 +328,8 @@ void
 SECTION
 __dbl_mp (double x, mp_no *y, int p)
 {
-  int i, n;
+  long i, n;
+  long p2 = p;
   double u;
 
   /* Sign.  */
@@ -347,7 +353,7 @@ __dbl_mp (double x, mp_no *y, int p)
     x *= RADIX;
 
   /* Digits.  */
-  n = MIN (p, 4);
+  n = MIN (p2, 4);
   for (i = 1; i <= n; i++)
     {
       u = (x + TWO52) - TWO52;
@@ -357,7 +363,7 @@ __dbl_mp (double x, mp_no *y, int p)
       x -= u;
       x *= RADIX;
     }
-  for (; i <= p; i++)
+  for (; i <= p2; i++)
     Y[i] = ZERO;
 }
 
@@ -369,53 +375,64 @@ static void
 SECTION
 add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  int i, j, k;
+  long i, j, k;
+  long p2 = p;
+  double zk;
 
   EZ = EX;
 
-  i = p;
-  j = p + EY - EX;
-  k = p + 1;
+  i = p2;
+  j = p2 + EY - EX;
+  k = p2 + 1;
 
-  if (j < 1)
+  if (__glibc_unlikely (j < 1))
     {
       __cpy (x, z, p);
       return;
     }
-  else
-    Z[k] = ZERO;
+
+  zk = ZERO;
 
   for (; j > 0; i--, j--)
     {
-      Z[k] += X[i] + Y[j];
-      if (Z[k] >= RADIX)
+      zk += X[i] + Y[j];
+      if (zk >= RADIX)
 	{
-	  Z[k] -= RADIX;
-	  Z[--k] = ONE;
+	  Z[k--] = zk - RADIX;
+	  zk = ONE;
 	}
       else
-	Z[--k] = ZERO;
+        {
+	  Z[k--] = zk;
+	  zk = ZERO;
+	}
     }
 
   for (; i > 0; i--)
     {
-      Z[k] += X[i];
-      if (Z[k] >= RADIX)
+      zk += X[i];
+      if (zk >= RADIX)
 	{
-	  Z[k] -= RADIX;
-	  Z[--k] = ONE;
+	  Z[k--] = zk - RADIX;
+	  zk = ONE;
 	}
       else
-	Z[--k] = ZERO;
+        {
+	  Z[k--] = zk;
+	  zk = ZERO;
+	}
     }
 
-  if (Z[1] == ZERO)
+  if (zk == ZERO)
     {
-      for (i = 1; i <= p; i++)
+      for (i = 1; i <= p2; i++)
 	Z[i] = Z[i + 1];
     }
   else
-    EZ += ONE;
+    {
+      Z[1] = zk;
+      EZ += ONE;
+    }
 }
 
 /* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
@@ -426,71 +443,70 @@ static void
 SECTION
 sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  int i, j, k;
+  long i, j, k;
+  long p2 = p;
+  double zk;
 
   EZ = EX;
+  i = p2;
+  j = p2 + EY - EX;
+  k = p2;
 
-  if (EX == EY)
+  /* Y is too small compared to X, copy X over to the result.  */
+  if (__glibc_unlikely (j < 1))
     {
-      i = j = k = p;
-      Z[k] = Z[k + 1] = ZERO;
+      __cpy (x, z, p);
+      return;
     }
-  else
+
+  /* The relevant least significant digit in Y is non-zero, so we factor it in
+     to enhance accuracy.  */
+  if (j < p2 && Y[j + 1] > ZERO)
     {
-      j = EX - EY;
-      if (j > p)
-	{
-	  __cpy (x, z, p);
-	  return;
-	}
-      else
-	{
-	  i = p;
-	  j = p + 1 - j;
-	  k = p;
-	  if (Y[j] > ZERO)
-	    {
-	      Z[k + 1] = RADIX - Y[j--];
-	      Z[k] = MONE;
-	    }
-	  else
-	    {
-	      Z[k + 1] = ZERO;
-	      Z[k] = ZERO;
-	      j--;
-	    }
-	}
+      Z[k + 1] = RADIX - Y[j + 1];
+      zk = MONE;
     }
+  else
+    zk = Z[k + 1] = ZERO;
 
+  /* Subtract and borrow.  */
   for (; j > 0; i--, j--)
     {
-      Z[k] += (X[i] - Y[j]);
-      if (Z[k] < ZERO)
+      zk += (X[i] - Y[j]);
+      if (zk < ZERO)
 	{
-	  Z[k] += RADIX;
-	  Z[--k] = MONE;
+	  Z[k--] = zk + RADIX;
+	  zk = MONE;
 	}
       else
-	Z[--k] = ZERO;
+        {
+	  Z[k--] = zk;
+	  zk = ZERO;
+	}
     }
 
+  /* We're done with digits from Y, so it's just digits in X.  */
   for (; i > 0; i--)
     {
-      Z[k] += X[i];
-      if (Z[k] < ZERO)
+      zk += X[i];
+      if (zk < ZERO)
 	{
-	  Z[k] += RADIX;
-	  Z[--k] = MONE;
+	  Z[k--] = zk + RADIX;
+	  zk = MONE;
 	}
       else
-	Z[--k] = ZERO;
+        {
+	  Z[k--] = zk;
+	  zk = ZERO;
+	}
     }
 
+  /* Normalize.  */
   for (i = 1; Z[i] == ZERO; i++);
   EZ = EZ - i + 1;
-  for (k = 1; i <= p + 1;)
+  for (k = 1; i <= p2 + 1;)
     Z[k++] = Z[i++];
-  for (; k <= p;)
+  for (; k <= p2;)
     Z[k++] = ZERO;
 }
 
@@ -602,8 +618,11 @@ void
 SECTION
 __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  int i, j, k, k2;
+  long i, j, k, ip, ip2;
+  long p2 = p;
   double u, zk;
+  const mp_no *a;
+  double *diag;
 
   /* Is z=0?  */
   if (__glibc_unlikely (X[0] * Y[0] == ZERO))
@@ -612,48 +631,238 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
       return;
     }
 
-  /* Multiply, add and carry.  */
-  k2 = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
-  zk = Z[k2] = ZERO;
+  /* We need not iterate through all X's and Y's since it's pointless to
+     multiply zeroes.  Here, both are zero...  */
+  for (ip2 = p2; ip2 > 0; ip2--)
+    if (X[ip2] != ZERO || Y[ip2] != ZERO)
+      break;
+
+  a = X[ip2] != ZERO ? y : x;
+
+  /* ... and here, at least one of them is still zero.  */
+  for (ip = ip2; ip > 0; ip--)
+    if (a->d[ip] != ZERO)
+      break;
+
+  /* The product looks like this for p = 3 (as an example):
+
+
+				a1    a2    a3
+		 x		b1    b2    b3
+		 -----------------------------
+			     a1*b3 a2*b3 a3*b3
+		       a1*b2 a2*b2 a3*b2
+		 a1*b1 a2*b1 a3*b1
+
+     So our K needs to ideally be P*2, but we're limiting ourselves to P + 3
+     for P >= 3.  We compute the above digits in two parts; the last P-1
+     digits and then the first P digits.  The last P-1 digits are a sum of
+     products of the input digits from P to P-k where K is 0 for the least
+     significant digit and increases as we go towards the left.  The product
+     term is of the form X[k]*X[P-k] as can be seen in the above example.
+
+     The first P digits are also a sum of products with the same product term,
+     except that the sum is from 1 to k.  This is also evident from the above
+     example.
+
+     Another thing that becomes evident is that only the most significant
+     ip+ip2 digits of the result are non-zero, where ip and ip2 are the
+     'internal precision' of the input numbers, i.e. digits after ip and ip2
+     are all ZERO.  */
+
+  k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3;
+
+  while (k > ip + ip2 + 1)
+    Z[k--] = ZERO;
 
-  for (k = k2; k > p; k--)
+  zk = ZERO;
+
+  /* Precompute sums of diagonal elements so that we can directly use them
+     later.  See the next comment to know we why need them.  */
+  diag = alloca (k * sizeof (double));
+  double d = ZERO;
+  for (i = 1; i <= ip; i++)
+    {
+      d += X[i] * Y[i];
+      diag[i] = d;
+    }
+  while (i < k)
+    diag[i++] = d;
+
+  while (k > p2)
     {
-      for (i = k - p, j = p; i < p + 1; i++, j--)
-	zk += X[i] * Y[j];
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+	/* We want to add this only once, but since we subtract it in the sum
+	   of products above, we add twice.  */
+	zk += 2 * X[lim] * Y[lim];
+
+      for (i = k - p2, j = p2; i < j; i++, j--)
+	zk += (X[i] + X[j]) * (Y[i] + Y[j]);
+
+      zk -= diag[k - 1];
 
       u = (zk + CUTTER) - CUTTER;
       if (u > zk)
 	u -= RADIX;
-      Z[k] = zk - u;
+      Z[k--] = zk - u;
       zk = u * RADIXI;
     }
 
+  /* The real deal.  Mantissa digit Z[k] is the sum of all X[i] * Y[j] where i
+     goes from 1 -> k - 1 and j goes the same range in reverse.  To reduce the
+     number of multiplications, we halve the range and if k is an even number,
+     add the diagonal element X[k/2]Y[k/2].  Through the half range, we compute
+     X[i] * Y[j] as (X[i] + X[j]) * (Y[i] + Y[j]) - X[i] * Y[i] - X[j] * Y[j].
+
+     This reduction tells us that we're summing two things, the first term
+     through the half range and the negative of the sum of the product of all
+     terms of X and Y in the full range.  i.e.
+
+     SUM(X[i] * Y[i]) for k terms.  This is precalculated above for each k in
+     a single loop so that it completes in O(n) time and can hence be directly
+     used in the loop below.  */
   while (k > 1)
     {
-      for (i = 1, j = k - 1; i < k; i++, j--)
-	zk += X[i] * Y[j];
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+	/* We want to add this only once, but since we subtract it in the sum
+	   of products above, we add twice.  */
+        zk += 2 * X[lim] * Y[lim];
+
+      for (i = 1, j = k - 1; i < j; i++, j--)
+	zk += (X[i] + X[j]) * (Y[i] + Y[j]);
+
+      zk -= diag[k - 1];
 
       u = (zk + CUTTER) - CUTTER;
       if (u > zk)
 	u -= RADIX;
-      Z[k] = zk - u;
+      Z[k--] = zk - u;
       zk = u * RADIXI;
-      k--;
     }
   Z[k] = zk;
 
-  EZ = EX + EY;
+  /* Get the exponent sum into an intermediate variable.  This is a subtle
+     optimization, where given enough registers, all operations on the exponent
+     happen in registers and the result is written out only once into EZ.  */
+  int e = EX + EY;
+
   /* Is there a carry beyond the most significant digit?  */
   if (__glibc_unlikely (Z[1] == ZERO))
     {
-      for (i = 1; i <= p; i++)
+      for (i = 1; i <= p2; i++)
 	Z[i] = Z[i + 1];
-      EZ--;
+      e--;
     }
 
+  EZ = e;
   Z[0] = X[0] * Y[0];
 }
 
+/* Square *X and store result in *Y.  X and Y may not overlap.  For P in
+   [1, 2, 3], the exact result is truncated to P digits.  In case P > 3 the
+   error is bounded by 1.001 ULP.  This is a faster special case of
+   multiplication.  */
+void
+SECTION
+__sqr (const mp_no *x, mp_no *y, int p)
+{
+  long i, j, k, ip;
+  double u, yk;
+
+  /* Is z=0?  */
+  if (__glibc_unlikely (X[0] == ZERO))
+    {
+      Y[0] = ZERO;
+      return;
+    }
+
+  /* We need not iterate through all X's since it's pointless to
+     multiply zeroes.  */
+  for (ip = p; ip > 0; ip--)
+    if (X[ip] != ZERO)
+      break;
+
+  k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
+
+  while (k > 2 * ip + 1)
+    Y[k--] = ZERO;
+
+  yk = ZERO;
+
+  while (k > p)
+    {
+      double yk2 = 0.0;
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+	yk += X[lim] * X[lim];
+
+      /* In __mul, this loop (and the one within the next while loop) run
+         between a range to calculate the mantissa as follows:
+
+         Z[k] = X[k] * Y[n] + X[k+1] * Y[n-1] ... + X[n-1] * Y[k+1]
+		+ X[n] * Y[k]
+
+         For X == Y, we can get away with summing halfway and doubling the
+	 result.  For cases where the range size is even, the mid-point needs
+	 to be added separately (above).  */
+      for (i = k - p, j = p; i < j; i++, j--)
+	yk2 += X[i] * X[j];
+
+      yk += 2.0 * yk2;
+
+      u = (yk + CUTTER) - CUTTER;
+      if (u > yk)
+	u -= RADIX;
+      Y[k--] = yk - u;
+      yk = u * RADIXI;
+    }
+
+  while (k > 1)
+    {
+      double yk2 = 0.0;
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+	yk += X[lim] * X[lim];
+
+      /* Likewise for this loop.  */
+      for (i = 1, j = k - 1; i < j; i++, j--)
+	yk2 += X[i] * X[j];
+
+      yk += 2.0 * yk2;
+
+      u = (yk + CUTTER) - CUTTER;
+      if (u > yk)
+	u -= RADIX;
+      Y[k--] = yk - u;
+      yk = u * RADIXI;
+    }
+  Y[k] = yk;
+
+  /* Squares are always positive.  */
+  Y[0] = 1.0;
+
+  /* Get the exponent sum into an intermediate variable.  This is a subtle
+     optimization, where given enough registers, all operations on the exponent
+     happen in registers and the result is written out only once into EZ.  */
+  int e = EX * 2;
+
+  /* Is there a carry beyond the most significant digit?  */
+  if (__glibc_unlikely (Y[1] == ZERO))
+    {
+      for (i = 1; i <= p; i++)
+	Y[i] = Y[i + 1];
+      e--;
+    }
+
+  EY = e;
+}
+
 /* Invert *X and store in *Y.  Relative error bound:
    - For P = 2: 1.001 * R ^ (1 - P)
    - For P = 3: 1.063 * R ^ (1 - P)
@@ -664,7 +873,7 @@ static void
 SECTION
 __inv (const mp_no *x, mp_no *y, int p)
 {
-  int i;
+  long i;
   double t;
   mp_no z, w;
   static const int np1[] =
diff --git a/libc/sysdeps/ieee754/dbl-64/mpa.h b/libc/sysdeps/ieee754/dbl-64/mpa.h
index 06343d46d..168b334ed 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpa.h
+++ b/libc/sysdeps/ieee754/dbl-64/mpa.h
@@ -115,6 +115,7 @@ void __dbl_mp (double, mp_no *, int);
 void __add (const mp_no *, const mp_no *, mp_no *, int);
 void __sub (const mp_no *, const mp_no *, mp_no *, int);
 void __mul (const mp_no *, const mp_no *, mp_no *, int);
+void __sqr (const mp_no *, mp_no *, int);
 void __dvd (const mp_no *, const mp_no *, mp_no *, int);
 
 extern void __mpatan (mp_no *, mp_no *, int);
diff --git a/libc/sysdeps/ieee754/dbl-64/mpatan.c b/libc/sysdeps/ieee754/dbl-64/mpatan.c
index db5868092..cc879d8ec 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpatan.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpatan.c
@@ -39,63 +39,78 @@
 
 #include "mpatan.h"
 
-void __mpsqrt(mp_no *, mp_no *, int);
-
 void
 SECTION
-__mpatan(mp_no *x, mp_no *y, int p) {
+__mpatan (mp_no *x, mp_no *y, int p)
+{
 
-  int i,m,n;
+  int i, m, n;
   double dx;
-  mp_no
-    mptwoim1 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
-		0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
-		0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
+  mp_no mptwoim1 =
+  {
+    0,
+    {
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+      0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+    }
+  };
 
-  mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3;
+  mp_no mps, mpsm, mpt, mpt1, mpt2, mpt3;
 
-		      /* Choose m and initiate mptwoim1 */
-    if      (EX>0) m=7;
-    else if (EX<0) m=0;
-    else {
-      __mp_dbl(x,&dx,p);  dx=ABS(dx);
-      for (m=6; m>0; m--)
-	{if (dx>__atan_xm[m].d) break;}
+  /* Choose m and initiate mptwoim1.  */
+  if (EX > 0)
+    m = 7;
+  else if (EX < 0)
+    m = 0;
+  else
+    {
+      __mp_dbl (x, &dx, p);
+      dx = ABS (dx);
+      for (m = 6; m > 0; m--)
+	{
+	  if (dx > __atan_xm[m].d)
+	    break;
+	}
     }
-    mptwoim1.e = 1;
-    mptwoim1.d[0] = ONE;
+  mptwoim1.e = 1;
+  mptwoim1.d[0] = ONE;
 
-				 /* Reduce x m times */
-    __mul(x,x,&mpsm,p);
-    if (m==0) __cpy(x,&mps,p);
-    else {
-      for (i=0; i<m; i++) {
-	__add(&mpone,&mpsm,&mpt1,p);
-	__mpsqrt(&mpt1,&mpt2,p);
-	__add(&mpt2,&mpt2,&mpt1,p);
-	__add(&mptwo,&mpsm,&mpt2,p);
-	__add(&mpt1,&mpt2,&mpt3,p);
-	__dvd(&mpsm,&mpt3,&mpt1,p);
-	__cpy(&mpt1,&mpsm,p);
-      }
-      __mpsqrt(&mpsm,&mps,p);    mps.d[0] = X[0];
+  /* Reduce x m times.  */
+  __sqr (x, &mpsm, p);
+  if (m == 0)
+    __cpy (x, &mps, p);
+  else
+    {
+      for (i = 0; i < m; i++)
+	{
+	  __add (&mpone, &mpsm, &mpt1, p);
+	  __mpsqrt (&mpt1, &mpt2, p);
+	  __add (&mpt2, &mpt2, &mpt1, p);
+	  __add (&mptwo, &mpsm, &mpt2, p);
+	  __add (&mpt1, &mpt2, &mpt3, p);
+	  __dvd (&mpsm, &mpt3, &mpt1, p);
+	  __cpy (&mpt1, &mpsm, p);
+	}
+      __mpsqrt (&mpsm, &mps, p);
+      mps.d[0] = X[0];
     }
 
-		    /* Evaluate a truncated power series for Atan(s) */
-    n=__atan_np[p];    mptwoim1.d[1] = __atan_twonm1[p].d;
-    __dvd(&mpsm,&mptwoim1,&mpt,p);
-    for (i=n-1; i>1; i--) {
+  /* Evaluate a truncated power series for Atan(s).  */
+  n = __atan_np[p];
+  mptwoim1.d[1] = __atan_twonm1[p].d;
+  __dvd (&mpsm, &mptwoim1, &mpt, p);
+  for (i = n - 1; i > 1; i--)
+    {
       mptwoim1.d[1] -= TWO;
-      __dvd(&mpsm,&mptwoim1,&mpt1,p);
-      __mul(&mpsm,&mpt,&mpt2,p);
-      __sub(&mpt1,&mpt2,&mpt,p);
+      __dvd (&mpsm, &mptwoim1, &mpt1, p);
+      __mul (&mpsm, &mpt, &mpt2, p);
+      __sub (&mpt1, &mpt2, &mpt, p);
     }
-    __mul(&mps,&mpt,&mpt1,p);
-    __sub(&mps,&mpt1,&mpt,p);
-
-			  /* Compute Atan(x) */
-    mptwoim1.d[1] = 1 << m;
-    __mul(&mptwoim1,&mpt,y,p);
+  __mul (&mps, &mpt, &mpt1, p);
+  __sub (&mps, &mpt1, &mpt, p);
 
-  return;
+  /* Compute Atan(x).  */
+  mptwoim1.d[1] = 1 << m;
+  __mul (&mptwoim1, &mpt, y, p);
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/mpatan2.c b/libc/sysdeps/ieee754/dbl-64/mpatan2.c
index c0b9aea1e..d29c2fbad 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpatan2.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpatan2.c
@@ -32,37 +32,36 @@
 /*                                                                */
 /******************************************************************/
 
-
-
 #include "mpa.h"
 
 #ifndef SECTION
 # define SECTION
 #endif
 
-void __mpsqrt(mp_no *, mp_no *, int);
-void __mpatan(mp_no *, mp_no *, int);
-
-/* Multi-Precision Atan2(y,x) function subroutine, for p >= 4.    */
-/* y=0 is not permitted if x<=0. No error messages are given.     */
+/* Multi-Precision Atan2 (y, x) function subroutine, for p >= 4.
+   y = 0 is not permitted if x <= 0. No error messages are given.  */
 void
 SECTION
-__mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
+__mpatan2 (mp_no *y, mp_no *x, mp_no *z, int p)
+{
+  mp_no mpt1, mpt2, mpt3;
 
-  mp_no mpt1,mpt2,mpt3;
-
-
-  if (X[0] <= ZERO) {
-    __dvd(x,y,&mpt1,p);          __mul(&mpt1,&mpt1,&mpt2,p);
-    if (mpt1.d[0] != ZERO)       mpt1.d[0] = ONE;
-    __add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p);
-    __add(&mpt1,&mpt2,&mpt3,p);  mpt3.d[0]=Y[0];
-    __mpatan(&mpt3,&mpt1,p);     __add(&mpt1,&mpt1,z,p);
-  }
+  if (X[0] <= ZERO)
+    {
+      __dvd (x, y, &mpt1, p);
+      __mul (&mpt1, &mpt1, &mpt2, p);
+      if (mpt1.d[0] != ZERO)
+	mpt1.d[0] = ONE;
+      __add (&mpt2, &mpone, &mpt3, p);
+      __mpsqrt (&mpt3, &mpt2, p);
+      __add (&mpt1, &mpt2, &mpt3, p);
+      mpt3.d[0] = Y[0];
+      __mpatan (&mpt3, &mpt1, p);
+      __add (&mpt1, &mpt1, z, p);
+    }
   else
-  { __dvd(y,x,&mpt1,p);
-    __mpatan(&mpt1,z,p);
-  }
-
-  return;
+    {
+      __dvd (y, x, &mpt1, p);
+      __mpatan (&mpt1, z, p);
+    }
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/mpexp.c b/libc/sysdeps/ieee754/dbl-64/mpexp.c
index 8d288ff9a..565c6c853 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpexp.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpexp.c
@@ -49,6 +49,7 @@ __mpexp (mp_no *x, mp_no *y, int p)
       0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6,
       6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8
     };
+
   static const int m1p[33] =
     {
       0, 0, 0, 0,
@@ -71,16 +72,7 @@ __mpexp (mp_no *x, mp_no *y, int p)
       {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63},
       {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54}
     };
-  mp_no mpk =
-    {
-      0,
-      {
-	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
-      }
-    };
-  mp_no mps, mpak, mpt1, mpt2;
+  mp_no mps, mpk, mpt1, mpt2;
 
   /* Choose m,n and compute a=2**(-m).  */
   n = np[p];
@@ -115,37 +107,52 @@ __mpexp (mp_no *x, mp_no *y, int p)
 	  break;
     }
 
-  /* Compute s=x*2**(-m). Put result in mps.  */
+  /* Compute s=x*2**(-m). Put result in mps.  This is the range-reduced input
+     that we will use to compute e^s.  For the final result, simply raise it
+     to 2^m.  */
   __pow_mp (-m, &mpt1, p);
   __mul (x, &mpt1, &mps, p);
 
-  /* Evaluate the polynomial. Put result in mpt2.  */
-  mpk.e = 1;
-  mpk.d[0] = ONE;
-  mpk.d[1] = n;
-  __dvd (&mps, &mpk, &mpt1, p);
-  __add (&mpone, &mpt1, &mpak, p);
-  for (k = n - 1; k > 1; k--)
+  /* Compute the Taylor series for e^s:
+
+         1 + x/1! + x^2/2! + x^3/3! ...
+
+     for N iterations.  We compute this as:
+
+         e^x = 1 + (x * n!/1! + x^2 * n!/2! + x^3 * n!/3!) / n!
+             = 1 + (x * (n!/1! + x * (n!/2! + x * (n!/3! + x ...)))) / n!
+
+     k! is computed on the fly as KF and at the end of the polynomial loop, KF
+     is n!, which can be used directly.  */
+  __cpy (&mps, &mpt2, p);
+
+  double kf = 1.0;
+
+  /* Evaluate the rest.  The result will be in mpt2.  */
+  for (k = n - 1; k > 0; k--)
     {
-      __mul (&mps, &mpak, &mpt1, p);
-      mpk.d[1] = k;
-      __dvd (&mpt1, &mpk, &mpt2, p);
-      __add (&mpone, &mpt2, &mpak, p);
+      /* n! / k! = n * (n - 1) ... * (n - k + 1) */
+      kf *= k + 1;
+
+      __dbl_mp (kf, &mpk, p);
+      __add (&mpt2, &mpk, &mpt1, p);
+      __mul (&mps, &mpt1, &mpt2, p);
     }
-  __mul (&mps, &mpak, &mpt1, p);
+  __dbl_mp (kf, &mpk, p);
+  __dvd (&mpt2, &mpk, &mpt1, p);
   __add (&mpone, &mpt1, &mpt2, p);
 
   /* Raise polynomial value to the power of 2**m. Put result in y.  */
   for (k = 0, j = 0; k < m;)
     {
-      __mul (&mpt2, &mpt2, &mpt1, p);
+      __sqr (&mpt2, &mpt1, p);
       k++;
       if (k == m)
 	{
 	  j = 1;
 	  break;
 	}
-      __mul (&mpt1, &mpt1, &mpt2, p);
+      __sqr (&mpt1, &mpt2, p);
       k++;
     }
   if (j)
diff --git a/libc/sysdeps/ieee754/dbl-64/mplog.c b/libc/sysdeps/ieee754/dbl-64/mplog.c
index e3d10846e..f8d5c1095 100644
--- a/libc/sysdeps/ieee754/dbl-64/mplog.c
+++ b/libc/sysdeps/ieee754/dbl-64/mplog.c
@@ -1,4 +1,3 @@
-
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
@@ -37,27 +36,30 @@
 #include "endian.h"
 #include "mpa.h"
 
-void __mpexp(mp_no *, mp_no *, int);
-
-void __mplog(mp_no *x, mp_no *y, int p) {
-  int i,m;
-  static const int mp[33] = {0,0,0,0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
-                             4,4,4,4,4,4,4,4,4,4,4,4,4,4};
-  mp_no mpt1,mpt2;
+void
+__mplog (mp_no *x, mp_no *y, int p)
+{
+  int i, m;
+  static const int mp[33] =
+    {
+      0, 0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
+      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
+    };
+  mp_no mpt1, mpt2;
 
-  /* Choose m */
+  /* Choose m.  */
   m = mp[p];
 
-  /* Perform m newton iterations to solve for y: exp(y)-x=0.     */
-  /* The iterations formula is:  y(n+1)=y(n)+(x*exp(-y(n))-1).   */
-  __cpy(y,&mpt1,p);
-  for (i=0; i<m; i++) {
-    mpt1.d[0]=-mpt1.d[0];
-    __mpexp(&mpt1,&mpt2,p);
-    __mul(x,&mpt2,&mpt1,p);
-    __sub(&mpt1,&mpone,&mpt2,p);
-    __add(y,&mpt2,&mpt1,p);
-    __cpy(&mpt1,y,p);
-  }
-  return;
+  /* Perform m newton iterations to solve for y: exp(y) - x = 0.  The
+     iterations formula is:  y(n + 1) = y(n) + (x * exp(-y(n)) - 1).   */
+  __cpy (y, &mpt1, p);
+  for (i = 0; i < m; i++)
+    {
+      mpt1.d[0] = -mpt1.d[0];
+      __mpexp (&mpt1, &mpt2, p);
+      __mul (x, &mpt2, &mpt1, p);
+      __sub (&mpt1, &mpone, &mpt2, p);
+      __add (y, &mpt2, &mpt1, p);
+      __cpy (&mpt1, y, p);
+    }
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/mpsqrt.c b/libc/sysdeps/ieee754/dbl-64/mpsqrt.c
index 65df9fd06..71ef5ce77 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpsqrt.c
@@ -45,33 +45,37 @@
 /* p as integer. Routine computes sqrt(*x) and stores result in *y          */
 /****************************************************************************/
 
-static double fastiroot(double);
+static double fastiroot (double);
 
 void
 SECTION
-__mpsqrt(mp_no *x, mp_no *y, int p) {
-  int i,m,ey;
-  double dx,dy;
-  static const mp_no
-    mphalf   = {0,{1.0,8388608.0 /* 2^23 */}},
-    mp3halfs = {1,{1.0,1.0,8388608.0 /* 2^23 */}};
-  mp_no mpxn,mpz,mpu,mpt1,mpt2;
+__mpsqrt (mp_no *x, mp_no *y, int p)
+{
+  int i, m, ey;
+  double dx, dy;
+  static const mp_no mphalf = {0, {1.0, 8388608.0 /* 2^23 */}};
+  static const mp_no mp3halfs = {1, {1.0, 1.0, 8388608.0 /* 2^23 */}};
+  mp_no mpxn, mpz, mpu, mpt1, mpt2;
 
-  ey=EX/2;     __cpy(x,&mpxn,p);    mpxn.e -= (ey+ey);
-  __mp_dbl(&mpxn,&dx,p);   dy=fastiroot(dx);    __dbl_mp(dy,&mpu,p);
-  __mul(&mpxn,&mphalf,&mpz,p);
+  ey = EX / 2;
+  __cpy (x, &mpxn, p);
+  mpxn.e -= (ey + ey);
+  __mp_dbl (&mpxn, &dx, p);
+  dy = fastiroot (dx);
+  __dbl_mp (dy, &mpu, p);
+  __mul (&mpxn, &mphalf, &mpz, p);
 
-  m=__mpsqrt_mp[p];
-  for (i=0; i<m; i++) {
-    __mul(&mpu,&mpu,&mpt1,p);
-    __mul(&mpt1,&mpz,&mpt2,p);
-    __sub(&mp3halfs,&mpt2,&mpt1,p);
-    __mul(&mpu,&mpt1,&mpt2,p);
-    __cpy(&mpt2,&mpu,p);
-  }
-  __mul(&mpxn,&mpu,y,p);  EY += ey;
-
-  return;
+  m = __mpsqrt_mp[p];
+  for (i = 0; i < m; i++)
+    {
+      __sqr (&mpu, &mpt1, p);
+      __mul (&mpt1, &mpz, &mpt2, p);
+      __sub (&mp3halfs, &mpt2, &mpt1, p);
+      __mul (&mpu, &mpt1, &mpt2, p);
+      __cpy (&mpt2, &mpu, p);
+    }
+  __mul (&mpxn, &mpu, y, p);
+  EY += ey;
 }
 
 /***********************************************************/
@@ -80,22 +84,28 @@ __mpsqrt(mp_no *x, mp_no *y, int p) {
 /***********************************************************/
 static double
 SECTION
-fastiroot(double x) {
-  union {int i[2]; double d;} p,q;
-  double y,z, t;
+fastiroot (double x)
+{
+  union
+  {
+    int i[2];
+    double d;
+  } p, q;
+  double y, z, t;
   int n;
-  static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
+  static const double c0 = 0.99674, c1 = -0.53380;
+  static const double c2 = 0.45472, c3 = -0.21553;
 
   p.d = x;
-  p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
+  p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF) | 0x3FE00000;
   q.d = x;
   y = p.d;
-  z = y -1.0;
-  n = (q.i[HIGH_HALF] - p.i[HIGH_HALF])>>1;
-  z = ((c3*z + c2)*z + c1)*z + c0;            /* 2**-7         */
-  z = z*(1.5 - 0.5*y*z*z);                    /* 2**-14        */
-  p.d = z*(1.5 - 0.5*y*z*z);                  /* 2**-28        */
+  z = y - 1.0;
+  n = (q.i[HIGH_HALF] - p.i[HIGH_HALF]) >> 1;
+  z = ((c3 * z + c2) * z + c1) * z + c0;	/* 2**-7         */
+  z = z * (1.5 - 0.5 * y * z * z);		/* 2**-14        */
+  p.d = z * (1.5 - 0.5 * y * z * z);		/* 2**-28        */
   p.i[HIGH_HALF] -= n;
-  t = x*p.d;
-  return p.d*(1.5 - 0.5*p.d*t);
+  t = x * p.d;
+  return p.d * (1.5 - 0.5 * p.d * t);
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/mptan.c b/libc/sysdeps/ieee754/dbl-64/mptan.c
index 234108e37..51b5718e7 100644
--- a/libc/sysdeps/ieee754/dbl-64/mptan.c
+++ b/libc/sysdeps/ieee754/dbl-64/mptan.c
@@ -40,23 +40,25 @@
 # define SECTION
 #endif
 
-int __mpranred(double, mp_no *, int);
-void __c32(mp_no *, mp_no *, mp_no *, int);
-
 void
 SECTION
-__mptan(double x, mp_no *mpy, int p) {
+__mptan (double x, mp_no *mpy, int p)
+{
 
   int n;
   mp_no mpw, mpc, mps;
 
-  n = __mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */
-  __c32(&mpw, &mpc, &mps, p);              /* computing sin(x) and cos(x) */
-  if (n)                     /* second or fourth quarter of unit circle */
-  { __dvd(&mpc,&mps,mpy,p);
-    mpy->d[0] *= MONE;
-  }                          /* tan is negative in this area */
-  else  __dvd(&mps,&mpc,mpy,p);
-
-  return;
+  /* Negative or positive result.  */
+  n = __mpranred (x, &mpw, p) & 0x00000001;
+  /* Computing sin(x) and cos(x).  */
+  __c32 (&mpw, &mpc, &mps, p);
+  /* Second or fourth quarter of unit circle.  */
+  if (n)
+    {
+      __dvd (&mpc, &mps, mpy, p);
+      mpy->d[0] *= MONE;
+    }
+  /* tan is negative in this area.  */
+  else
+    __dvd (&mps, &mpc, mpy, p);
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/sincos32.c b/libc/sysdeps/ieee754/dbl-64/sincos32.c
index 6c5ffded5..3d2b2914a 100644
--- a/libc/sysdeps/ieee754/dbl-64/sincos32.c
+++ b/libc/sysdeps/ieee754/dbl-64/sincos32.c
@@ -57,17 +57,10 @@ SECTION
 ss32(mp_no *x, mp_no *y, int p) {
   int i;
   double a;
-#if 0
-  double b;
-  static const mp_no mpone = {1,{1.0,1.0}};
-#endif
   mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}};
-#if 0
-  mp_no mpt2;
-#endif
   for (i=1;i<=p;i++) mpk.d[i]=0;
 
-  __mul(x,x,&x2,p);
+  __sqr(x,&x2,p);
   __cpy(&oofac27,&gor,p);
   __cpy(&gor,&sum,p);
   for (a=27.0;a>1.0;a-=2.0) {
@@ -89,17 +82,10 @@ SECTION
 cc32(mp_no *x, mp_no *y, int p) {
   int i;
   double a;
-#if 0
-  double b;
-  static const mp_no mpone = {1,{1.0,1.0}};
-#endif
   mp_no mpt1,x2,gor,sum ,mpk={1,{1.0}};
-#if 0
-  mp_no mpt2;
-#endif
   for (i=1;i<=p;i++) mpk.d[i]=0;
 
-  __mul(x,x,&x2,p);
+  __sqr(x,&x2,p);
   mpk.d[1]=27.0;
   __mul(&oofac27,&mpk,&gor,p);
   __cpy(&gor,&sum,p);
@@ -119,7 +105,6 @@ cc32(mp_no *x, mp_no *y, int p) {
 void
 SECTION
 __c32(mp_no *x, mp_no *y, mp_no *z, int p) {
-  static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}};
   mp_no u,t,t1,t2,c,s;
   int i;
   __cpy(x,&u,p);
@@ -130,11 +115,11 @@ __c32(mp_no *x, mp_no *y, mp_no *z, int p) {
     __mul(&c,&s,&t,p);
     __sub(&s,&t,&t1,p);
     __add(&t1,&t1,&s,p);
-    __sub(&mpt,&c,&t1,p);
+    __sub(&mptwo,&c,&t1,p);
     __mul(&t1,&c,&t2,p);
     __add(&t2,&t2,&c,p);
   }
-  __sub(&one,&c,y,p);
+  __sub(&mpone,&c,y,p);
   __cpy(&s,z,p);
 }
 
@@ -251,7 +236,6 @@ __mpranred(double x, mp_no *y, int p)
   number v;
   double t,xn;
   int i,k,n;
-  static const mp_no one = {1,{1.0,1.0}};
   mp_no a,b,c;
 
   if (ABS(x) < 2.8e14) {
@@ -280,7 +264,7 @@ __mpranred(double x, mp_no *y, int p)
     c.e=0;
     if (c.d[1] >=  8388608.0)
     { t +=1.0;
-      __sub(&c,&one,&b,p);
+      __sub(&c,&mpone,&b,p);
       __mul(&b,&hp,y,p);
     }
     else __mul(&c,&hp,y,p);
diff --git a/libc/sysdeps/ieee754/dbl-64/slowexp.c b/libc/sysdeps/ieee754/dbl-64/slowexp.c
index 34ca3275e..c423fc311 100644
--- a/libc/sysdeps/ieee754/dbl-64/slowexp.c
+++ b/libc/sysdeps/ieee754/dbl-64/slowexp.c
@@ -34,38 +34,36 @@
 # define SECTION
 #endif
 
-void __mpexp(mp_no *x, mp_no *y, int p);
+void __mpexp (mp_no *x, mp_no *y, int p);
 
 /*Converting from double precision to Multi-precision and calculating  e^x */
 double
 SECTION
-__slowexp(double x) {
-  double w,z,res,eps=3.0e-26;
-#if 0
-  double y;
-#endif
+__slowexp (double x)
+{
+  double w, z, res, eps = 3.0e-26;
   int p;
-#if 0
-  int orig,i;
-#endif
-  mp_no mpx, mpy, mpz,mpw,mpeps,mpcor;
+  mp_no mpx, mpy, mpz, mpw, mpeps, mpcor;
 
-  p=6;
-  __dbl_mp(x,&mpx,p); /* Convert a double precision number  x               */
-		    /* into a multiple precision number mpx with prec. p. */
-  __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */
-  __dbl_mp(eps,&mpeps,p);
-  __mul(&mpeps,&mpy,&mpcor,p);
-  __add(&mpy,&mpcor,&mpw,p);
-  __sub(&mpy,&mpcor,&mpz,p);
-  __mp_dbl(&mpw, &w, p);
-  __mp_dbl(&mpz, &z, p);
-  if (w == z) return w;
-  else  {                   /* if calculating is not exactly   */
-    p = 32;
-    __dbl_mp(x,&mpx,p);
-    __mpexp(&mpx, &mpy, p);
-    __mp_dbl(&mpy, &res, p);
-    return res;
-  }
+  /* Use the multiple precision __MPEXP function to compute the exponential
+     First at 144 bits and if it is not accurate enough, at 768 bits.  */
+  p = 6;
+  __dbl_mp (x, &mpx, p);
+  __mpexp (&mpx, &mpy, p);
+  __dbl_mp (eps, &mpeps, p);
+  __mul (&mpeps, &mpy, &mpcor, p);
+  __add (&mpy, &mpcor, &mpw, p);
+  __sub (&mpy, &mpcor, &mpz, p);
+  __mp_dbl (&mpw, &w, p);
+  __mp_dbl (&mpz, &z, p);
+  if (w == z)
+    return w;
+  else
+    {
+      p = 32;
+      __dbl_mp (x, &mpx, p);
+      __mpexp (&mpx, &mpy, p);
+      __mp_dbl (&mpy, &res, p);
+      return res;
+    }
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/slowpow.c b/libc/sysdeps/ieee754/dbl-64/slowpow.c
index c303eaa5a..cccc7e32c 100644
--- a/libc/sysdeps/ieee754/dbl-64/slowpow.c
+++ b/libc/sysdeps/ieee754/dbl-64/slowpow.c
@@ -38,42 +38,59 @@
 # define SECTION
 #endif
 
-void __mpexp(mp_no *x, mp_no *y, int p);
-void __mplog(mp_no *x, mp_no *y, int p);
-double ulog(double);
-double __halfulp(double x,double y);
+void __mpexp (mp_no *x, mp_no *y, int p);
+void __mplog (mp_no *x, mp_no *y, int p);
+double ulog (double);
+double __halfulp (double x, double y);
 
 double
 SECTION
-__slowpow(double x, double y, double z) {
-  double res,res1;
-  mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1;
-  static const mp_no eps = {-3,{1.0,4.0}};
+__slowpow (double x, double y, double z)
+{
+  double res, res1;
+  mp_no mpx, mpy, mpz, mpw, mpp, mpr, mpr1;
+  static const mp_no eps = {-3, {1.0, 4.0}};
   int p;
 
-  res = __halfulp(x,y);        /* halfulp() returns -10 or x^y             */
-  if (res >= 0) return res;  /* if result was really computed by halfulp */
-		  /*  else, if result was not really computed by halfulp */
-  p = 10;         /*  p=precision   */
-  __dbl_mp(x,&mpx,p);
-  __dbl_mp(y,&mpy,p);
-  __dbl_mp(z,&mpz,p);
-  __mplog(&mpx, &mpz, p);     /* log(x) = z   */
-  __mul(&mpy,&mpz,&mpw,p);    /*  y * z =w    */
-  __mpexp(&mpw, &mpp, p);     /*  e^w =pp     */
-  __add(&mpp,&eps,&mpr,p);    /*  pp+eps =r   */
-  __mp_dbl(&mpr, &res, p);
-  __sub(&mpp,&eps,&mpr1,p);   /*  pp -eps =r1 */
-  __mp_dbl(&mpr1, &res1, p);  /*  converting into double precision */
-  if (res == res1) return res;
+  /* __HALFULP returns -10 or X^Y.  */
+  res = __halfulp (x, y);
 
-  p = 32;     /* if we get here result wasn't calculated exactly, continue */
-  __dbl_mp(x,&mpx,p);                          /* for more exact calculation */
-  __dbl_mp(y,&mpy,p);
-  __dbl_mp(z,&mpz,p);
-  __mplog(&mpx, &mpz, p);   /* log(c)=z  */
-  __mul(&mpy,&mpz,&mpw,p);  /* y*z =w    */
-  __mpexp(&mpw, &mpp, p);   /* e^w=pp    */
-  __mp_dbl(&mpp, &res, p);  /* converting into double precision */
+  /* Return if the result was computed by __HALFULP.  */
+  if (res >= 0)
+    return res;
+
+  /* Or else, calculate using multiple precision.  P = 10 implies accuracy of
+     240 bits accuracy, since MP_NO has a radix of 2^24.  */
+  p = 10;
+  __dbl_mp (x, &mpx, p);
+  __dbl_mp (y, &mpy, p);
+  __dbl_mp (z, &mpz, p);
+
+  /* z = x ^ y
+     log (z) = y * log (x)
+     z = exp (y * log (x))  */
+  __mplog (&mpx, &mpz, p);
+  __mul (&mpy, &mpz, &mpw, p);
+  __mpexp (&mpw, &mpp, p);
+
+  /* Add and subtract EPS to ensure that the result remains unchanged, i.e. we
+     have last bit accuracy.  */
+  __add (&mpp, &eps, &mpr, p);
+  __mp_dbl (&mpr, &res, p);
+  __sub (&mpp, &eps, &mpr1, p);
+  __mp_dbl (&mpr1, &res1, p);
+  if (res == res1)
+    return res;
+
+  /* If we don't, then we repeat using a higher precision.  768 bits of
+     precision ought to be enough for anybody.  */
+  p = 32;
+  __dbl_mp (x, &mpx, p);
+  __dbl_mp (y, &mpy, p);
+  __dbl_mp (z, &mpz, p);
+  __mplog (&mpx, &mpz, p);
+  __mul (&mpy, &mpz, &mpw, p);
+  __mpexp (&mpw, &mpp, p);
+  __mp_dbl (&mpp, &res, p);
   return res;
 }
diff --git a/libc/sysdeps/ieee754/dbl-64/x2y2m1.c b/libc/sysdeps/ieee754/dbl-64/x2y2m1.c
index 0b73f0a2e..d36a950e3 100644
--- a/libc/sysdeps/ieee754/dbl-64/x2y2m1.c
+++ b/libc/sysdeps/ieee754/dbl-64/x2y2m1.c
@@ -37,7 +37,7 @@ add_split (double *hi, double *lo, double x, double y)
    given that the values are small enough that no overflow occurs and
    large enough (or zero) that no underflow occurs.  */
 
-static inline void
+static void
 mul_split (double *hi, double *lo, double x, double y)
 {
 #ifdef __FP_FAST_FMA