[Gmp-commit] /var/hg/gmp-proj/mini-gmp: 3 new changesets

mercurial at gmplib.org mercurial at gmplib.org
Wed Jan 4 21:26:36 CET 2012


details:   /var/hg/gmp-proj/mini-gmp/rev/4bde2db4ad2f
changeset: 43:4bde2db4ad2f
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Wed Jan 04 19:53:34 2012 +0100
description:
Implemented mpz_invert and mpz_powm.

details:   /var/hg/gmp-proj/mini-gmp/rev/f776169353f3
changeset: 44:f776169353f3
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Wed Jan 04 19:56:16 2012 +0100
description:
Tests for powm, and logops (forgotten in earlier commit).

details:   /var/hg/gmp-proj/mini-gmp/rev/b753004ece2c
changeset: 45:b753004ece2c
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Wed Jan 04 21:25:37 2012 +0100
description:
Improved preinv division, and use it in mpz_powm.
struct gmp_div_qr_1_inverse: generalized and renamed to gmp_div_inverse.
mpn_div_qr_2_invert, mpn_div_qr_invert: New functions.
mpn_div_qr_2_preinv: New function.
mpn_div_qr_preinv: New function. For dn > 2, requires that caller
normalizes d. Still shifts n as needed.
mpn_div_qr: Rewritten to use mpn_div_qr_preinv.
mpz_powm: Use mpn_div_qr_preinv.

diffstat:

 .hgignore           |    1 +
 mini-gmp.c          |  310 +++++++++++++++++++++++++++++++++++++++++++--------
 mini-gmp.h          |    3 +
 tests/Makefile      |    2 +-
 tests/hex-random.c  |  118 +++++++++++++------
 tests/hex-random.h  |   11 +-
 tests/mini-random.c |   18 +-
 tests/mini-random.h |    4 +-
 tests/t-div.c       |    2 +-
 tests/t-logops.c    |   75 ++++++++++++
 tests/t-powm.c      |   54 +++++++++
 11 files changed, 492 insertions(+), 106 deletions(-)

diffs (truncated from 856 to 300 lines):

diff -r 0a5e174f762f -r b753004ece2c .hgignore
--- a/.hgignore	Tue Jan 03 22:31:57 2012 +0100
+++ b/.hgignore	Wed Jan 04 21:25:37 2012 +0100
@@ -10,5 +10,6 @@
 ^tests/t-div$
 ^tests/t-double$
 ^tests/t-gcd$
+^tests/t-powm$
 ^tests/t-logops$
 ^tests/t-str$
diff -r 0a5e174f762f -r b753004ece2c mini-gmp.c
--- a/mini-gmp.c	Tue Jan 03 22:31:57 2012 +0100
+++ b/mini-gmp.c	Wed Jan 04 21:25:37 2012 +0100
@@ -674,31 +674,83 @@
   return m;
 }
 
-struct gmp_div_qr_1_inverse
+struct gmp_div_inverse
 {
   /* Normalization shift count. */
   unsigned shift;
-  /* Normalized divisor */
-  mp_limb_t d;
-  /* Inverse */
+  /* Normalized divisor (d0 unused for mpn_div_qr_1) */
+  mp_limb_t d1, d0;
+  /* Inverse, for 2/1 or 3/2. */
   mp_limb_t di;
 };
 
 static void
-mpn_div_qr_1_invert (struct gmp_div_qr_1_inverse *inv, mp_limb_t d)
+mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
 {
   unsigned shift;
+
+  assert (d > 0);
   gmp_clz (shift, d);
   inv->shift = shift;
-  inv->d = d << shift;
-  inv->di = mpn_invert_limb (inv->d);
+  inv->d1 = d << shift;
+  inv->di = mpn_invert_limb (inv->d1);
+}
+
+static void
+mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
+		     mp_limb_t d1, mp_limb_t d0)
+{
+  unsigned shift;
+
+  assert (d1 > 0);
+  gmp_clz (shift, d1);
+  inv->shift = shift;
+  if (shift > 0)
+    {
+      d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
+      d0 <<= shift;
+    }
+  inv->d1 = d1;
+  inv->d0 = d0;
+  inv->di = mpn_invert_3by2 (d1, d0);
+}
+
+static void
+mpn_div_qr_invert (struct gmp_div_inverse *inv,
+		   mp_srcptr dp, mp_size_t dn)
+{
+  assert (dn > 0);
+  
+  if (dn == 1)
+    mpn_div_qr_1_invert (inv, dp[0]);
+  else if (dn == 2)
+    mpn_div_qr_2_invert (inv, dp[1], dp[0]);
+  else
+    {
+      unsigned shift;
+      mp_limb_t d1, d0;
+
+      d1 = dp[dn-1];
+      d0 = dp[dn-2];
+      assert (d1 > 0);
+      gmp_clz (shift, d1);
+      inv->shift = shift;
+      if (shift > 0)
+	{
+	  d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
+	  d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
+	}
+      inv->d1 = d1;
+      inv->d0 = d0;
+      inv->di = mpn_invert_3by2 (d1, d0);
+    }
 }
 
 /* Not matching current public gmp interface, rather corresponding to
    the sbpi1_div_* functions. */
 static mp_limb_t
 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
