Partial backports of:
=====================
 
commit c5d5d574cbfa96d0f6c1db24d1e072c472627e41
Author: Ondřej Bílka <neleai@seznam.cz>
Date:   Thu Oct 17 16:03:24 2013 +0200
 
    Format floating routines.
 
commit da08f647d58d674db08cdb3e61c8826c89470e2e
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Fri Dec 21 09:15:10 2012 +0530
 
    Move more constants into static variables
    
    Code cleanup.
 
commit f93a8d15699ee699282465dc1e03e033f3fabb52
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Wed Jan 16 16:06:48 2013 +0530
 
    Consolidate constant defines into mpa.h
 
commit caa99d06e7f1403887294442af520b0f8c6f3de0
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Fri Jan 18 11:18:13 2013 +0530
 
    Simplify calculation of 2^-m in __mpexp
 
commit 107a5bf085f5c4ef8c28266a34d476724cfc3475
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Tue Nov 18 15:40:56 2014 +0000
 
    Fix libm mpone, mptwo namespace (bug 17616).
 
To provided __mptwo for __inv.
 
Full backports of the following:
================================
 
commit 44e0d4c20ce5bf3825897e5d4b7caae94016214d
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Wed Jan 2 11:44:13 2013 +0530
 
    Split mantissa calculation loop and add branch prediction
 
commit f8af25d218202ff2f5d167b8e44e4b79f91d147f
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Fri Jan 4 15:09:33 2013 +0530
 
    Remove commented declarations
 
commit a9e48ab40e230c7fe34e4892bec8af4f3f975a20
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Fri Jan 4 15:42:09 2013 +0530
 
    Clean up comment for MP_NO
 
commit fffb407f4668b40b3df1eb8ee3ae3bc64ee79e20
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Fri Jan 4 22:52:12 2013 +0530
 
    Remove unused __cr and __cpymn
 
commit 950c99ca9094e7dc6394e90395f51e12093393aa
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Wed Jan 9 19:07:15 2013 +0530
 
    Update comments in mpa.c
    
    Fixed comment style and clearer wording in some cases.
 
commit 1066a53440d2744566e97c59bcd0d422186b3e90
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Mon Jan 14 21:31:25 2013 +0530
 
    Fix code formatting in mpa.c
    
    This includes the overridden mpa.c in power4.
 
commit 2a91b5735ac1bc65ce5c2a3646d75ba7208e26e9
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Mon Jan 14 21:36:58 2013 +0530
 
    Minor tweak to mp multiplication
    
    Add a local variable to remove extra copies to/from memory in the Z
    array.
 
ommit 45f058844c33f670475bd02f266942746bcb332b
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Tue Feb 26 21:28:16 2013 +0530
 
    Another tweak to the multiplication algorithm
    
    Reduce the formula to calculate mantissa so that we reduce the net
    number of multiplications performed.
 
commit bab8a695ee79a5a6e9b2b699938952b006fcbbec
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Thu Feb 21 14:29:18 2013 +0530
 
    Fix whitespace differences between generic and powerpc mpa.c
 
 
commit 2d0e0f29f85036d1189231cb7c1f19f27ba14a89
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Fri Feb 15 23:56:20 2013 +0530
 
    Fix determination of lower precision in __mul
 
commit 909279a5cfa938c989c9b01c8f48a0276291ec45
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Wed Feb 13 14:16:23 2013 +0530
 
    Optimized mp multiplication
    
    Don't bother multiplying zeroes since that only wastes cycles.
 
commit bdf028142eb77d6ae59500db875068fa5d7b059d
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Wed Feb 13 13:55:29 2013 +0530
 
    Clean up add_magnitudes and sub_magnitudes
 
commit d6752ccd696c71d23cd3df8fb9cc60b61c32e65a
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Thu Feb 14 10:31:09 2013 +0530
 
    New __sqr function as a faster special case of __mul
 
commit 4709fe7602b56e9f6ee1ab6afb4067409a784f29
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Sat Feb 16 00:09:29 2013 +0530
 
    Use intermediate variable to compute exponent in __mul
 
commit 20cd7fb3ae63795ae7c9a464abf5ed19b364ade0
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Wed Feb 20 18:56:20 2013 +0530
 
    Copy comment about inner loop from powerpc mpa.c to the default one
 
commit e69804d14e43f14c3c65dc570afdbfb822c9838b
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Mon Feb 25 16:43:02 2013 +0530
 
    Use long wherever possible in mpa.c
    
    Using long throughout like powerpc does is beneficial since it reduces
    the need to switch to 32-bit instructions.  It gives a very minor
    performance improvement.
 
commit 82a9811d29c00161c7c8ea7f70e2cc30988e192e
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Thu Mar 7 12:23:29 2013 +0530
 
    Use generic mpa.c code for everything except __mul and __sqr
 
commit 6f2e90e78f151bab153c2b38492505ae2012db06
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Tue Mar 26 19:28:50 2013 +0530
 
    Make mantissa type of mp_no configurable
 
    The mantissa of mp_no is intended to take only integral values.  This
    is a relatively good choice for powerpc due to its 4 fpus, but not for
    other architectures, which suffer due to this choice.  This change
    makes the default mantissa a long integer and allows powerpc to
    override it.  Additionally, some operations have been optimized for
    integer manipulation, resulting in a significant improvement in
    performance.
 
commit 5739f705eed5cf58e7b439e5983542e06d7fc2da
Author: Siddhesh Poyarekar <siddhesh@redhat.com>
Date:   Tue Mar 26 20:24:04 2013 +0530
 
    Use integral constants
    
    The compiler is smart enough to convert those into double for powerpc,
    but if we put them as doubles, it adds overhead by performing those
    operations in floating point mode.
 
commit 89f3b6e18c6e7833438789746fcfc2e7189f7cac
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Thu May 21 23:05:45 2015 +0000
 
    Fix sysdeps/ieee754/dbl-64/mpa.c for -Wuninitialized.
    
    If you remove the "override CFLAGS += -Wno-uninitialized" in
    math/Makefile, one of the errors you get is:
    
    ../sysdeps/ieee754/dbl-64/mpa.c: In function '__mp_dbl.part.0':
    ../sysdeps/ieee754/dbl-64/mpa.c:183:5: error: 'c' may be used uninitialized in this function [-Werror=maybe-uninitialized]
       c *= X[0];
    
    The problem is that the p < 5 case initializes c if p is 1, 2, 3 or 4
    but not otherwise, and in fact p is positive for all calls to this
    function so the uninitialized case can't actually occur.  This patch
    replaces the "if (p == 4)" last case with a comment so the compiler
    can see that all paths do initialize c.
 
    Tested for x86_64.
 
            * sysdeps/ieee754/dbl-64/mpa.c (norm): Remove if condition on
            (p == 4) case.
 
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpa.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.c
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation
+ * Copyright (C) 2001-2017 Free Software Foundation, Inc.
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -22,9 +22,7 @@
 /*  FUNCTIONS:                                                          */
 /*               mcr                                                    */
 /*               acr                                                    */
-/*               cr                                                     */
 /*               cpy                                                    */
-/*               cpymn                                                  */
 /*               norm                                                   */
 /*               denorm                                                 */
 /*               mp_dbl                                                 */
@@ -44,479 +42,868 @@
 
 #include "endian.h"
 #include "mpa.h"
-#include "mpa2.h"
-#include <sys/param.h>	/* For MIN() */
+#include <sys/param.h>
+#include <alloca.h>
 
 #ifndef SECTION
 # define SECTION
 #endif
 
+#ifndef NO__CONST
+/* TODO: With only a partial backport of the constant cleanup from
+	 upstream we limit the definition of these two constants to
+	 this file.  */
+static const mp_no __mpone = { 1, { 1.0, 1.0 } };
+static const mp_no __mptwo = { 1, { 1.0, 2.0 } };
+#endif
+
 #ifndef NO___ACR
-/* mcr() compares the sizes of the mantissas of two multiple precision  */
-/* numbers. Mantissas are compared regardless of the signs of the       */
-/* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also     */
-/* disregarded.                                                         */
+/* Compare mantissa of two multiple precision numbers regardless of the sign
+   and exponent of the numbers.  */
 static int
