[libmath-prime-util-perl] 14/55: Improved code for mont_prod64

[Partha P. Mukherjee]

This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.41
in repository libmath-prime-util-perl.

commit 9733e53ae29ad9b43188e38781ae9750c42b6998
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu May 1 14:24:44 2014 -0700

    Improved code for mont_prod64
---
 Changes     |  3 ++-
 primality.c | 32 +++++++++++++++++++++++++++-----
 2 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/Changes b/Changes
index b78ae51..b9285c0 100644
--- a/Changes
+++ b/Changes
@@ -9,7 +9,8 @@ Revision history for Perl module Math::Prime::Util
     [FUNCTIONALITY AND PERFORMANCE]
 
     - New modular inverse from W. Izykowski, from Arazi (1994).  This is a
-      substantial (10-20%) speedup for 64-bit primality testing.
+      substantial (20-40%) speedup for primality testing in range 2^32 to 2^64.
+      This affects functions like next_prime, prev_prime, etc.
 
     - Small improvement to twin_prime_count_approx and nth_twin_prime_approx.
 
diff --git a/primality.c b/primality.c
index d6f3962..02b0d40 100644
--- a/primality.c
+++ b/primality.c
@@ -46,16 +46,14 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 static INLINE uint64_t mont_prod64(uint64_t a, uint64_t b, uint64_t n, uint64_t npi)
 {
   uint64_t t_hi, t_lo, m, mn_hi, mn_lo, u;
-  int carry;
   /* t_hi * 2^64 + t_lo = a*b */
   asm("mulq %3" : "=a"(t_lo), "=d"(t_hi) : "a"(a), "rm"(b));
   m = t_lo * npi;
   /* mn_hi * 2^64 + mn_lo = m*n */
   asm("mulq %3" : "=a"(mn_lo), "=d"(mn_hi) : "a"(m), "rm"(n));
-  carry = t_lo + mn_lo < t_lo ? 1 : 0;
-  u = t_hi + mn_hi + carry;
-  if (u < t_hi) return u-n;
-  return u >= n ? u-n : u;
+  u = t_hi + mn_hi;
+  if (t_lo + mn_lo < t_lo) u++;
+  return (u < t_hi || u >= n)  ?  u-n  :  u;
 }
 #define mont_square64(a, n, npi)  mont_prod64(a, a, n, npi)
 static INLINE UV mont_powmod64(uint64_t a, uint64_t k, uint64_t one, uint64_t n, uint64_t npi)
@@ -72,6 +70,7 @@ static INLINE UV mont_powmod64(uint64_t a, uint64_t k, uint64_t one, uint64_t n,
  * Divisionless Operations", Computer Journal (1994) 37 (3): 219-222, Proc 5 */
 static INLINE uint64_t modular_inverse64(const uint64_t a)
 {
+#if 0
   uint64_t S = 1, J = 0;
   int i;
   for (i = 0; i < 64; i++) {
@@ -81,6 +80,29 @@ static INLINE uint64_t modular_inverse64(const uint64_t a)
     }
     S >>= 1;
   }
+#else
+  /* 2 at a time. */
+  uint64_t S = 1, J = 0;
+  const uint64_t a1p = (a>>1)+1, a2 = a>>2, a2p = a2+1;
+  int i;
+  if (a & 2) {
+    const uint64_t tab[4] = { 0,  a2p,  a1p,  a1p+a2p };
+    for (i = 0; i < 32; i++) {
+      uint64_t S2 = S >> 2;
+      int index = S & 3;
+      J |= (uint64_t)index << (2 * i);
+      S = S2 + tab[index];
+    }
+  } else {
+    const uint64_t tab[4] = { 0,  a1p+a2,  a1p,  a2p };
+    for (i = 0; i < 32; i++) {
+      uint64_t S2 = S >> 2;
+      int index = S & 3;
+      J |= (uint64_t)((4-index)&3) << (2 * i);
+      S = S2 + tab[index];
+    }
+  }
+#endif
   return J;
 }
 static INLINE uint64_t compute_modn64(const uint64_t n)

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git