Simplify calculation of 2^-m in __mpexp

3 files changed, 39 insertions, 29 deletions

diff --git a/ChangeLog b/ChangeLog
index 1d17c4d428..52a95426d1 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,10 @@
 2013-01-18  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
+	* sysdeps/ieee754/dbl-64/mpa.h (__pow_mp): New function to get an
+	mp_no from a power of two.
+	* sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): Remove
+	__mpexp_twomm1.  Use __pow_mp.
+
 	* sysdeps/ieee754/dbl-64/mpexp.c (__mpexp): Remove unnecessary
 	multiplication.
 
diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index debb3b2a96..06343d46d1 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -123,3 +123,33 @@ extern void __mpsqrt (mp_no *, mp_no *, 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] = ONE;
+  Y[1] = 1 << rem;
+
+  /* Everything else is ZERO.  */
+  for (i = 2; i <= p; i++)
+    Y[i] = ZERO;
+}
diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c
index 811fb465a1..8d288ff9a1 100644
--- a/sysdeps/ieee754/dbl-64/mpexp.c
+++ b/sysdeps/ieee754/dbl-64/mpexp.c
@@ -43,7 +43,7 @@ SECTION
 __mpexp (mp_no *x, mp_no *y, int p)
 {
   int i, j, k, m, m1, m2, n;
-  double a, b;
+  double b;
   static const int np[33] =
     {
       0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6,
@@ -61,19 +61,6 @@ __mpexp (mp_no *x, mp_no *y, int p)
       70, 73, 76, 78,
       81
     };
-  /* Stored values for 2^-m, where values of m are defined in M1P above.   */
-  static const double __mpexp_twomm1[33] =
-    {
-      0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0,
-      0x1.0p-17, 0x1.0p-23, 0x1.0p-23, 0x1.0p-28,
-      0x1.0p-27, 0x1.0p-38, 0x1.0p-42, 0x1.0p-39,
-      0x1.0p-43, 0x1.0p-47, 0x1.0p-43, 0x1.0p-47,
-      0x1.0p-50, 0x1.0p-54, 0x1.0p-57, 0x1.0p-60,
-      0x1.0p-64, 0x1.0p-67, 0x1.0p-71, 0x1.0p-74,
-      0x1.0p-68, 0x1.0p-71, 0x1.0p-74, 0x1.0p-77,
-      0x1.0p-70, 0x1.0p-73, 0x1.0p-76, 0x1.0p-78,
-      0x1.0p-81
-    };
   static const int m1np[7][18] =
     {
       {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
@@ -98,18 +85,10 @@ __mpexp (mp_no *x, mp_no *y, int p)
   /* Choose m,n and compute a=2**(-m).  */
   n = np[p];
   m1 = m1p[p];
-  a = __mpexp_twomm1[p];
-  for (i = 0; i < EX; i++)
-    a *= RADIXI;
-  for (; i > EX; i--)
-    a *= RADIX;
   b = X[1];
   m2 = 24 * EX;
   for (; b < HALFRAD; m2--)
-    {
-      a *= TWO;
-      b *= TWO;
-    }
+    b *= TWO;
   if (b == HALFRAD)
     {
       for (i = 2; i <= p; i++)
@@ -118,10 +97,7 @@ __mpexp (mp_no *x, mp_no *y, int p)
 	    break;
 	}
       if (i == p + 1)
-	{
-	  m2--;
-	  a *= TWO;
-	}
+	m2--;
     }
 
   m = m1 + m2;
@@ -134,14 +110,13 @@ __mpexp (mp_no *x, mp_no *y, int p)
 	 than 2^-55.  */
       assert (p < 18);
       m = 0;
-      a = ONE;
       for (i = n - 1; i > 0; i--, n--)
 	if (m1np[i][p] + m2 > 0)
 	  break;
     }
 
   /* Compute s=x*2**(-m). Put result in mps.  */
-  __dbl_mp (a, &mpt1, p);
+  __pow_mp (-m, &mpt1, p);
   __mul (x, &mpt1, &mps, p);
 
   /* Evaluate the polynomial. Put result in mpt2.  */