-mcr(const mp_no *x, const mp_no *y, int p) {
-  int i;
-  for (i=1; i<=p; i++) {
-    if      (X[i] == Y[i])  continue;
-    else if (X[i] >  Y[i])  return  1;
-    else                    return -1; }
+mcr (const mp_no *x, const mp_no *y, int p)
+{
+  long i;
+  long p2 = p;
+  for (i = 1; i <= p2; i++)
+    {
+      if (X[i] == Y[i])
+	continue;
+      else if (X[i] > Y[i])
+	return 1;
+      else
+	return -1;
+    }
   return 0;
 }
 
-
-/* acr() compares the absolute values of two multiple precision numbers */
+/* Compare the absolute values of two multiple precision numbers.  */
 int
-__acr(const mp_no *x, const mp_no *y, int p) {
-  int i;
-
-  if      (X[0] == ZERO) {
-    if    (Y[0] == ZERO) i= 0;
-    else                 i=-1;
-  }
-  else if (Y[0] == ZERO) i= 1;
-  else {
-    if      (EX >  EY)   i= 1;
-    else if (EX <  EY)   i=-1;
-    else                 i= mcr(x,y,p);
-  }
-
-  return i;
-}
-#endif
-
+__acr (const mp_no *x, const mp_no *y, int p)
+{
+  long i;
 
-#if 0
-/* cr() compares the values of two multiple precision numbers           */
-static int  __cr(const mp_no *x, const mp_no *y, int p) {
-  int i;
-
-  if      (X[0] > Y[0])  i= 1;
-  else if (X[0] < Y[0])  i=-1;
-  else if (X[0] < ZERO ) i= __acr(y,x,p);
-  else                   i= __acr(x,y,p);
+  if (X[0] == 0)
+    {
+      if (Y[0] == 0)
+	i = 0;
+      else
+	i = -1;
+    }
+  else if (Y[0] == 0)
+    i = 1;
+  else
+    {
+      if (EX > EY)
+	i = 1;
+      else if (EX < EY)
+	i = -1;
+      else
+	i = mcr (x, y, p);
+    }
 
   return i;
 }
 #endif
 
-
 #ifndef NO___CPY
-/* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
-void __cpy(const mp_no *x, mp_no *y, int p) {
+/* Copy multiple precision number X into Y.  They could be the same
+   number.  */
+void
+__cpy (const mp_no *x, mp_no *y, int p)
+{
+  long i;
+
   EY = EX;
-  for (int i=0; i <= p; i++)    Y[i] = X[i];
+  for (i = 0; i <= p; i++)
+    Y[i] = X[i];
 }
 #endif
 
+#ifndef NO___MP_DBL
+/* Convert a multiple precision number *X into a double precision
+   number *Y, normalized case (|x| >= 2**(-1022))).  X has precision
+   P, which is positive.  */
+static void
+norm (const mp_no *x, double *y, int p)
+{
+# define R RADIXI
+  long i;
+  double c;
+  mantissa_t a, u, v, z[5];
+  if (p < 5)
+    {
+      if (p == 1)
+	c = X[1];
+      else if (p == 2)
+	c = X[1] + R * X[2];
+      else if (p == 3)
+	c = X[1] + R * (X[2] + R * X[3]);
+      else /* p == 4.  */
+	c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
+    }
+  else
+    {
+      for (a = 1, z[1] = X[1]; z[1] < TWO23; )
+	{
+	  a *= 2;
+	  z[1] *= 2;
+	}
 
-#if 0
-/* Copy a multiple precision number x of precision m into a */
-/* multiple precision number y of precision n. In case n>m, */
-/* the digits of y beyond the m'th are set to zero. In case */
-/* n<m, the digits of x beyond the n'th are ignored.        */
-/* x=y is permissible.                                      */
-
-static void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
-
-  int i,k;
-
-  EY = EX;     k=MIN(m,n);
-  for (i=0; i <= k; i++)    Y[i] = X[i];
-  for (   ; i <= n; i++)    Y[i] = ZERO;
-}
-#endif
+      for (i = 2; i < 5; i++)
+	{
+	  mantissa_store_t d, r;
+	  d = X[i] * (mantissa_store_t) a;
+	  DIV_RADIX (d, r);
+	  z[i] = r;
+	  z[i - 1] += d;
+	}
 
+      u = ALIGN_DOWN_TO (z[3], TWO19);
+      v = z[3] - u;
 
-#ifndef NO___MP_DBL
-/* Convert a multiple precision number *x into a double precision */
-/* number *y, normalized case  (|x| >= 2**(-1022))) */
-static void norm(const mp_no *x, double *y, int p)
-{
-  #define R  radixi.d
-  int i;
-#if 0
-  int k;
-#endif
-  double a,c,u,v,z[5];
-  if (p<5) {
-    if      (p==1) c = X[1];
-    else if (p==2) c = X[1] + R* X[2];
-    else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
-    else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
-  }
-  else {
-    for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
-	{a *= TWO;   z[1] *= TWO; }
-
-    for (i=2; i<5; i++) {
-      z[i] = X[i]*a;
-      u = (z[i] + CUTTER)-CUTTER;
-      if  (u > z[i])  u -= RADIX;
-      z[i] -= u;
-      z[i-1] += u*RADIXI;
-    }
-
-    u = (z[3] + TWO71) - TWO71;
-    if (u > z[3])   u -= TWO19;
-    v = z[3]-u;
-
-    if (v == TWO18) {
-      if (z[4] == ZERO) {
-	for (i=5; i <= p; i++) {
-	  if (X[i] == ZERO)   continue;
-	  else                {z[3] += ONE;   break; }
+      if (v == TWO18)
+	{
+	  if (z[4] == 0)
+	    {
+	      for (i = 5; i <= p; i++)
+		{
+		  if (X[i] == 0)
+		    continue;
+		  else
+		    {
+		      z[3] += 1;
+		      break;
+		    }
+		}
+	    }
+	  else
+	    z[3] += 1;
 	}
-      }
-      else              z[3] += ONE;
-    }
 
-    c = (z[1] + R *(z[2] + R * z[3]))/a;
-  }
+      c = (z[1] + R * (z[2] + R * z[3])) / a;
+    }
 
   c *= X[0];
 
-  for (i=1; i<EX; i++)   c *= RADIX;
-  for (i=1; i>EX; i--)   c *= RADIXI;
+  for (i = 1; i < EX; i++)
+    c *= RADIX;
+  for (i = 1; i > EX; i--)
+    c *= RADIXI;
 
   *y = c;
-#undef R
+# undef R
 }
 
-/* Convert a multiple precision number *x into a double precision */
-/* number *y, denormalized case  (|x| < 2**(-1022))) */
-static void denorm(const mp_no *x, double *y, int p)
+/* Convert a multiple precision number *X into a double precision
+   number *Y, Denormal case  (|x| < 2**(-1022))).  */
+static void
+denorm (const mp_no *x, double *y, int p)
 {
-  int i,k;
-  double c,u,z[5];
-#if 0
-  double a,v;
-#endif
+  long i, k;
+  long p2 = p;
+  double c;
+  mantissa_t u, z[5];
+
+# define R RADIXI
+  if (EX < -44 || (EX == -44 && X[1] < TWO5))
+    {
+      *y = 0;
+      return;
+    }
 
-#define R  radixi.d
-  if (EX<-44 || (EX==-44 && X[1]<TWO5))
-     { *y=ZERO; return; }
-
-  if      (p==1) {
-    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
-    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
-    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
-  }
-  else if (p==2) {
-    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
-    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
-    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
-  }
-  else {
-    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
-    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
-    else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
-    z[3] = X[k];
-  }
-
-  u = (z[3] + TWO57) - TWO57;
-  if  (u > z[3])   u -= TWO5;
-
-  if (u==z[3]) {
-    for (i=k+1; i <= p; i++) {
-      if (X[i] == ZERO)   continue;
-      else {z[3] += ONE;   break; }
-    }
-  }
-
-  c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
-
-  *y = c*TWOM1032;
-#undef R
-}
-
-/* Convert a multiple precision number *x into a double precision number *y. */
-/* The result is correctly rounded to the nearest/even. *x is left unchanged */
-
-void __mp_dbl(const mp_no *x, double *y, int p) {
-#if 0
-  int i,k;
-  double a,c,u,v,z[5];
-#endif
+  if (p2 == 1)
+    {
+      if (EX == -42)
+	{
+	  z[1] = X[1] + TWO10;
+	  z[2] = 0;
+	  z[3] = 0;
+	  k = 3;
+	}
+      else if (EX == -43)
+	{
+	  z[1] = TWO10;
+	  z[2] = X[1];
+	  z[3] = 0;
+	  k = 2;
+	}
+      else
+	{
+	  z[1] = TWO10;
+	  z[2] = 0;
+	  z[3] = X[1];
+	  k = 1;
+	}
+    }
+  else if (p2 == 2)
+    {
+      if (EX == -42)
+	{
+	  z[1] = X[1] + TWO10;
+	  z[2] = X[2];
+	  z[3] = 0;
+	  k = 3;
+	}
+      else if (EX == -43)
+	{
+	  z[1] = TWO10;
+	  z[2] = X[1];
+	  z[3] = X[2];
+	  k = 2;
+	}
+      else
+	{
+	  z[1] = TWO10;
+	  z[2] = 0;
+	  z[3] = X[1];
+	  k = 1;
+	}
+    }
+  else
+    {
+      if (EX == -42)
+	{
+	  z[1] = X[1] + TWO10;
+	  z[2] = X[2];
+	  k = 3;
+	}
+      else if (EX == -43)
+	{
+	  z[1] = TWO10;
+	  z[2] = X[1];
+	  k = 2;
+	}
+      else
+	{
+	  z[1] = TWO10;
+	  z[2] = 0;
+	  k = 1;
+	}
+      z[3] = X[k];
+    }
 
-  if (X[0] == ZERO)  {*y = ZERO;  return; }
+  u = ALIGN_DOWN_TO (z[3], TWO5);
 
-  if      (EX> -42)                 norm(x,y,p);
-  else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
-  else                              denorm(x,y,p);
+  if (u == z[3])
+    {
+      for (i = k + 1; i <= p2; i++)
+	{
+	  if (X[i] == 0)
+	    continue;
+	  else
+	    {
+	      z[3] += 1;
+	      break;
+	    }
+	}
+    }
+
+  c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
+
+  *y = c * TWOM1032;
+# undef R
 }
-#endif
 
+/* Convert multiple precision number *X into double precision number *Y.  The
+   result is correctly rounded to the nearest/even.  */
+void
+__mp_dbl (const mp_no *x, double *y, int p)
+{
+  if (X[0] == 0)
+    {
+      *y = 0;
+      return;
+    }
 
-/* dbl_mp() converts a double precision number x into a multiple precision  */
-/* number *y. If the precision p is too small the result is truncated. x is */
-/* left unchanged.                                                          */
+  if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
+    norm (x, y, p);
+  else
+    denorm (x, y, p);
+}
+#endif
 
+/* Get the multiple precision equivalent of X into *Y.  If the precision is too
+   small, the result is truncated.  */
 void
 SECTION
-__dbl_mp(double x, mp_no *y, int p) {
+__dbl_mp (double x, mp_no *y, int p)
+{
+  long i, n;
+  long p2 = p;
 
-  int i,n;
-  double u;
+  /* Sign.  */
+  if (x == 0)
+    {
+      Y[0] = 0;
+      return;
+    }
+  else if (x > 0)
+    Y[0] = 1;
+  else
+    {
+      Y[0] = -1;
+      x = -x;
+    }
 
-  /* Sign */
-  if      (x == ZERO)  {Y[0] = ZERO;  return; }
-  else if (x >  ZERO)   Y[0] = ONE;
-  else                 {Y[0] = MONE;  x=-x;   }
-
-  /* Exponent */
-  for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
-  for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;
-
-  /* Digits */
-  n=MIN(p,4);
-  for (i=1; i<=n; i++) {
-    u = (x + TWO52) - TWO52;
-    if (u>x)   u -= ONE;
-    Y[i] = u;     x -= u;    x *= RADIX; }
-  for (   ; i<=p; i++)     Y[i] = ZERO;
+  /* Exponent.  */
+  for (EY = 1; x >= RADIX; EY += 1)
+    x *= RADIXI;
+  for (; x < 1; EY -= 1)
+    x *= RADIX;
+
+  /* Digits.  */
+  n = MIN (p2, 4);
+  for (i = 1; i <= n; i++)
+    {
+      INTEGER_OF (x, Y[i]);
+      x *= RADIX;
+    }
+  for (; i <= p2; i++)
+    Y[i] = 0;
 }
 
-
-/*  add_magnitudes() adds the magnitudes of *x & *y assuming that           */
-/*  abs(*x) >= abs(*y) > 0.                                                 */
-/* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
-/* No guard digit is used. The result equals the exact sum, truncated.      */
-/* *x & *y are left unchanged.                                              */
-
+/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0.  The
+   sign of the sum *Z is not changed.  X and Y may overlap but not X and Z or
+   Y and Z.  No guard digit is used.  The result equals the exact sum,
+   truncated.  */
 static void
 SECTION
-add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
-  int i,j,k;
+add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+  long i, j, k;
+  long p2 = p;
+  mantissa_t zk;
 
   EZ = EX;
 
-  i=p;    j=p+ EY - EX;    k=p+1;
+  i = p2;
+  j = p2 + EY - EX;
+  k = p2 + 1;
+
+  if (__glibc_unlikely (j < 1))
+    {
+      __cpy (x, z, p);
+      return;
+    }
 
-  if (j<1)
-     {__cpy(x,z,p);  return; }
-  else   Z[k] = ZERO;
-
-  for (; j>0; i--,j--) {
-    Z[k] += X[i] + Y[j];
-    if (Z[k] >= RADIX) {
-      Z[k]  -= RADIX;
-      Z[--k] = ONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  for (; i>0; i--) {
-    Z[k] += X[i];
-    if (Z[k] >= RADIX) {
-      Z[k]  -= RADIX;
-      Z[--k] = ONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  if (Z[1] == ZERO) {
-    for (i=1; i<=p; i++)    Z[i] = Z[i+1]; }
-  else   EZ += ONE;
-}
+  zk = 0;
 
+  for (; j > 0; i--, j--)
+    {
+      zk += X[i] + Y[j];
+      if (zk >= RADIX)
+	{
+	  Z[k--] = zk - RADIX;
+	  zk = 1;
+	}
+      else
+	{
+	  Z[k--] = zk;
+	  zk = 0;
+	}
+    }
 
-/*  sub_magnitudes() subtracts the magnitudes of *x & *y assuming that      */
-/*  abs(*x) > abs(*y) > 0.                                                  */
-/* The sign of the difference *z is undefined. x&y may overlap but not x&z  */
-/* or y&z. One guard digit is used. The error is less than one ulp.         */
-/* *x & *y are left unchanged.                                              */
+  for (; i > 0; i--)
+    {
+      zk += X[i];
+      if (zk >= RADIX)
+	{
+	  Z[k--] = zk - RADIX;
+	  zk = 1;
+	}
+      else
+	{
+	  Z[k--] = zk;
+	  zk = 0;
+	}
+    }
 
+  if (zk == 0)
+    {
+      for (i = 1; i <= p2; i++)
+	Z[i] = Z[i + 1];
+    }
+  else
+    {
+      Z[1] = zk;
+      EZ += 1;
+    }
+}
+
+/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
+   The sign of the difference *Z is not changed.  X and Y may overlap but not X
+   and Z or Y and Z.  One guard digit is used.  The error is less than one
+   ULP.  */
 static void
 SECTION
-sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
-  int i,j,k;
+sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+  long i, j, k;
+  long p2 = p;
+  mantissa_t zk;
 
   EZ = EX;
+  i = p2;
+  j = p2 + EY - EX;
+  k = p2;
+
+  /* Y is too small compared to X, copy X over to the result.  */
+  if (__glibc_unlikely (j < 1))
+    {
+      __cpy (x, z, p);
+      return;
+    }
+
+  /* The relevant least significant digit in Y is non-zero, so we factor it in
+     to enhance accuracy.  */
+  if (j < p2 && Y[j + 1] > 0)
+    {
+      Z[k + 1] = RADIX - Y[j + 1];
+      zk = -1;
+    }
+  else
+    zk = Z[k + 1] = 0;
 
-  if (EX == EY) {
-    i=j=k=p;
-    Z[k] = Z[k+1] = ZERO; }
-  else {
-    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--;}
-    }
-  }
-
-  for (; j>0; i--,j--) {
-    Z[k] += (X[i] - Y[j]);
-    if (Z[k] < ZERO) {
-      Z[k]  += RADIX;
-      Z[--k] = MONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  for (; i>0; i--) {
-    Z[k] += X[i];
-    if (Z[k] < ZERO) {
-      Z[k]  += RADIX;
-      Z[--k] = MONE; }
-    else
-      Z[--k] = ZERO;
-  }
+  /* Subtract and borrow.  */
+  for (; j > 0; i--, j--)
+    {
+      zk += (X[i] - Y[j]);
+      if (zk < 0)
+	{
+	  Z[k--] = zk + RADIX;
+	  zk = -1;
+	}
+      else
+	{
+	  Z[k--] = zk;
+	  zk = 0;
+	}
+    }
 
-  for (i=1; Z[i] == ZERO; i++) ;
+  /* We're done with digits from Y, so it's just digits in X.  */
+  for (; i > 0; i--)
+    {
+      zk += X[i];
+      if (zk < 0)
+	{
+	  Z[k--] = zk + RADIX;
+	  zk = -1;
+	}
+      else
+	{
+	  Z[k--] = zk;
+	  zk = 0;
+	}
+    }
+
+  /* Normalize.  */
+  for (i = 1; Z[i] == 0; i++)
+    ;
   EZ = EZ - i + 1;
-  for (k=1; i <= p+1; )
+  for (k = 1; i <= p2 + 1; )
     Z[k++] = Z[i++];
-  for (; k <= p; )
-    Z[k++] = ZERO;
+  for (; k <= p2; )
+    Z[k++] = 0;
 }
 
-
-/* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap  */
-/* but not x&z or y&z. One guard digit is used. The error is less than    */
-/* one ulp. *x & *y are left unchanged.                                   */
-
+/* Add *X and *Y and store the result in *Z.  X and Y may overlap, but not X
+   and Z or Y and Z.  One guard digit is used.  The error is less than one
+   ULP.  */
 void
 SECTION
-__add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
+__add (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
   int n;
 
-  if      (X[0] == ZERO)     {__cpy(y,z,p);  return; }
-  else if (Y[0] == ZERO)     {__cpy(x,z,p);  return; }
+  if (X[0] == 0)
+    {
+      __cpy (y, z, p);
+      return;
+    }
+  else if (Y[0] == 0)
+    {
+      __cpy (x, z, p);
+      return;
+    }
 
-  if (X[0] == Y[0])   {
-    if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] = X[0]; }
-    else                     {add_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
-  }
-  else                       {
-    if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] = X[0]; }
-    else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
-    else                      Z[0] = ZERO;
-  }
+  if (X[0] == Y[0])
+    {
+      if (__acr (x, y, p) > 0)
+	{
+	  add_magnitudes (x, y, z, p);
+	  Z[0] = X[0];
+	}
+      else
+	{
+	  add_magnitudes (y, x, z, p);
+	  Z[0] = Y[0];
+	}
+    }
+  else
+    {
+      if ((n = __acr (x, y, p)) == 1)
+	{
+	  sub_magnitudes (x, y, z, p);
+	  Z[0] = X[0];
+	}
+      else if (n == -1)
+	{
+	  sub_magnitudes (y, x, z, p);
+	  Z[0] = Y[0];
+	}
+      else
+	Z[0] = 0;
+    }
 }
 
+/* Subtract *Y from *X and return the result in *Z.  X and Y may overlap but
+   not X and Z or Y and Z.  One guard digit is used.  The error is less than
+   one ULP.  */
+void
+SECTION
+__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+  int n;
+
+  if (X[0] == 0)
+    {
+      __cpy (y, z, p);
+      Z[0] = -Z[0];
+      return;
+    }
+  else if (Y[0] == 0)
+    {
+      __cpy (x, z, p);
+      return;
+    }
 
-/* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
-/* overlap but not x&z or y&z. One guard digit is used. The error is      */
-/* less than one ulp. *x & *y are left unchanged.                         */
+  if (X[0] != Y[0])
+    {
+      if (__acr (x, y, p) > 0)
+	{
+	  add_magnitudes (x, y, z, p);
+	  Z[0] = X[0];
+	}
+      else
+	{
+	  add_magnitudes (y, x, z, p);
+	  Z[0] = -Y[0];
+	}
+    }
+  else
+    {
+      if ((n = __acr (x, y, p)) == 1)
+	{
+	  sub_magnitudes (x, y, z, p);
+	  Z[0] = X[0];
+	}
+      else if (n == -1)
+	{
+	  sub_magnitudes (y, x, z, p);
+	  Z[0] = -Y[0];
+	}
+      else
+	Z[0] = 0;
+    }
+}
 
+#ifndef NO__MUL
+/* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
+   and Z or Y and Z.  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.  */
 void
 SECTION
-__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+__mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
+  long i, j, k, ip, ip2;
+  long p2 = p;
+  mantissa_store_t zk;
+  const mp_no *a;
+  mantissa_store_t *diag;
+
+  /* Is z=0?  */
+  if (__glibc_unlikely (X[0] * Y[0] == 0))
+    {
+      Z[0] = 0;
+      return;
+    }
 
-  int n;
+  /* 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] != 0 || Y[ip2] != 0)
+      break;
+
+  a = X[ip2] != 0 ? y : x;
+
+  /* ... and here, at least one of them is still zero.  */
+  for (ip = ip2; ip > 0; ip--)
+    if (a->d[ip] != 0)
+      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 0.  */
+
+  k = (__glibc_unlikely (p2 < 3)) ? p2 + p2 : p2 + 3;
+
+  while (k > ip + ip2 + 1)
+    Z[k--] = 0;
+
+  zk = 0;
+
+  /* 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 (mantissa_store_t));
+  mantissa_store_t d = 0;
+  for (i = 1; i <= ip; i++)
+    {
+      d += X[i] * (mantissa_store_t) Y[i];
+      diag[i] = d;
+    }
+  while (i < k)
+    diag[i++] = d;
 
-  if      (X[0] == ZERO)     {__cpy(y,z,p);  Z[0] = -Z[0];  return; }
-  else if (Y[0] == ZERO)     {__cpy(x,z,p);                 return; }
+  while (k > p2)
+    {
+      long lim = k / 2;
 
-  if (X[0] != Y[0])    {
-    if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
-    else                     {add_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
-  }
-  else                       {
-    if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
-    else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
-    else                      Z[0] = ZERO;
-  }
-}
+      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] * (mantissa_store_t) Y[lim];
 
+      for (i = k - p2, j = p2; i < j; i++, j--)
+	zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
 
-/* Multiply two multiple precision numbers. *z is set to *x * *y. x&y      */
-/* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is     */
-/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp.   */
-/* *x & *y are left unchanged.                                             */
+      zk -= diag[k - 1];
 
-void
-SECTION
-__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+      DIV_RADIX (zk, Z[k]);
+      k--;
+    }
 
