[Gmp-commit] /var/hg/gmp: 2 new changesets

[mercurial at gmplib.org]

details:   /var/hg/gmp/rev/875b15a6160d
changeset: 16797:875b15a6160d
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 31 18:31:26 2015 +0200
description:
mpn/generic/powlo.c: Use mpn_sqrlo.

details:   /var/hg/gmp/rev/d8428a0f1b2c
changeset: 16798:d8428a0f1b2c
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 31 18:31:59 2015 +0200
description:
ChangeLog

diffstat:

 ChangeLog           |   2 ++
 mpn/generic/powlo.c |  52 ++++++++++++++++++++++++++--------------------------
 2 files changed, 28 insertions(+), 26 deletions(-)

diffs (137 lines):

diff -r 4d02fcf88613 -r d8428a0f1b2c ChangeLog
--- a/ChangeLog	Mon Aug 31 10:40:37 2015 +0200
+++ b/ChangeLog	Mon Aug 31 18:31:59 2015 +0200
@@ -7,6 +7,8 @@
 	* tests/mpq/t-cmp_z.c: New file to test mpq_cmp_z (from t-cmp.c).
 	* tests/mpq/Makefile.am (check_PROGRAMS): Add t-cmp_z.
 
+	* mpn/generic/powlo.c: Use mpn_sqrlo.
+
 2015-08-29  Torbjörn Granlund  <torbjorng at google.com>
 
 	* mpn/x86_64/fat/fat.c (fake_cpuid_table): Add CPU aliases.
diff -r 4d02fcf88613 -r d8428a0f1b2c mpn/generic/powlo.c
--- a/mpn/generic/powlo.c	Mon Aug 31 10:40:37 2015 +0200
+++ b/mpn/generic/powlo.c	Mon Aug 31 18:31:59 2015 +0200
@@ -1,6 +1,6 @@
 /* mpn_powlo -- Compute R = U^E mod B^n, where B is the limb base.
 
-Copyright 2007-2009, 2012 Free Software Foundation, Inc.
+Copyright 2007-2009, 2012, 2015 Free Software Foundation, Inc.
 
 This file is part of the GNU MP Library.
 
@@ -66,7 +66,8 @@
 {
   int k;
   static mp_bitcnt_t x[] = {1,7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0};
-  for (k = 0; eb > x[k]; k++)
+  ASSERT (eb > 1);
+  for (k = 1; eb > x[k]; ++k)
     ;
   return k;
 }
@@ -74,6 +75,8 @@
 /* rp[n-1..0] = bp[n-1..0] ^ ep[en-1..0] mod B^n, B is the limb base.
    Requires that ep[en-1] is non-zero.
    Uses scratch space tp[3n-1..0], i.e., 3n words.  */
+/* We only use n words in the scratch space, we should pass tp + n to
+   mullo/sqrlo as a temporary area, it is needed. */
 void
 mpn_powlo (mp_ptr rp, mp_srcptr bp,
 	   mp_srcptr ep, mp_size_t en,
@@ -84,7 +87,6 @@
   int windowsize, this_windowsize;
   mp_limb_t expbits;
   mp_limb_t *pp, *this_pp, *last_pp;
-  mp_limb_t *b2p;
   long i;
   TMP_DECL;
 
@@ -95,47 +97,42 @@
   MPN_SIZEINBASE_2EXP(ebi, ep, en, 1);
 
   windowsize = win_size (ebi);
+  ASSERT (windowsize < ebi);
 
-  pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)) + n); /* + n is for mullo ign part */
+  pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)));
 
   this_pp = pp;
 
   MPN_COPY (this_pp, bp, n);
 
-  b2p = tp + 2*n;
-
-  /* Store b^2 in b2.  */
-  mpn_sqr (tp, bp, n);	/* FIXME: Use "mpn_sqrlo" */
-  MPN_COPY (b2p, tp, n);
+  /* Store b^2 in tp.  */
+  mpn_sqrlo (tp, bp, n);
 
   /* Precompute odd powers of b and put them in the temporary area at pp.  */
   for (i = (1 << (windowsize - 1)) - 1; i > 0; i--)
     {
       last_pp = this_pp;
       this_pp += n;
-      mpn_mullo_n (this_pp, last_pp, b2p, n);
+      mpn_mullo_n (this_pp, last_pp, tp, n);
     }
 
   expbits = getbits (ep, ebi, windowsize);
-  if (ebi < windowsize)
-    ebi = 0;
-  else
-    ebi -= windowsize;
 
+  /* FIXME: for even expbits, we can init with a mullo. */
   count_trailing_zeros (cnt, expbits);
+  ebi -= windowsize;
   ebi += cnt;
   expbits >>= cnt;
 
   MPN_COPY (rp, pp + n * (expbits >> 1), n);
 
-  while (ebi != 0)
+  do
     {
       while (getbit (ep, ebi) == 0)
 	{
-	  mpn_sqr (tp, rp, n);	/* FIXME: Use "mpn_sqrlo" */
+	  mpn_sqrlo (tp, rp, n);
 	  MPN_COPY (rp, tp, n);
-	  ebi--;
-	  if (ebi == 0)
+	  if (--ebi == 0)
 	    goto done;
 	}
 
@@ -157,17 +154,20 @@
       ebi += cnt;
       expbits >>= cnt;
 
-      do
+      while (this_windowsize > 1)
 	{
-	  mpn_sqr (tp, rp, n);
-	  MPN_COPY (rp, tp, n);
-	  this_windowsize--;
+	  mpn_sqrlo (tp, rp, n);
+	  mpn_sqrlo (rp, tp, n);
+	  this_windowsize -= 2;
 	}
-      while (this_windowsize != 0);
 
-      mpn_mullo_n (tp, rp, pp + n * (expbits >> 1), n);
-      MPN_COPY (rp, tp, n);
-    }
+      if (this_windowsize != 0)
+	mpn_sqrlo (tp, rp, n);
+      else
+	MPN_COPY (tp, rp, n);
+      
+      mpn_mullo_n (rp, tp, pp + n * (expbits >> 1), n);
+    } while (ebi != 0);
 
  done:
   TMP_FREE;