This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[PATCH] New function to generate mp_no from powers of two


Hi,

The __mpexp_twomm1 array (and a lot of unnecessary computation around
it) is required in mpexp because there is currently no way to generate
mp_no numbers that are powers of two.  These are actually very simple
numbers since they just have the exponent, sign and the first mantissa
digit set and hence are a good special case to have (although I
haven't done a code search to see yet if it's useful in any other
functions).  Attached patch implements an inline function __pow_mp
(inline since it is not a very big or complicated function) and
implements it in its first user, which is the __mpexp function.  I
have verified that this does not break any test cases on x86_64.

Performance:

With the patch:

Total:26237847075, Fastest:255700, Slowest:893136, Avg:262378.470750

Without the patch:

Total:26245174617, Fastest:256384, Slowest:895106, Avg:262451.746170

That's a miniscule improvement in performance (0.2% in the fastest
case and barely 0.02% in the average case) but the real advantage is
that the code is cleaner.  The numbers are on top of the fast
multiplication patch[1], which is why they're so low.

OK to commit?

Siddhesh

[1] http://sourceware.org/ml/libc-alpha/2013-01/msg00614.html


	* 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.


diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h
index debb3b2..70a3614 100644
--- a/sysdeps/ieee754/dbl-64/mpa.h
+++ b/sysdeps/ieee754/dbl-64/mpa.h
@@ -123,3 +123,25 @@ 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;
+
+  EY = pow / 24;
+  rem = pow - EY * 24;
+  EY++;
+
+  if (rem < 0)
+    {
+      EY--;
+      rem += 24;
+    }
+  Y[0] = ONE;
+  Y[1] = 1 << rem;
+
+  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 811fb46..8d288ff 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.  */


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]