-  int i, i1, i2, j, k, k2;
-  double u;
+  /* 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)
+    {
+      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] * (mantissa_store_t) Y[lim];
 
-		      /* Is z=0? */
-  if (X[0]*Y[0]==ZERO)
-     { Z[0]=ZERO;  return; }
-
-		       /* Multiply, add and carry */
-  k2 = (p<3) ? p+p : p+3;
-  Z[k2]=ZERO;
-  for (k=k2; k>1; ) {
-    if (k > p)  {i1=k-p; i2=p+1; }
-    else        {i1=1;   i2=k;   }
-    for (i=i1,j=i2-1; i<i2; i++,j--)  Z[k] += X[i]*Y[j];
-
-    u = (Z[k] + CUTTER)-CUTTER;
-    if  (u > Z[k])  u -= RADIX;
-    Z[k]  -= u;
-    Z[--k] = u*RADIXI;
-  }
-
-		 /* Is there a carry beyond the most significant digit? */
-  if (Z[1] == ZERO) {
-    for (i=1; i<=p; i++)  Z[i]=Z[i+1];
-    EZ = EX + EY - 1; }
-  else
-    EZ = EX + EY;
+      for (i = 1, j = k - 1; i < j; i++, j--)
+	zk += (X[i] + X[j]) * (mantissa_store_t) (Y[i] + Y[j]);
+
+      zk -= diag[k - 1];
+
+      DIV_RADIX (zk, Z[k]);
+      k--;
+    }
+  Z[k] = zk;
+
+  /* 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] == 0))
+    {
+      for (i = 1; i <= p2; i++)
+	Z[i] = Z[i + 1];
+      e--;
+    }
 
+  EZ = e;
   Z[0] = X[0] * Y[0];
 }
+#endif
+
+#ifndef NO__SQR
+/* 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;
+  mantissa_store_t yk;
 
+  /* Is z=0?  */
+  if (__glibc_unlikely (X[0] == 0))
+    {
+      Y[0] = 0;
+      return;
+    }
 
-/* Invert a multiple precision number. Set *y = 1 / *x.                     */
-/* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3,   */
-/* 2.001*r**(1-p) for p>3.                                                  */
-/* *x=0 is not permissible. *x is left unchanged.                           */
+  /* We need not iterate through all X's since it's pointless to
+     multiply zeroes.  */
+  for (ip = p; ip > 0; ip--)
+    if (X[ip] != 0)
+      break;
 
-static
-SECTION
-void __inv(const mp_no *x, mp_no *y, int p) {
-  int i;
-#if 0
-  int l;
+  k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
+
+  while (k > 2 * ip + 1)
+    Y[k--] = 0;
+
+  yk = 0;
+
+  while (k > p)
+    {
+      mantissa_store_t yk2 = 0;
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+	yk += X[lim] * (mantissa_store_t) 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] * (mantissa_store_t) X[j];
+
+      yk += 2 * yk2;
+
+      DIV_RADIX (yk, Y[k]);
+      k--;
+    }
+
+  while (k > 1)
+    {
+      mantissa_store_t yk2 = 0;
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+	yk += X[lim] * (mantissa_store_t) X[lim];
+
+      /* Likewise for this loop.  */
+      for (i = 1, j = k - 1; i < j; i++, j--)
+	yk2 += X[i] * (mantissa_store_t) X[j];
+
+      yk += 2 * yk2;
+
+      DIV_RADIX (yk, Y[k]);
+      k--;
+    }
+  Y[k] = yk;
+
+  /* Squares are always positive.  */
+  Y[0] = 1;
+
+  /* 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] == 0))
+    {
+      for (i = 1; i <= p; i++)
+	Y[i] = Y[i + 1];
+      e--;
+    }
+
+  EY = e;
+}
 #endif
+
+/* 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)
+   - For P > 3: 2.001 * R ^ (1 - P)
+
+   *X = 0 is not permissible.  */
+static void
+SECTION
+__inv (const mp_no *x, mp_no *y, int p)
+{
+  long i;
   double t;
-  mp_no z,w;
-  static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
-			    4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
-  const mp_no mptwo = {1,{1.0,2.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}};
-
-  __cpy(x,&z,p);  z.e=0;  __mp_dbl(&z,&t,p);
-  t=ONE/t;   __dbl_mp(t,y,p);    EY -= EX;
-
-  for (i=0; i<np1[p]; i++) {
-    __cpy(y,&w,p);
-    __mul(x,&w,y,p);
-    __sub(&mptwo,y,&z,p);
-    __mul(&w,&z,y,p);
-  }
+  mp_no z, w;
+  static const int np1[] =
+    { 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
+    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
+  };
+
+  __cpy (x, &z, p);
+  z.e = 0;
+  __mp_dbl (&z, &t, p);
+  t = 1 / t;
+  __dbl_mp (t, y, p);
+  EY -= EX;
+
+  for (i = 0; i < np1[p]; i++)
+    {
+      __cpy (y, &w, p);
+      __mul (x, &w, y, p);
+      __sub (&__mptwo, y, &z, p);
+      __mul (&w, &z, y, p);
+    }
 }
 