-		     const struct gmp_div_qr_1_inverse *inv)
+		     const struct gmp_div_inverse *inv)
 {
   mp_limb_t d, di;
   mp_limb_t r;
@@ -713,8 +765,7 @@
   else
     r = 0;
 
-
-  d = inv->d;
+  d = inv->d1;
   di = inv->di;
   while (nn-- > 0)
     {
@@ -735,40 +786,42 @@
 {
   assert (d > 0);
 
+  /* Special case for powers of two. */
   if (d > 1 && (d & (d-1)) == 0)
     {
-      /* Power of two */
       unsigned shift;
       mp_limb_t r = np[0] & (d-1);
       gmp_ctz (shift, d);
-      mpn_rshift (qp, np, nn, shift);
+      if (qp)
+	mpn_rshift (qp, np, nn, shift);
 
       return r;
     }
   else
     {
-      struct gmp_div_qr_1_inverse inv;
+      struct gmp_div_inverse inv;
       mpn_div_qr_1_invert (&inv, d);
       return mpn_div_qr_1_preinv (qp, np, nn, &inv);
     }
 }
 
 static void
-mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
-	      mp_limb_t d1, mp_limb_t d0)
+mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
+		     const struct gmp_div_inverse *inv)
 {
   unsigned shift;
   mp_size_t i;
-  mp_limb_t di, r1, r0;
+  mp_limb_t d1, d0, di, r1, r0;
   mp_ptr tp;
 
   assert (nn >= 2);
-
-  gmp_clz (shift, d1);
+  shift = inv->shift;
+  d1 = inv->d1;
+  d0 = inv->d0;
+  di = inv->di;
+  
   if (shift > 0)
     {
-      d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
-      d0 <<= shift;
       tp = gmp_xalloc_limbs (nn);
       r1 = mpn_lshift (tp, np, nn, shift);
       np = tp;
@@ -776,9 +829,6 @@
   else
     r1 = 0;
 
-  di = mpn_invert_3by2 (d1, d0);
-
-  /* np += nn - 1; */
   r0 = np[nn - 1];
 
   for (i = nn - 2; i >= 0; i--)
@@ -804,6 +854,19 @@
   rp[0] = r0;
 }
 
+#if 0
+static void
+mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
+	      mp_limb_t d1, mp_limb_t d0)
+{
+  struct gmp_div_inverse inv;
+  assert (nn >= 2);
+
+  mpn_div_qr_2_invert (&inv, d1, d0);
+  mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
+}
+#endif
+
 static void
 mpn_div_qr_pi1 (mp_ptr qp,
 		mp_ptr np, mp_size_t nn, mp_limb_t n1,
@@ -866,48 +929,61 @@
 }
 
 static void
-mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
+mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
+		   mp_srcptr dp, mp_size_t dn,
+		   const struct gmp_div_inverse *inv)
 {
   assert (dn > 0);
   assert (nn >= dn);
 
   if (dn == 1)
-    np[0] = mpn_div_qr_1 (qp, np, nn, dp[0]);
+    np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
   else if (dn == 2)
-    mpn_div_qr_2 (qp, np, np, nn, dp[1], dp[0]);
+    mpn_div_qr_2_preinv (qp, np, np, nn, inv);
   else
     {
-      mp_ptr tp = NULL;
-      mp_limb_t d1, d0, di;
       mp_limb_t nh;
       unsigned shift;
 
-      d1 = dp[dn - 1];
-
-      gmp_clz (shift, d1);
+      assert (dp[dn-1] & GMP_LIMB_HIGHBIT);
+
+      shift = inv->shift;
       if (shift > 0)
-	{
-	  tp = gmp_xalloc_limbs (dn);
-	  gmp_assert_nocarry (mpn_lshift (tp, dp, dn, shift));
-	  nh = mpn_lshift (np, np, nn, shift);
-	  dp = tp;
-	  d1 = dp[dn - 1];
-	}
+	nh = mpn_lshift (np, np, nn, shift);
       else
 	nh = 0;
 
-      d0 = dp[dn-2];
-      di = mpn_invert_3by2 (d1, d0);
-      mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, di);
+      assert (inv->d1 == dp[dn-1]);
+      assert (inv->d0 == dp[dn-2]);
+      
+      mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
 
       if (shift > 0)
-	{
-	  free (tp);
-	  gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
-	}
+	gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
     }
 }
 
+static void
+mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
+{
+  struct gmp_div_inverse inv;
+  mp_ptr tp = NULL;
+
+  assert (dn > 0);
+  assert (nn >= dn);
+
+  mpn_div_qr_invert (&inv, dp, dn);
+  if (dn > 2 && inv.shift > 0)
+    {
+      tp = gmp_xalloc_limbs (dn);
+      gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
+      dp = tp;
+    }
+  mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
+  if (tp)
+    free (tp);
+}
+
 static unsigned
 mpn_base_power_of_two_p (unsigned b)
 {
@@ -992,7 +1068,7 @@
    the end. */
 static size_t
 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
-	      const struct gmp_div_qr_1_inverse *binv)
+		  const struct gmp_div_inverse *binv)
 {
   mp_size_t i;
   for (i = 0; w > 0; i++)
@@ -1002,7 +1078,7 @@
       h = w >> (GMP_LIMB_BITS - binv->shift);
       l = w << binv->shift;
 
-      gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d, binv->di);
+      gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
       assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
       r >>= binv->shift;
 


More information about the gmp-commit mailing list