[libmath-prime-util-perl] 12/55: New monty modinverse, nice primality speedup

[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 59d0b4ff22884dc918f9099860496815f3bde669
Author: Dana Jacobsen <dana at acm.org>
Date:   Wed Apr 30 11:37:02 2014 -0700

    New monty modinverse, nice primality speedup
---
 Changes     |  3 +++
 primality.c | 40 ++++++++++++++++------------------------
 2 files changed, 19 insertions(+), 24 deletions(-)

diff --git a/Changes b/Changes
index 57a5773..b78ae51 100644
--- a/Changes
+++ b/Changes
@@ -8,6 +8,9 @@ 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.
+
     - Small improvement to twin_prime_count_approx and nth_twin_prime_approx.
 
     - Better AKS testing in xt/primality-aks.pl.
diff --git a/primality.c b/primality.c
index 38d3e7f..d6f3962 100644
--- a/primality.c
+++ b/primality.c
@@ -18,7 +18,7 @@ static const UV mr_bases_const2[1] = {2};
   Code inside USE_MONT_PRIMALITY is Montgomery math and efficient M-R from
   Wojciech Izykowski.  See:  https://github.com/wizykowski/miller-rabin
 
-Copyright (c) 2013, Wojciech Izykowski
+Copyright (c) 2013-2014, Wojciech Izykowski
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
@@ -68,28 +68,20 @@ static INLINE UV mont_powmod64(uint64_t a, uint64_t k, uint64_t one, uint64_t n,
   }
   return t;
 }
+/* Returns -a^-1 mod 2^64.  From B. Arazi "On Primality Testing Using Purely
+ * Divisionless Operations", Computer Journal (1994) 37 (3): 219-222, Proc 5 */
 static INLINE uint64_t modular_inverse64(const uint64_t a)
 {
-  uint64_t u,x,w,z,q;
-
-  x = 1; z = a;
-
-  q = (-z)/z + 1;  /* = 2^64 / z                 */
-  u = - q;         /* = -q * x                   */
-  w = - q * z;     /* = b - q * z = 2^64 - q * z */
-
-  /* after first iteration all variables are 64-bit */
-
-  while (w) {
-    if (w < z) {
-      q = u; u = x; x = q; /* swap(u, x) */
-      q = w; w = z; z = q; /* swap(w, z) */
+  uint64_t S = 1, J = 0;
+  int i;
+  for (i = 0; i < 64; i++) {
+    if (S & 1) {
+      J |= (1ULL << i);
+      S += a;
     }
-    q = w / z;
-    u -= q * x;
-    w -= q * z;
+    S >>= 1;
   }
-  return x;
+  return J;
 }
 static INLINE uint64_t compute_modn64(const uint64_t n)
 {
@@ -118,7 +110,7 @@ static INLINE uint64_t compute_2_65_mod_n(const uint64_t n, const uint64_t modn)
 static int monty_mr64(const uint64_t n, const UV* bases, int cnt)
 {
   int i, j, t;
-  const uint64_t npi = modular_inverse64(-((int64_t)n));
+  const uint64_t npi = modular_inverse64(n);
   const uint64_t r = compute_modn64(n);
   uint64_t u = n - 1;
   const uint64_t nr = n - r;
@@ -239,7 +231,7 @@ int _XS_BPSW(UV const n)
     return    _XS_miller_rabin(n, mr_bases_const2, 1)
            && _XS_is_almost_extra_strong_lucas_pseudoprime(n,1);
   } else {
-    const uint64_t npi = modular_inverse64(-((int64_t)n));
+    const uint64_t npi = modular_inverse64(n);
     const uint64_t montr = compute_modn64(n);
     const uint64_t mont2 = compute_2_65_mod_n(n, montr);
     uint64_t u = n-1;
@@ -445,7 +437,7 @@ int _XS_is_lucas_pseudoprime(UV n, int strength)
 
 #if USE_MONT_PRIMALITY
   if (n > UVCONST(4294967295)) {
-    const uint64_t npi = modular_inverse64(-((int64_t)n));
+    const uint64_t npi = modular_inverse64(n);
     const uint64_t mont1 = compute_modn64(n);
     const uint64_t mont2 = compute_2_65_mod_n(n, mont1);
     const uint64_t montP = (P == 1) ? mont1
@@ -604,7 +596,7 @@ int _XS_is_almost_extra_strong_lucas_pseudoprime(UV n, UV increment)
 
 #if USE_MONT_PRIMALITY
   if (n > UVCONST(4294967295)) {
-    const uint64_t npi = modular_inverse64(-((int64_t)n));
+    const uint64_t npi = modular_inverse64(n);
     const uint64_t montr = compute_modn64(n);
     const uint64_t mont2 = compute_2_65_mod_n(n, montr);
     const uint64_t montP = compute_a_times_2_64_mod_n(P, n, montr);
@@ -689,7 +681,7 @@ int _XS_is_frobenius_underwood_pseudoprime(UV n)
 
 #if USE_MONT_PRIMALITY
   if (n > UVCONST(4294967295)) {
-    const uint64_t npi = modular_inverse64(-((int64_t)n));
+    const uint64_t npi = modular_inverse64(n);
     const uint64_t mont1 = compute_modn64(n);
     const uint64_t mont2 = compute_2_65_mod_n(n, mont1);
     const uint64_t mont5 = compute_a_times_2_64_mod_n(5, n, mont1);

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