+/* Divide *X by *Y and store result in *Z.  X and Y may overlap but not X and Z
+   or Y and Z.  Relative error bound:
+   - For P = 2: 2.001 * R ^ (1 - P)
+   - For P = 3: 2.063 * R ^ (1 - P)
+   - For P > 3: 3.001 * R ^ (1 - P)
 
-/* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
-/* are left unchanged. x&y may overlap but not x&z or y&z.                   */
-/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3     */
-/* and 3.001*r**(1-p) for p>3. *y=0 is not permissible.                      */
-
+   *X = 0 is not permissible.  */
 void
 SECTION
-__dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
+__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p)
+{
   mp_no w;
 
-  if (X[0] == ZERO)    Z[0] = ZERO;
-  else                {__inv(y,&w,p);   __mul(x,&w,z,p);}
+  if (X[0] == 0)
+    Z[0] = 0;
+  else
+    {
+      __inv (y, &w, p);
+      __mul (x, &w, z, p);
+    }
 }
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpa.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa.h
@@ -1,7 +1,7 @@
 /*
  * IBM Accurate Mathematical Library
  * Written by International Business Machines Corp.
- * Copyright (C) 2001, 2011 Free Software Foundation, Inc.
+ * Copyright (C) 2001-2017 Free Software Foundation, Inc.
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -23,36 +23,58 @@
 /*  FUNCTIONS:                                                          */
 /*               mcr                                                    */
 /*               acr                                                    */
-/*               cr                                                     */
 /*               cpy                                                    */
-/*               cpymn                                                  */
 /*               mp_dbl                                                 */
 /*               dbl_mp                                                 */
 /*               add                                                    */
 /*               sub                                                    */
 /*               mul                                                    */
-/*               inv                                                    */
 /*               dvd                                                    */
 /*                                                                      */
 /* Arithmetic functions for multiple precision numbers.                 */
 /* Common types and definition                                          */
 /************************************************************************/
 
+#include <mpa-arch.h>
 
-typedef struct {/* This structure holds the details of a multi-precision     */
-  int e;        /* floating point number, x: d[0] holds its sign (-1,0 or 1) */
-  double d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and            */
-} mp_no;        /* d[1]...d[p] hold its mantissa digits. The value of x is,  */
-		/* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p).  */
-		/* Here   r = 2**24,   0 <= d[i] < r  and  1 <= p <= 32.     */
-		/* p is a global variable. A multi-precision number is       */
-		/* always normalized. Namely, d[1] > 0. An exception is      */
-		/* a zero which is characterized by d[0] = 0. The terms      */
-		/* d[p+1], d[p+2], ... of a none zero number have no         */
-		/* significance and so are the terms e, d[1],d[2],...        */
-		/* of a zero.                                                */
+/* The mp_no structure holds the details of a multi-precision floating point
+   number.
 
-typedef union { int i[2]; double d; } number;
+   - The radix of the number (R) is 2 ^ 24.
+
+   - E: The exponent of the number.
+
+   - D[0]: The sign (-1, 1) or 0 if the value is 0.  In the latter case, the
+     values of the remaining members of the structure are ignored.
+
+   - D[1] - D[p]: The mantissa of the number where:
+
+	0 <= D[i] < R and
+	P is the precision of the number and 1 <= p <= 32
+
+     D[p+1] ... D[39] have no significance.
+
+   - The value of the number is:
+
+	D[1] * R ^ (E - 1) + D[2] * R ^ (E - 2) ... D[p] * R ^ (E - p)
+
+   */
+typedef struct
+{
+  int e;
+  mantissa_t d[40];
+} mp_no;
+
+typedef union
+{
+  int i[2];
+  double d;
+} number;
+
+/* TODO: With only a partial backport of the constant cleanup we don't
+         define __mpone or __mptwo here for other code to use.  */
+/* extern const mp_no __mpone;
+extern const mp_no __mptwo;  */
 
 #define  X   x->d
 #define  Y   y->d
@@ -63,21 +85,73 @@ typedef union { int i[2]; double d; } nu
 
 #define ABS(x)   ((x) <  0  ? -(x) : (x))
 
-int __acr(const mp_no *, const mp_no *, int);
-// int  __cr(const mp_no *, const mp_no *, int);
-void __cpy(const mp_no *, mp_no *, int);
-// void __cpymn(const mp_no *, int, mp_no *, int);
-void __mp_dbl(const mp_no *, double *, int);
-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 __inv(const mp_no *, mp_no *, int);
-void __dvd(const mp_no *, const mp_no *, mp_no *, int);
+#ifndef RADIXI
+# define  RADIXI    0x1.0p-24		/* 2^-24   */
+#endif
+
+#ifndef TWO52
+# define  TWO52     0x1.0p52		/* 2^52    */
+#endif
+
+#define  TWO5      TWOPOW (5)		/* 2^5     */
+#define  TWO8      TWOPOW (8)		/* 2^52    */
+#define  TWO10     TWOPOW (10)		/* 2^10    */
+#define  TWO18     TWOPOW (18)		/* 2^18    */
+#define  TWO19     TWOPOW (19)		/* 2^19    */
+#define  TWO23     TWOPOW (23)		/* 2^23    */
+
+#define  TWO57     0x1.0p57            /* 2^57    */
+#define  TWO71     0x1.0p71            /* 2^71    */
+#define  TWOM1032  0x1.0p-1032         /* 2^-1032 */
+#define  TWOM1022  0x1.0p-1022         /* 2^-1022 */
+
+#define  HALF      0x1.0p-1            /* 1/2 */
+#define  MHALF     -0x1.0p-1           /* -1/2 */
+#define  HALFRAD   0x1.0p23            /* 2^23 */
+
+int __acr (const mp_no *, const mp_no *, int);
+void __cpy (const mp_no *, mp_no *, int);
+void __mp_dbl (const mp_no *, double *, int);
+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);
 extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
 extern void __mpsqrt (mp_no *, mp_no *, int);
-extern void __mpexp (mp_no *, mp_no *__y, int);
+extern void __mpexp (mp_no *, mp_no *, int);
 extern void __c32 (mp_no *, mp_no *, mp_no *, int);
 extern int __mpranred (double, mp_no *, int);
