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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits