Clean up add_magnitudes and sub_magnitudes

diff --git a/ChangeLog b/ChangeLog
index b46fb85..9b224a4 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,9 @@
+2013-02-13  Siddhesh Poyarekar  <siddhesh@redhat.com>
+
+       * sysdeps/ieee754/dbl-64/mpa.c (add_magnitudes): Use ZK to
+       minimize writes to Z.
+       (sub_magnitudes): Simplify code a bit.
+
 2013-02-12  Roland McGrath  <roland@hack.frob.com>
 
        * include/libc-symbols.h (stub_warning): Remove "warning: " prefix

diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c
index ede8ed1..005aa0c 100644 (file)
--- a/sysdeps/ieee754/dbl-64/mpa.c
+++ b/sysdeps/ieee754/dbl-64/mpa.c
@@ -370,6 +370,7 @@ SECTION
 add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   int i, j, k;
+  double zk;
 
   EZ = EX;
 
@@ -377,45 +378,54 @@ add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
   j = p + EY - EX;
   k = p + 1;
 
-  if (j < 1)
+  if (__glibc_unlikely (j < 1))
     {
       __cpy (x, z, p);
       return;
     }
-  else
-    Z[k] = ZERO;
+
+  zk = ZERO;
 
   for (; j > 0; i--, j--)
     {
-      Z[k] += X[i] + Y[j];
-      if (Z[k] >= RADIX)
+      zk += X[i] + Y[j];
+      if (zk >= RADIX)
        {
-         Z[k] -= RADIX;
-         Z[--k] = ONE;
+         Z[k--] = zk - RADIX;
+         zk = ONE;
        }
       else
-       Z[--k] = ZERO;
+        {
+         Z[k--] = zk;
+         zk = ZERO;
+       }
     }
 
   for (; i > 0; i--)
     {
-      Z[k] += X[i];
-      if (Z[k] >= RADIX)
+      zk += X[i];
+      if (zk >= RADIX)
        {
-         Z[k] -= RADIX;
-         Z[--k] = ONE;
+         Z[k--] = zk - RADIX;
+         zk = ONE;
        }
       else
-       Z[--k] = ZERO;
+        {
+         Z[k--] = zk;
+         zk = ZERO;
+       }
     }
 
-  if (Z[1] == ZERO)
+  if (zk == ZERO)
     {
       for (i = 1; i <= p; i++)
        Z[i] = Z[i + 1];
     }
   else
-    EZ += ONE;
+    {
+      Z[1] = zk;
+      EZ += ONE;
+    }
 }
 
 /* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
@@ -427,65 +437,63 @@ SECTION
 sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
 {
   int i, j, k;
+  double zk;
 
   EZ = EX;
+  i = p;
+  j = p + EY - EX;
+  k = p;
 
-  if (EX == EY)
+  /* Y is too small compared to X, copy X over to the result.  */
+  if (__glibc_unlikely (j < 1))
     {
-      i = j = k = p;
-      Z[k] = Z[k + 1] = ZERO;
+      __cpy (x, z, p);
+      return;
     }
-  else
+
+  /* The relevant least significant digit in Y is non-zero, so we factor it in
+     to enhance accuracy.  */
+  if (j < p && Y[j + 1] > ZERO)
     {
-      j = EX - EY;
-      if (j > p)
-       {
-         __cpy (x, z, p);
-         return;
-       }
-      else
-       {
-         i = p;
-         j = p + 1 - j;
-         k = p;
-         if (Y[j] > ZERO)
-           {
-             Z[k + 1] = RADIX - Y[j--];
-             Z[k] = MONE;
-           }
-         else
-           {
-             Z[k + 1] = ZERO;
-             Z[k] = ZERO;
-             j--;
-           }
-       }
+      Z[k + 1] = RADIX - Y[j + 1];
+      zk = MONE;
     }
+  else
+    zk = Z[k + 1] = ZERO;
 
+  /* Subtract and borrow.  */
   for (; j > 0; i--, j--)
     {
-      Z[k] += (X[i] - Y[j]);
-      if (Z[k] < ZERO)
+      zk += (X[i] - Y[j]);
+      if (zk < ZERO)
        {
-         Z[k] += RADIX;
-         Z[--k] = MONE;
+         Z[k--] = zk + RADIX;
+         zk = MONE;
        }
       else
-       Z[--k] = ZERO;
+        {
+         Z[k--] = zk;
+         zk = ZERO;
+       }
     }
 
+  /* We're done with digits from Y, so it's just digits in X.  */
   for (; i > 0; i--)
     {
-      Z[k] += X[i];
-      if (Z[k] < ZERO)
+      zk += X[i];
+      if (zk < ZERO)
        {
-         Z[k] += RADIX;
-         Z[--k] = MONE;
+         Z[k--] = zk + RADIX;
+         zk = MONE;
        }
       else
-       Z[--k] = ZERO;
+        {
+         Z[k--] = zk;
+         zk = ZERO;
+       }
     }
 
+  /* Normalize.  */
   for (i = 1; Z[i] == ZERO; i++);
   EZ = EZ - i + 1;
   for (k = 1; i <= p + 1;)