+
+/* Given a power POW, build a multiprecision number 2^POW.  */
+static inline void
+__pow_mp (int pow, mp_no *y, int p)
+{
+  int i, rem;
+
+  /* The exponent is E such that E is a factor of 2^24.  The remainder (of the
+     form 2^x) goes entirely into the first digit of the mantissa as it is
+     always less than 2^24.  */
+  EY = pow / 24;
+  rem = pow - EY * 24;
+  EY++;
+
+  /* If the remainder is negative, it means that POW was negative since
+     |EY * 24| <= |pow|.  Adjust so that REM is positive and still less than
+     24 because of which, the mantissa digit is less than 2^24.  */
+  if (rem < 0)
+    {
+      EY--;
+      rem += 24;
+    }
+  /* The sign of any 2^x is always positive.  */
+  Y[0] = 1;
+  Y[1] = 1 << rem;
+
+  /* Everything else is 0.  */
+  for (i = 2; i <= p; i++)
+    Y[i] = 0;
+}
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpexp.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.c
@@ -69,13 +69,13 @@ __mpexp(mp_no *x, mp_no *y, int p) {
   for (i=0; i<EX; i++)  a *= RADIXI;
   for (   ; i>EX; i--)  a *= RADIX;
   b = X[1]*RADIXI;   m2 = 24*EX;
-  for (; b<HALF; m2--)  { a *= TWO;   b *= TWO; }
+  for (; b<HALF; m2--)  { a *= 2;   b *= 2; }
   if (b == HALF) {
-    for (i=2; i<=p; i++) { if (X[i]!=ZERO)  break; }
-    if (i==p+1)  { m2--;  a *= TWO; }
+    for (i=2; i<=p; i++) { if (X[i]!=0)  break; }
+    if (i==p+1)  { m2--;  a *= 2; }
   }
   if ((m=m1+m2) <= 0) {
-    m=0;  a=ONE;
+    m=0;  a=1;
     for (i=n-1; i>0; i--,n--) { if (m1np[i][p]+m2>0)  break; }
   }
 
@@ -84,8 +84,8 @@ __mpexp(mp_no *x, mp_no *y, int p) {
   __mul(x,&mpt1,&mps,p);
 
   /* Evaluate the polynomial. Put result in mpt2 */
-  mpone.e=1;   mpone.d[0]=ONE;   mpone.d[1]=ONE;
-  mpk.e = 1;   mpk.d[0] = ONE;   mpk.d[1]=__mpexp_nn[n].d;
+  mpone.e=1;   mpone.d[0]=1;   mpone.d[1]=1;
+  mpk.e = 1;   mpk.d[0] = 1;   mpk.d[1]=__mpexp_nn[n].d;
   __dvd(&mps,&mpk,&mpt1,p);
   __add(&mpone,&mpt1,&mpak,p);
   for (k=n-1; k>1; k--) {
@@ -99,9 +99,9 @@ __mpexp(mp_no *x, mp_no *y, int 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);  k++;
+    __sqr (&mpt2, &mpt1, p);  k++;
     if (k==m)  { j=1;  break; }
-    __mul(&mpt1,&mpt1,&mpt2,p);  k++;
+    __sqr (&mpt1, &mpt2, p);  k++;
   }
   if (j)  __cpy(&mpt1,y,p);
   else    __cpy(&mpt2,y,p);
Index: glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-avx.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/x86_64/fpu/multiarch/mpa-avx.c
+++ glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-avx.c
@@ -1,5 +1,6 @@
 #define __add __add_avx
 #define __mul __mul_avx
+#define __sqr __sqr_avx
 #define __sub __sub_avx
 #define __dbl_mp __dbl_mp_avx
 #define __dvd __dvd_avx
Index: glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c
+++ glibc-2.17-c758a686/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c
@@ -1,5 +1,6 @@
 #define __add __add_fma4
 #define __mul __mul_fma4
+#define __sqr __sqr_fma4
 #define __sub __sub_fma4
 #define __dbl_mp __dbl_mp_fma4
 #define __dvd __dvd_fma4
Index: glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/powerpc/power4/fpu/mpa.c
+++ glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa.c
@@ -2,7 +2,7 @@
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
- * Copyright (C) 2001, 2006 Free Software Foundation
+ * Copyright (C) 2001-2017 Free Software Foundation, Inc.
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU Lesser General Public License as published by
@@ -17,532 +17,198 @@
  * You should have received a copy of the GNU Lesser General Public License
  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  */
-/************************************************************************/
-/*  MODULE_NAME: mpa.c                                                  */
-/*                                                                      */
-/*  FUNCTIONS:                                                          */
-/*               mcr                                                    */
-/*               acr                                                    */
-/*               cr                                                     */
-/*               cpy                                                    */
-/*               cpymn                                                  */
-/*               norm                                                   */
-/*               denorm                                                 */
-/*               mp_dbl                                                 */
-/*               dbl_mp                                                 */
-/*               add_magnitudes                                         */
-/*               sub_magnitudes                                         */
-/*               add                                                    */
-/*               sub                                                    */
-/*               mul                                                    */
-/*               inv                                                    */
-/*               dvd                                                    */
-/*                                                                      */
-/* Arithmetic functions for multiple precision numbers.                 */
-/* Relative errors are bounded                                          */
-/************************************************************************/
-
-
-#include "endian.h"
-#include "mpa.h"
-#include "mpa2.h"
-#include <sys/param.h>	/* For MIN() */
-/* mcr() compares the sizes of the mantissas of two multiple precision  */
-/* numbers. Mantissas are compared regardless of the signs of the       */
-/* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also     */
-/* disregarded.                                                         */
-static int mcr(const mp_no *x, const mp_no *y, int p) {
-  long i;
-  long p2 = p;
-  for (i=1; i<=p2; i++) {
-    if      (X[i] == Y[i])  continue;
-    else if (X[i] >  Y[i])  return  1;
-    else                    return -1; }
-  return 0;
-}
-
-
-
-/* acr() compares the absolute values of two multiple precision numbers */
-int __acr(const mp_no *x, const mp_no *y, int p) {
-  long i;
-
-  if      (X[0] == ZERO) {
-    if    (Y[0] == ZERO) i= 0;
-    else                 i=-1;
-  }
-  else if (Y[0] == ZERO) i= 1;
-  else {
-    if      (EX >  EY)   i= 1;
-    else if (EX <  EY)   i=-1;
-    else                 i= mcr(x,y,p);
-  }
-
-  return i;
-}
-
-
-/* cr90 compares the values of two multiple precision numbers           */
-int  __cr(const mp_no *x, const mp_no *y, int p) {
-  int i;
-
-  if      (X[0] > Y[0])  i= 1;
-  else if (X[0] < Y[0])  i=-1;
-  else if (X[0] < ZERO ) i= __acr(y,x,p);
-  else                   i= __acr(x,y,p);
-
-  return i;
-}
-
-
-/* Copy a multiple precision number. Set *y=*x. x=y is permissible.      */
-void __cpy(const mp_no *x, mp_no *y, int p) {
-  long i;
-
-  EY = EX;
-  for (i=0; i <= p; i++)    Y[i] = X[i];
-
-  return;
-}
-
-
-/* Copy a multiple precision number x of precision m into a */
-/* multiple precision number y of precision n. In case n>m, */
-/* the digits of y beyond the m'th are set to zero. In case */
-/* n<m, the digits of x beyond the n'th are ignored.        */
-/* x=y is permissible.                                      */
-
-void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
-
-  long i,k;
-  long n2 = n;
-  long m2 = m;
-
-  EY = EX;     k=MIN(m2,n2);
-  for (i=0; i <= k; i++)    Y[i] = X[i];
-  for (   ; i <= n2; i++)    Y[i] = ZERO;
-
-  return;
-}
-
-/* Convert a multiple precision number *x into a double precision */
-/* number *y, normalized case  (|x| >= 2**(-1022))) */
-static void norm(const mp_no *x, double *y, int p)
-{
-  #define R  radixi.d
-  long i;
-#if 0
-  int k;
-#endif
-  double a,c,u,v,z[5];
-  if (p<5) {
-    if      (p==1) c = X[1];
-    else if (p==2) c = X[1] + R* X[2];
-    else if (p==3) c = X[1] + R*(X[2]  +   R* X[3]);
-    else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
-  }
-  else {
-    for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
-        {a *= TWO;   z[1] *= TWO; }
-
-    for (i=2; i<5; i++) {
-      z[i] = X[i]*a;
-      u = (z[i] + CUTTER)-CUTTER;
-      if  (u > z[i])  u -= RADIX;
-      z[i] -= u;
-      z[i-1] += u*RADIXI;
-    }
-
-    u = (z[3] + TWO71) - TWO71;
-    if (u > z[3])   u -= TWO19;
-    v = z[3]-u;
-
-    if (v == TWO18) {
-      if (z[4] == ZERO) {
-        for (i=5; i <= p; i++) {
-          if (X[i] == ZERO)   continue;
-          else                {z[3] += ONE;   break; }
-        }
-      }
-      else              z[3] += ONE;
-    }
-
-    c = (z[1] + R *(z[2] + R * z[3]))/a;
-  }
 
-  c *= X[0];
-
-  for (i=1; i<EX; i++)   c *= RADIX;
-  for (i=1; i>EX; i--)   c *= RADIXI;
-
-  *y = c;
-  return;
-#undef R
-}
-
-/* Convert a multiple precision number *x into a double precision */
-/* number *y, denormalized case  (|x| < 2**(-1022))) */
-static void denorm(const mp_no *x, double *y, int p)
+/* Define __mul and __sqr and use the rest from generic code.  */
+#define NO__MUL
+#define NO__SQR
+
+#include <sysdeps/ieee754/dbl-64/mpa.c>
+
+/* Multiply *X and *Y and store result in *Z.  X and Y may overlap but not X
+   and Z or Y and Z.  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.  */
+void
+__mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
-  long i,k;
+  long i, i1, i2, j, k, k2;
   long p2 = p;
-  double c,u,z[5];
-#if 0
-  double a,v;
-#endif
+  double u, zk, zk2;
 
-#define R  radixi.d
-  if (EX<-44 || (EX==-44 && X[1]<TWO5))
-     { *y=ZERO; return; }
-
-  if      (p2==1) {
-    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=ZERO;  z[3]=ZERO;  k=3;}
-    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=ZERO;  k=2;}
-    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
-  }
-  else if (p2==2) {
-    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  z[3]=ZERO;  k=3;}
-    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  z[3]=X[2];  k=2;}
-    else              {z[1]=     TWO10;  z[2]=ZERO;  z[3]=X[1];  k=1;}
-  }
-  else {
-    if      (EX==-42) {z[1]=X[1]+TWO10;  z[2]=X[2];  k=3;}
-    else if (EX==-43) {z[1]=     TWO10;  z[2]=X[1];  k=2;}
-    else              {z[1]=     TWO10;  z[2]=ZERO;  k=1;}
-    z[3] = X[k];
-  }
-
-  u = (z[3] + TWO57) - TWO57;
-  if  (u > z[3])   u -= TWO5;
-
-  if (u==z[3]) {
-    for (i=k+1; i <= p2; i++) {
-      if (X[i] == ZERO)   continue;
-      else {z[3] += ONE;   break; }
+  /* Is z=0?  */
+  if (__glibc_unlikely (X[0] * Y[0] == 0))
+    {
+      Z[0] = 0;
+      return;
     }
-  }
-
-  c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
 
-  *y = c*TWOM1032;
-  return;
-
-#undef R
-}
-
-/* Convert a multiple precision number *x into a double precision number *y. */
-/* The result is correctly rounded to the nearest/even. *x is left unchanged */
-
-void __mp_dbl(const mp_no *x, double *y, int p) {
-#if 0
-  int i,k;
-  double a,c,u,v,z[5];
+  /* Multiply, add and carry */
+  k2 = (p2 < 3) ? p2 + p2 : p2 + 3;
+  zk = Z[k2] = 0;
+  for (k = k2; k > 1;)
+    {
+      if (k > p2)
+	{
+	  i1 = k - p2;
+	  i2 = p2 + 1;
+	}
+      else
+	{
+	  i1 = 1;
+	  i2 = k;
+	}
+#if 1
+      /* Rearrange this inner loop to allow the fmadd instructions to be
+         independent and execute in parallel on processors that have
+         dual symmetrical FP pipelines.  */
+      if (i1 < (i2 - 1))
+	{
+	  /* Make sure we have at least 2 iterations.  */
+	  if (((i2 - i1) & 1L) == 1L)
+	    {
+	      /* Handle the odd iterations case.  */
+	      zk2 = x->d[i2 - 1] * y->d[i1];
+	    }
+	  else
+	    zk2 = 0.0;
+	  /* Do two multiply/adds per loop iteration, using independent
+	     accumulators; zk and zk2.  */
+	  for (i = i1, j = i2 - 1; i < i2 - 1; i += 2, j -= 2)
+	    {
+	      zk += x->d[i] * y->d[j];
+	      zk2 += x->d[i + 1] * y->d[j - 1];
+	    }
+	  zk += zk2;		/* Final sum.  */
+	}
+      else
+	{
+	  /* Special case when iterations is 1.  */
+	  zk += x->d[i1] * y->d[i1];
+	}
+#else
+      /* The original code.  */
+      for (i = i1, j = i2 - 1; i < i2; i++, j--)
+	zk += X[i] * Y[j];
 #endif
 
-  if (X[0] == ZERO)  {*y = ZERO;  return; }
-
-  if      (EX> -42)                 norm(x,y,p);
-  else if (EX==-42 && X[1]>=TWO10)  norm(x,y,p);
-  else                              denorm(x,y,p);
-}
-
-
-/* dbl_mp() converts a double precision number x into a multiple precision  */
-/* number *y. If the precision p is too small the result is truncated. x is */
-/* left unchanged.                                                          */
-
-void __dbl_mp(double x, mp_no *y, int p) {
+      u = (zk + CUTTER) - CUTTER;
+      if (u > zk)
+	u -= RADIX;
+      Z[k] = zk - u;
+      zk = u * RADIXI;
+      --k;
+    }
+  Z[k] = zk;
 
-  long i,n;
-  long p2 = p;
-  double u;
+  int e = EX + EY;
+  /* Is there a carry beyond the most significant digit?  */
+  if (Z[1] == 0)
+    {
+      for (i = 1; i <= p2; i++)
+	Z[i] = Z[i + 1];
+      e--;
+    }
 
-  /* Sign */
-  if      (x == ZERO)  {Y[0] = ZERO;  return; }
-  else if (x >  ZERO)   Y[0] = ONE;
-  else                 {Y[0] = MONE;  x=-x;   }
-
-  /* Exponent */
-  for (EY=ONE; x >= RADIX; EY += ONE)   x *= RADIXI;
-  for (      ; x <  ONE;   EY -= ONE)   x *= RADIX;
-
-  /* Digits */
-  n=MIN(p2,4);
-  for (i=1; i<=n; i++) {
-    u = (x + TWO52) - TWO52;
-    if (u>x)   u -= ONE;
-    Y[i] = u;     x -= u;    x *= RADIX; }
-  for (   ; i<=p2; i++)     Y[i] = ZERO;
-  return;
+  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
+__sqr (const mp_no *x, mp_no *y, int p)
+{
+  long i, j, k, ip;
+  double u, yk;
 
-/*  add_magnitudes() adds the magnitudes of *x & *y assuming that           */
-/*  abs(*x) >= abs(*y) > 0.                                                 */
-/* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
-/* No guard digit is used. The result equals the exact sum, truncated.      */
-/* *x & *y are left unchanged.                                              */
+  /* Is z=0?  */
+  if (__glibc_unlikely (X[0] == 0))
+    {
+      Y[0] = 0;
+      return;
+    }
 
-static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+  /* We need not iterate through all X's since it's pointless to
+     multiply zeroes.  */
+  for (ip = p; ip > 0; ip--)
+    if (X[ip] != 0)
+      break;
 
-  long i,j,k;
-  long p2 = p;
+  k = (__glibc_unlikely (p < 3)) ? p + p : p + 3;
 
-  EZ = EX;
+  while (k > 2 * ip + 1)
+    Y[k--] = 0;
 
-  i=p2;    j=p2+ EY - EX;    k=p2+1;
+  yk = 0;
 
-  if (j<1)
-     {__cpy(x,z,p);  return; }
-  else   Z[k] = ZERO;
-
-  for (; j>0; i--,j--) {
-    Z[k] += X[i] + Y[j];
-    if (Z[k] >= RADIX) {
-      Z[k]  -= RADIX;
-      Z[--k] = ONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  for (; i>0; i--) {
-    Z[k] += X[i];
-    if (Z[k] >= RADIX) {
-      Z[k]  -= RADIX;
-      Z[--k] = ONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  if (Z[1] == ZERO) {
-    for (i=1; i<=p2; i++)    Z[i] = Z[i+1]; }
-  else   EZ += ONE;
-}
+  while (k > p)
+    {
+      double yk2 = 0.0;
+      long lim = k / 2;
 
+      if (k % 2 == 0)
+        {
+	  yk += X[lim] * X[lim];
+	  lim--;
+	}
 
-/*  sub_magnitudes() subtracts the magnitudes of *x & *y assuming that      */
-/*  abs(*x) > abs(*y) > 0.                                                  */
-/* The sign of the difference *z is undefined. x&y may overlap but not x&z  */
-/* or y&z. One guard digit is used. The error is less than one ulp.         */
-/* *x & *y are left unchanged.                                              */
+      /* In __mul, this loop (and the one within the next while loop) run
+         between a range to calculate the mantissa as follows:
 
-static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+         Z[k] = X[k] * Y[n] + X[k+1] * Y[n-1] ... + X[n-1] * Y[k+1]
+		+ X[n] * Y[k]
 
-  long i,j,k;
-  long p2 = p;
+         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 <= lim; i++, j--)
+	yk2 += X[i] * X[j];
 
-  EZ = EX;
+      yk += 2.0 * yk2;
 
-  if (EX == EY) {
-    i=j=k=p2;
-    Z[k] = Z[k+1] = ZERO; }
-  else {
-    j= EX - EY;
-    if (j > p2)  {__cpy(x,z,p);  return; }
-    else {
-      i=p2;   j=p2+1-j;   k=p2;
-      if (Y[j] > ZERO) {
-        Z[k+1] = RADIX - Y[j--];
-        Z[k]   = MONE; }
-      else {
-        Z[k+1] = ZERO;
-        Z[k]   = ZERO;   j--;}
+      u = (yk + CUTTER) - CUTTER;
+      if (u > yk)
+	u -= RADIX;
+      Y[k--] = yk - u;
+      yk = u * RADIXI;
     }
-  }
-
-  for (; j>0; i--,j--) {
-    Z[k] += (X[i] - Y[j]);
-    if (Z[k] < ZERO) {
-      Z[k]  += RADIX;
-      Z[--k] = MONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  for (; i>0; i--) {
-    Z[k] += X[i];
-    if (Z[k] < ZERO) {
-      Z[k]  += RADIX;
-      Z[--k] = MONE; }
-    else
-      Z[--k] = ZERO;
-  }
-
-  for (i=1; Z[i] == ZERO; i++) ;
-  EZ = EZ - i + 1;
-  for (k=1; i <= p2+1; )
-    Z[k++] = Z[i++];
-  for (; k <= p2; )
-    Z[k++] = ZERO;
-
-  return;
-}
-
-
-/* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap  */
-/* but not x&z or y&z. One guard digit is used. The error is less than    */
-/* one ulp. *x & *y are left unchanged.                                   */
-
-void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
-  int n;
-
-  if      (X[0] == ZERO)     {__cpy(y,z,p);  return; }
-  else if (Y[0] == ZERO)     {__cpy(x,z,p);  return; }
-
-  if (X[0] == Y[0])   {
-    if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] = X[0]; }
-    else                     {add_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
-  }
-  else                       {
-    if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] = X[0]; }
-    else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = Y[0]; }
-    else                      Z[0] = ZERO;
-  }
-  return;
-}
-
-
-/* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
-/* overlap but not x&z or y&z. One guard digit is used. The error is      */
-/* less than one ulp. *x & *y are left unchanged.                         */
-
-void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
-  int n;
-
-  if      (X[0] == ZERO)     {__cpy(y,z,p);  Z[0] = -Z[0];  return; }
-  else if (Y[0] == ZERO)     {__cpy(x,z,p);                 return; }
-
-  if (X[0] != Y[0])    {
-    if (__acr(x,y,p) > 0)      {add_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
-    else                     {add_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
-  }
-  else                       {
-    if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p);  Z[0] =  X[0]; }
-    else if (n == -1)        {sub_magnitudes(y,x,z,p);  Z[0] = -Y[0]; }
-    else                      Z[0] = ZERO;
-  }
-  return;
-}
-
-
-/* Multiply two multiple precision numbers. *z is set to *x * *y. x&y      */
-/* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is     */
-/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp.   */
-/* *x & *y are left unchanged.                                             */
 
-void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
-  long i, i1, i2, j, k, k2;
-  long p2 = p;
-  double u, zk, zk2;
-
-                      /* Is z=0? */
-  if (X[0]*Y[0]==ZERO)
-     { Z[0]=ZERO;  return; }
-
-                       /* Multiply, add and carry */
-  k2 = (p2<3) ? p2+p2 : p2+3;
-  zk = Z[k2]=ZERO;
-  for (k=k2; k>1; ) {
-    if (k > p2)  {i1=k-p2; i2=p2+1; }
-    else        {i1=1;   i2=k;   }
-#if 1
-    /* rearange this inner loop to allow the fmadd instructions to be
-       independent and execute in parallel on processors that have
-       dual symetrical FP pipelines.  */
-    if (i1 < (i2-1))
+  while (k > 1)
     {
-	/* make sure we have at least 2 iterations */
-	if (((i2 - i1) & 1L) == 1L)
-	{
-                /* Handle the odd iterations case.  */
-		zk2 = x->d[i2-1]*y->d[i1];
-	}
-	else
-		zk2 = zero.d;
-	/* Do two multiply/adds per loop iteration, using independent
-           accumulators; zk and zk2.  */
-	for (i=i1,j=i2-1; i<i2-1; i+=2,j-=2) 
-	{
-		zk += x->d[i]*y->d[j];
-		zk2 += x->d[i+1]*y->d[j-1];
+      double yk2 = 0.0;
+      long lim = k / 2;
+
+      if (k % 2 == 0)
+        {
+	  yk += X[lim] * X[lim];
+	  lim--;
 	}
-	zk += zk2; /* final sum.  */
-    }
-    else
+
+      /* Likewise for this loop.  */
+      for (i = 1, j = k - 1; i <= lim; 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;
+
+  int e = EX * 2;
+  /* Is there a carry beyond the most significant digit?  */
+  if (__glibc_unlikely (Y[1] == 0))
     {
-        /* Special case when iterations is 1.  */
-	zk += x->d[i1]*y->d[i1];
+      for (i = 1; i <= p; i++)
+	Y[i] = Y[i + 1];
+      e--;
     }
-#else
-    /* The orginal code.  */
-    for (i=i1,j=i2-1; i<i2; i++,j--)  zk += X[i]*Y[j];
-#endif
-
-    u = (zk + CUTTER)-CUTTER;
-    if  (u > zk)  u -= RADIX;
-    Z[k]  = zk - u;
-    zk = u*RADIXI;
-    --k;
-  }
-  Z[k] = zk;
-
-                 /* Is there a carry beyond the most significant digit? */
-  if (Z[1] == ZERO) {
-    for (i=1; i<=p2; i++)  Z[i]=Z[i+1];
-    EZ = EX + EY - 1; }
-  else
-    EZ = EX + EY;
-
-  Z[0] = X[0] * Y[0];
-  return;
-}
-
-
-/* Invert a multiple precision number. Set *y = 1 / *x.                     */
-/* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3,   */
-/* 2.001*r**(1-p) for p>3.                                                  */
-/* *x=0 is not permissible. *x is left unchanged.                           */
-
-void __inv(const mp_no *x, mp_no *y, int p) {
-  long i;
-#if 0
-  int l;
-#endif
-  double t;
-  mp_no z,w;
-  static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
-                            4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
-  const mp_no mptwo = {1,{1.0,2.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}};
-
-  __cpy(x,&z,p);  z.e=0;  __mp_dbl(&z,&t,p);
-  t=ONE/t;   __dbl_mp(t,y,p);    EY -= EX;
-
-  for (i=0; i<np1[p]; i++) {
-    __cpy(y,&w,p);
-    __mul(x,&w,y,p);
-    __sub(&mptwo,y,&z,p);
-    __mul(&w,&z,y,p);
-  }
-  return;
-}
-
-
-/* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
-/* are left unchanged. x&y may overlap but not x&z or y&z.                   */
-/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3     */
-/* and 3.001*r**(1-p) for p>3. *y=0 is not permissible.                      */
-
-void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
-
-  mp_no w;
-
-  if (X[0] == ZERO)    Z[0] = ZERO;
-  else                {__inv(y,&w,p);   __mul(x,&w,z,p);}
-  return;
+  EY = e;
 }
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/atnat.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat.h
@@ -138,8 +138,6 @@
 #endif
 #endif
 
-#define  ZERO      zero.d
-#define  ONE       one.d
 #define  A         a.d
 #define  B         b.d
 #define  C         c.d
@@ -160,7 +158,5 @@
 #define  U6        u6.d
 #define  U7        u7.d
 #define  U8        u8.d
-#define  TWO8      two8.d
-#define  TWO52     two52.d
 
 #endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat2.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/atnat2.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/atnat2.h
@@ -174,11 +174,4 @@
 #endif
 #endif
 
-#define  ZERO      zero.d
-#define  MZERO     mzero.d
-#define  ONE       one.d
-#define  TWO8      two8.d
-#define  TWO52     two52.d
-#define  TWOM1022  twom1022.d
-
 #endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa2.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpa2.h
+++ /dev/null
@@ -1,94 +0,0 @@
-
-/*
- * IBM Accurate Mathematical Library
- * Written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation, Inc.
- *
- * This program 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.
- *
- * This program 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 this program; if not, see <http://www.gnu.org/licenses/>.
- */
-
-/**************************************************************************/
-/*                                                                        */
-/* MODULE_NAME:mpa2.h                                                     */
-/*                                                                        */
-/*                                                                        */
-/*   variables prototype and definition   according to type of processor  */
-/*   types definition                                                     */
-/**************************************************************************/
-
-#ifndef MPA2_H
-#define MPA2_H
-
-
-#ifdef BIG_ENDI
-static const number
-/**/ radix          = {{0x41700000, 0x00000000} }, /* 2**24  */
-/**/ radixi         = {{0x3e700000, 0x00000000} }, /* 2**-24 */
-/**/ cutter         = {{0x44b00000, 0x00000000} }, /* 2**76  */
-/**/ zero           = {{0x00000000, 0x00000000} }, /*  0     */
-/**/ one            = {{0x3ff00000, 0x00000000} }, /*  1     */
-/**/ mone           = {{0xbff00000, 0x00000000} }, /* -1     */
-/**/ two            = {{0x40000000, 0x00000000} }, /*  2     */
-/**/ two5           = {{0x40400000, 0x00000000} }, /* 2**5   */
-/**/ two10          = {{0x40900000, 0x00000000} }, /* 2**10  */
-/**/ two18          = {{0x41100000, 0x00000000} }, /* 2**18  */
-/**/ two19          = {{0x41200000, 0x00000000} }, /* 2**19  */
-/**/ two23          = {{0x41600000, 0x00000000} }, /* 2**23  */
-/**/ two52          = {{0x43300000, 0x00000000} }, /* 2**52  */
-/**/ two57          = {{0x43800000, 0x00000000} }, /* 2**57  */
-/**/ two71          = {{0x44600000, 0x00000000} }, /* 2**71  */
-/**/ twom1032       = {{0x00000400, 0x00000000} }; /* 2**-1032 */
-
-#else
-#ifdef LITTLE_ENDI
-static const number
-/**/ radix          = {{0x00000000, 0x41700000} }, /* 2**24  */
-/**/ radixi         = {{0x00000000, 0x3e700000} }, /* 2**-24 */
-/**/ cutter         = {{0x00000000, 0x44b00000} }, /* 2**76  */
-/**/ zero           = {{0x00000000, 0x00000000} }, /*  0     */
-/**/ one            = {{0x00000000, 0x3ff00000} }, /*  1     */
-/**/ mone           = {{0x00000000, 0xbff00000} }, /* -1     */
-/**/ two            = {{0x00000000, 0x40000000} }, /*  2     */
-/**/ two5           = {{0x00000000, 0x40400000} }, /* 2**5   */
-/**/ two10          = {{0x00000000, 0x40900000} }, /* 2**10  */
-/**/ two18          = {{0x00000000, 0x41100000} }, /* 2**18  */
-/**/ two19          = {{0x00000000, 0x41200000} }, /* 2**19  */
-/**/ two23          = {{0x00000000, 0x41600000} }, /* 2**23  */
-/**/ two52          = {{0x00000000, 0x43300000} }, /* 2**52  */
-/**/ two57          = {{0x00000000, 0x43800000} }, /* 2**57  */
-/**/ two71          = {{0x00000000, 0x44600000} }, /* 2**71  */
-/**/ twom1032       = {{0x00000000, 0x00000400} }; /* 2**-1032 */
-
-#endif
-#endif
-
-#define  RADIX     radix.d
-#define  RADIXI    radixi.d
-#define  CUTTER    cutter.d
-#define  ZERO      zero.d
-#define  ONE       one.d
-#define  MONE      mone.d
-#define  TWO       two.d
-#define  TWO5      two5.d
-#define  TWO10     two10.d
-#define  TWO18     two18.d
-#define  TWO19     two19.d
-#define  TWO23     two23.d
-#define  TWO52     two52.d
-#define  TWO57     two57.d
-#define  TWO71     two71.d
-#define  TWOM1032  twom1032.d
-
-
-#endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpatan.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.h
@@ -177,6 +177,3 @@ __atan_twonm1[33] = {
 
 #endif
 #endif
-
-#define  ONE       __atan_one.d
-#define  TWO       __atan_two.d
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan2.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpatan2.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan2.c
@@ -49,18 +49,16 @@ void
 SECTION
 __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) {
 
-  static const double ZERO = 0.0, ONE = 1.0;
-
   mp_no mpone = {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 mpt1,mpt2,mpt3;
 
 
-  if (X[0] <= ZERO) {
-    mpone.e = 1;                 mpone.d[0] = mpone.d[1] = ONE;
+  if (X[0] <= 0) {
+    mpone.e = 1;                 mpone.d[0] = mpone.d[1] = 1;
     __dvd(x,y,&mpt1,p);          __mul(&mpt1,&mpt1,&mpt2,p);
-    if (mpt1.d[0] != ZERO)       mpt1.d[0] = ONE;
+    if (mpt1.d[0] != 0)       mpt1.d[0] = 1;
     __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);
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpexp.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpexp.h
@@ -159,11 +159,4 @@ extern const number __mpexp_half attribu
 #endif
 #endif
 
-#define  RADIX     __mpexp_radix.d
-#define  RADIXI    __mpexp_radixi.d
-#define  ZERO      __mpexp_zero.d
-#define  ONE       __mpexp_one.d
-#define  TWO       __mpexp_two.d
-#define  HALF      __mpexp_half.d
-
 #endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpsqrt.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.h
@@ -51,7 +51,4 @@ extern const int __mpsqrt_mp[33] attribu
 			     4,4,4,4,4,4,4,4,4};
 #endif
 
-#define  ONE       __mpsqrt_one.d
-#define  HALFRAD   __mpsqrt_halfrad.d
-
 #endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mptan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mptan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mptan.c
@@ -47,8 +47,6 @@ void
 SECTION
 __mptan(double x, mp_no *mpy, int p) {
 
-  static const double MONE = -1.0;
-
   int n;
   mp_no mpw, mpc, mps;
 
@@ -56,7 +54,7 @@ __mptan(double x, mp_no *mpy, int p) {
   __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;
+    mpy->d[0] *= -1;
   }                          /* tan is negative in this area */
   else  __dvd(&mps,&mpc,mpy,p);
 
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/ulog.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/ulog.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/ulog.h
@@ -181,10 +181,6 @@
 #endif
 #endif
 
-#define  ZERO      zero.d
-#define  ONE       one.d
-#define  HALF      half.d
-#define  MHALF     mhalf.d
 #define  SQRT_2    sqrt_2.d
 #define  DEL_U     delu.d
 #define  DEL_V     delv.d
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/utan.h
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/utan.h
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/utan.h
@@ -270,10 +270,4 @@
 #endif
 #endif
 
-
-#define  ZERO      zero.d
-#define  ONE       one.d
-#define  MONE      mone.d
-#define  TWO8      two8.d
-
 #endif
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa-arch.h
===================================================================
--- /dev/null
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpa-arch.h
@@ -0,0 +1,47 @@
+/* Overridable constants and operations.
+   Copyright (C) 2013 Free Software Foundation, Inc.
+
+   This program 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.
+
+   This program 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 this program; if not, see <http://www.gnu.org/licenses/>.  */
+
+#include <stdint.h>
+
+typedef long mantissa_t;
+typedef int64_t mantissa_store_t;
+
+#define TWOPOW(i) (1L << i)
+
+#define RADIX_EXP 24
+#define  RADIX TWOPOW (RADIX_EXP)		/* 2^24    */
+
+/* Divide D by RADIX and put the remainder in R.  D must be a non-negative
+   integral value.  */
+#define DIV_RADIX(d, r) \
+  ({									      \
+    r = d & (RADIX - 1);						      \
+    d >>= RADIX_EXP;							      \
+  })
+
+/* Put the integer component of a double X in R and retain the fraction in
+   X.  This is used in extracting mantissa digits for MP_NO by using the
+   integer portion of the current value of the number as the current mantissa
+   digit and then scaling by RADIX to get the next mantissa digit in the same
+   manner.  */
+#define INTEGER_OF(x, i) \
+  ({									      \
+    i = (mantissa_t) x;							      \
+    x -= i;								      \
+  })
+
+/* Align IN down to F.  The code assumes that F is a power of two.  */
+#define ALIGN_DOWN_TO(in, f) ((in) & -(f))
Index: glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa-arch.h
===================================================================
--- /dev/null
+++ glibc-2.17-c758a686/sysdeps/powerpc/power4/fpu/mpa-arch.h
@@ -0,0 +1,56 @@
+/* Overridable constants and operations.
+   Copyright (C) 2013-2017 Free Software Foundation, Inc.
+
+   This program 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.
+
+   This program 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 this program; if not, see <http://www.gnu.org/licenses/>.  */
+
+typedef double mantissa_t;
+typedef double mantissa_store_t;
+
+#define TWOPOW(i) (0x1.0p##i)
+
+#define RADIX TWOPOW (24)		/* 2^24    */
+#define CUTTER TWOPOW (76)		/* 2^76    */
+#define RADIXI 0x1.0p-24		/* 2^-24 */
+#define TWO52 TWOPOW (52)		/* 2^52 */
+
+/* Divide D by RADIX and put the remainder in R.  */
+#define DIV_RADIX(d,r) \
+  ({									      \
+    double u = ((d) + CUTTER) - CUTTER;					      \
+    if (u > (d))							      \
+      u -= RADIX;							      \
+    r = (d) - u;							      \
+    (d) = u * RADIXI;							      \
+  })
+
+/* Put the integer component of a double X in R and retain the fraction in
+   X.  */
+#define INTEGER_OF(x, r) \
+  ({									      \
+    double u = ((x) + TWO52) - TWO52;					      \
+    if (u > (x))							      \
+      u -= 1;								      \
+    (r) = u;								      \
+    (x) -= u;								      \
+  })
+
+/* Align IN down to a multiple of F, where F is a power of two.  */
+#define ALIGN_DOWN_TO(in, f) \
+  ({									      \
+    double factor = f * TWO52;						      \
+    double u = (in + factor) - factor;					      \
+    if (u > in)								      \
+      u -= f;								      \
+    u;									      \
+  })
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_atan2.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/e_atan2.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -98,15 +98,15 @@ __ieee754_atan2(double y,double x) {
   /* y=+-0 */
   if      (uy==0x00000000) {
     if    (dy==0x00000000) {
-      if  ((ux&0x80000000)==0x00000000)  return ZERO;
+      if  ((ux&0x80000000)==0x00000000)  return 0;
       else                               return opi.d; } }
   else if (uy==0x80000000) {
     if    (dy==0x00000000) {
-      if  ((ux&0x80000000)==0x00000000)  return MZERO;
+      if  ((ux&0x80000000)==0x00000000)  return -0.0;
       else                               return mopi.d;} }
 
   /* x=+-0 */
-  if (x==ZERO) {
+  if (x==0) {
     if ((uy&0x80000000)==0x00000000)     return hpi.d;
     else                                 return mhpi.d; }
 
@@ -118,8 +118,8 @@ __ieee754_atan2(double y,double x) {
       else if (uy==0xfff00000) {
 	if    (dy==0x00000000)  return mqpi.d; }
       else {
-	if    ((uy&0x80000000)==0x00000000)  return ZERO;
-	else                                 return MZERO; }
+	if    ((uy&0x80000000)==0x00000000)  return 0;
+	else                                 return -0.0; }
     }
   }
   else if     (ux==0xfff00000) {
@@ -141,14 +141,14 @@ __ieee754_atan2(double y,double x) {
     if    (dy==0x00000000)  return mhpi.d; }
 
   /* either x/y or y/x is very close to zero */
-  ax = (x<ZERO) ? -x : x;    ay = (y<ZERO) ? -y : y;
+  ax = (x<0) ? -x : x;    ay = (y<0) ? -y : y;
   de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
-  if      (de>=ep)  { return ((y>ZERO) ? hpi.d : mhpi.d); }
+  if      (de>=ep)  { return ((y>0) ? hpi.d : mhpi.d); }
   else if (de<=em)  {
-    if    (x>ZERO)  {
+    if    (x>0)  {
       if  ((z=ay/ax)<TWOM1022)  return normalized(ax,ay,y,z);
       else                      return signArctan2(y,z); }
-    else            { return ((y>ZERO) ? opi.d : mopi.d); } }
+    else            { return ((y>0) ? opi.d : mopi.d); } }
 
   /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
   if (ax<twom500.d || ay<twom500.d) { ax*=two500.d;  ay*=two500.d; }
@@ -170,7 +170,7 @@ __ieee754_atan2(double y,double x) {
     EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
     du=((ax-v)-vv)/ay; }
 
-  if (x>ZERO) {
+  if (x>0) {
 
     /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
     if (ay<ax) {
@@ -180,7 +180,7 @@ __ieee754_atan2(double y,double x) {
 
 	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(f9.d,ff9.d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -212,7 +212,7 @@ __ieee754_atan2(double y,double x) {
 	EADD(t1,du,v,vv)
 	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 	   v*(hij[i][14].d+v* hij[i][15].d))));
-	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -237,7 +237,7 @@ __ieee754_atan2(double y,double x) {
 
 	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(f9.d,ff9.d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -265,7 +265,7 @@ __ieee754_atan2(double y,double x) {
 	EADD(t1,du,v,vv)
 	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 	   v*(hij[i][14].d+v* hij[i][15].d))));
-	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -293,7 +293,7 @@ __ieee754_atan2(double y,double x) {
 
 	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(f9.d,ff9.d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -321,7 +321,7 @@ __ieee754_atan2(double y,double x) {
 	EADD(t1,du,v,vv)
 	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 	   v*(hij[i][14].d+v* hij[i][15].d))));
-	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -347,7 +347,7 @@ __ieee754_atan2(double y,double x) {
 
 	MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(f9.d,ff9.d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -375,7 +375,7 @@ __ieee754_atan2(double y,double x) {
 	EADD(t1,du,v,vv)
 	s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
 	   v*(hij[i][14].d+v* hij[i][15].d))));
-	ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_log.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/e_log.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/e_log.c
@@ -78,9 +78,9 @@ __ieee754_log(double x) {
   n=0;
   if (__builtin_expect(ux < 0x00100000, 0)) {
     if (__builtin_expect(((ux & 0x7fffffff) | dx) == 0, 0))
-      return MHALF/ZERO; /* return -INF */
+      return MHALF/0; /* return -INF */
     if (__builtin_expect(ux < 0, 0))
-      return (x-x)/ZERO;                         /* return NaN  */
+      return (x-x)/0;                         /* return NaN  */
     n -= 54;    x *= two54.d;                              /* scale x     */
     num.d = x;
   }
@@ -89,7 +89,7 @@ __ieee754_log(double x) {
 
   /* Regular values of x */
 
-  w = x-ONE;
+  w = x-1;
   if (__builtin_expect(ABS(w) > U03, 1)) { goto case_03; }
 
 
@@ -113,25 +113,25 @@ __ieee754_log(double x) {
 	    w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
   EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
   ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
   ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
-  MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-  MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
-  ADD2(w,ZERO,    s3,ss3, b, bb,t1,t2)
+  MUL2(w,0,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(w,0,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
+  ADD2(w,0,    s3,ss3, b, bb,t1,t2)
 
   /* End stage II, case abs(x-1) < 0.03 */
   if ((y=b+(bb+b*E4)) == b+(bb-b*E4))  return y;
@@ -155,7 +155,7 @@ __ieee754_log(double x) {
   j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
 
   /* Compute w=(u-ui*vj)/(ui*vj) */
-  p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
+  p0=(1+(i-75)*DEL_U)*(1+(j-180)*DEL_V);
   q=u-p0;   r0=Iu[i].d*Iv[j].d;   w=q*r0;
 
   /* Evaluate polynomial I */
@@ -178,11 +178,11 @@ __ieee754_log(double x) {
 
   /* Improve the accuracy of r0 */
   EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
-  t=r0*((ONE-sa)-sb);
+  t=r0*((1-sa)-sb);
   EADD(r0,t,ra,rb)
 
   /* Compute w */
-  MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
+  MUL2(q,0,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
 
   EADD(A,B0,a0,aa0)
 
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_atan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/s_atan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_atan.c
@@ -82,7 +82,7 @@ double atan(double x) {
     return x+x;
 
   /* Regular values of x, including denormals +-0 and +-INF */
-  u = (x<ZERO) ? -x : x;
+  u = (x<0) ? -x : x;
   if (u<C) {
     if (u<B) {
       if (u<A) {                                           /* u < A */
@@ -93,7 +93,7 @@ double atan(double x) {
 
 	EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
 	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(f9.d,ff9.d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -101,8 +101,8 @@ double atan(double x) {
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
-	MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
-	ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
+	MUL2(x,0,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+	ADD2(x,0,s2,ss2,s1,ss1,t1,t2)
 	if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
 
 	return atanMp(x,pr);
@@ -124,14 +124,14 @@ double atan(double x) {
       z=u-hij[i][0].d;
       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
 	 z*(hij[i][14].d+z* hij[i][15].d))));
-      ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
-      MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+      ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
+      MUL2(z,0,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
-      MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+      MUL2(z,0,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
-      MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+      MUL2(z,0,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
-      MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+      MUL2(z,0,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
       if ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2))  return __signArctan(x,y);
 
@@ -140,9 +140,9 @@ double atan(double x) {
   }
   else {
     if (u<D) { /* C <= u < D */
-      w=ONE/u;
+      w=1/u;
       EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
-      ww=w*((ONE-t1)-t2);
+      ww=w*((1-t1)-t2);
       i=(TWO52+TWO8*w)-TWO52;  i-=16;
       z=(w-cij[i][0].d)+ww;
       yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
@@ -152,12 +152,12 @@ double atan(double x) {
       else        u3=U32;  /* w >= 1/2 */
       if ((y=t1+(yy-u3)) == t1+(yy+u3))  return __signArctan(x,y);
 
-      DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+      DIV2(1,0,u,0,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
       t1=w-hij[i][0].d;
       EADD(t1,ww,z,zz)
       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
 	 z*(hij[i][14].d+z* hij[i][15].d))));
-      ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+      ADD2(hij[i][9].d,hij[i][10].d,s1,0,s2,ss2,t1,t2)
       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
@@ -173,18 +173,18 @@ double atan(double x) {
     }
     else {
       if (u<E) { /* D <= u < E */
-	w=ONE/u;   v=w*w;
+	w=1/u;   v=w*w;
 	EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
 	yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
-	ww=w*((ONE-t1)-t2);
+	ww=w*((1-t1)-t2);
 	ESUB(HPI,w,t3,cor)
 	yy=((HPI1+cor)-ww)-yy;
 	if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
 
-	DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+	DIV2(1,0,u,0,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 	MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 	s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
-	ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+	ADD2(f9.d,ff9.d,s1,0,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 	ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 	MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_tan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/s_tan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/s_tan.c
@@ -84,7 +84,7 @@ tan(double x) {
     goto ret;
   }
 
-  w=(x<ZERO) ? -x : x;
+  w=(x<0) ? -x : x;
 
   /* (I) The case abs(x) <= 1.259e-8 */
   if (w<=g1.d) { retval = x; goto ret; }
@@ -125,11 +125,11 @@ tan(double x) {
 
     /* First stage */
     i = ((int) (mfftnhf.d+TWO8*w));
-    z = w-xfg[i][0].d;  z2 = z*z;   s = (x<ZERO) ? MONE : ONE;
+    z = w-xfg[i][0].d;  z2 = z*z;   s = (x<0) ? -1 : 1;
     pz = z+z*z2*(e0.d+z2*e1.d);
     fi = xfg[i][1].d;   gi = xfg[i][2].d;   t2 = pz*(gi+fi)/(gi-pz);
     if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d)) { retval = (s*y); goto ret; }
-    t3 = (t2<ZERO) ? -t2 : t2;
+    t3 = (t2<0) ? -t2 : t2;
     t4 = fi*ua3.d+t3*ub3.d;
     if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (s*y); goto ret; }
 
@@ -165,8 +165,8 @@ tan(double x) {
     da = xn*mp3.d;
     a=t1-da;
     da = (t1-a)-da;
-    if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
-    else         {ya= a;  yya= da;  sy= ONE;}
+    if (a<0)  {ya=-a;  yya=-da;  sy=-1;}
+    else         {ya= a;  yya= da;  sy= 1;}
 
     /* (IV),(V) The case 0.787 < abs(x) <= 25,    abs(y) <= 1e-7 */
     if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
@@ -240,14 +240,14 @@ tan(double x) {
       /* -cot */
       t2 = pz*(fi+gi)/(fi+pz);
       if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d)) { retval = (-sy*y); goto ret; }
-      t3 = (t2<ZERO) ? -t2 : t2;
+      t3 = (t2<0) ? -t2 : t2;
       t4 = gi*ua10.d+t3*ub10.d;
       if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
     else   {
       /* tan */
       t2 = pz*(gi+fi)/(gi-pz);
       if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d)) { retval = (sy*y); goto ret; }
-      t3 = (t2<ZERO) ? -t2 : t2;
+      t3 = (t2<0) ? -t2 : t2;
       t4 = fi*ua9.d+t3*ub9.d;
       if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
 
@@ -295,8 +295,8 @@ tan(double x) {
     a = t - t1;
     da = ((t-a)-t1)+da;
     EADD(a,da,t1,t2)   a=t1;  da=t2;
-    if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
-    else         {ya= a;  yya= da;  sy= ONE;}
+    if (a<0)  {ya=-a;  yya=-da;  sy=-1;}
+    else         {ya= a;  yya= da;  sy= 1;}
 
     /* (+++) The case 25 < abs(x) <= 1e8,    abs(y) <= 1e-7 */
     if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
@@ -355,14 +355,14 @@ tan(double x) {
       /* -cot */
       t2 = pz*(fi+gi)/(fi+pz);
       if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d)) { retval = (-sy*y); goto ret; }
-      t3 = (t2<ZERO) ? -t2 : t2;
+      t3 = (t2<0) ? -t2 : t2;
       t4 = gi*ua18.d+t3*ub18.d;
       if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
     else   {
       /* tan */
       t2 = pz*(gi+fi)/(gi-pz);
       if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d)) { retval = (sy*y); goto ret; }
-      t3 = (t2<ZERO) ? -t2 : t2;
+      t3 = (t2<0) ? -t2 : t2;
       t4 = fi*ua17.d+t3*ub17.d;
       if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
 
@@ -398,8 +398,8 @@ tan(double x) {
   /* Range reduction by algorithm iii */
   n = (__branred(x,&a,&da)) & 0x00000001;
   EADD(a,da,t1,t2)   a=t1;  da=t2;
-  if (a<ZERO)  {ya=-a;  yya=-da;  sy=MONE;}
-  else         {ya= a;  yya= da;  sy= ONE;}
+  if (a<0)  {ya=-a;  yya=-da;  sy=-1;}
+  else         {ya= a;  yya= da;  sy= 1;}
 
   /* (+++) The case 1e8 < abs(x) < 2**1024,    abs(y) <= 1e-7 */
   if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
@@ -463,14 +463,14 @@ tan(double x) {
     /* -cot */
     t2 = pz*(fi+gi)/(fi+pz);
     if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d)) { retval = (-sy*y); goto ret; }
-    t3 = (t2<ZERO) ? -t2 : t2;
+    t3 = (t2<0) ? -t2 : t2;
     t4 = gi*ua26.d+t3*ub26.d;
     if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
   else   {
     /* tan */
     t2 = pz*(gi+fi)/(gi-pz);
     if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d)) { retval = (sy*y); goto ret; }
-    t3 = (t2<ZERO) ? -t2 : t2;
+    t3 = (t2<0) ? -t2 : t2;
     t4 = fi*ua25.d+t3*ub25.d;
     if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
 
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpatan.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpatan.c
@@ -69,8 +69,8 @@ __mpatan(mp_no *x, mp_no *y, int p) {
 	{if (dx>__atan_xm[m].d) break;}
     }
     mpone.e    = mptwo.e    = mptwoim1.e = 1;
-    mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE;
-    mptwo.d[1] = TWO;
+    mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = 1;
+    mptwo.d[1] = 2;
 
 				 /* Reduce x m times */
     __mul(x,x,&mpsm,p);
@@ -92,7 +92,7 @@ __mpatan(mp_no *x, mp_no *y, int p) {
     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;
+      mptwoim1.d[1] -= 2;
       __dvd(&mpsm,&mptwoim1,&mpt1,p);
       __mul(&mpsm,&mpt,&mpt2,p);
       __sub(&mpt1,&mpt2,&mpt,p);
Index: glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.c
===================================================================
--- glibc-2.17-c758a686.orig/sysdeps/ieee754/dbl-64/mpsqrt.c
+++ glibc-2.17-c758a686/sysdeps/ieee754/dbl-64/mpsqrt.c
@@ -62,8 +62,8 @@ __mpsqrt(mp_no *x, mp_no *y, int p) {
   mp_no mpxn,mpz,mpu,mpt1,mpt2;
 
   /* Prepare multi-precision 1/2 and 3/2 */
-  mphalf.e  =0;  mphalf.d[0]  =ONE;  mphalf.d[1]  =HALFRAD;
-  mp3halfs.e=1;  mp3halfs.d[0]=ONE;  mp3halfs.d[1]=ONE;  mp3halfs.d[2]=HALFRAD;
+  mphalf.e  =0;  mphalf.d[0]  =1;  mphalf.d[1]  =HALFRAD;
+  mp3halfs.e=1;  mp3halfs.d[0]=1;  mp3halfs.d[1]=1;  mp3halfs.d[2]=HALFRAD;
 
   ey=EX/2;     __cpy(x,&mpxn,p);    mpxn.e -= (ey+ey);
   __mp_dbl(&mpxn,&dx,p);   dy=fastiroot(dx);    __dbl_mp(dy,&mpu,p);