From 5e058bdad3be9d97cf173f50ce6eb432c4f63519 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Sun, 28 Feb 2021 03:32:30 -0600 Subject: [PATCH 1/9] Use trial division and pohlard's P-1 algorithm to speed up ismersenne when false --- src/Primes.jl | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/src/Primes.jl b/src/Primes.jl index 2a0e075..6dc980a 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -459,9 +459,70 @@ function ismersenneprime(M::Integer; check::Bool = true) throw(ArgumentError("The argument given is not a valid Mersenne Number (`M = 2^p - 1`).")) end M < 7 && return M == 3 + mersennehassmallfactor(M) && return false + p_minus1(M) && return false return ll_primecheck(M) end +# Returns whether M=2^d-1 has a small factor +# For each factor, checks if 2^d = 1 \pmod factor. +# Uses the fact that all factors must be of the form +# 2kd+1 for integer k +function has_small_factor_mersenne(M::Integer) + d = ndigits(M, base=2) + #limit = round(BigInt, 2^(22.85*d^(1/16))/(2*d)) + limit = round(BigInt, 7.3e-6 * d^2.32) + step = 2 * d + wheel = Iterators.cycle(d % 4 == 1 ? + (false,false,true,true) : + (true,false,false,true)) + for (factor, w) in zip(1+step : step : step*limit, wheel) + !w && continue + factor & 7 in (3,5) && continue + !isprime(factor, 4) && continue + if powermod(2, d, factor) == 1 + return true + end + end + return false +end + +# Calculates n % 2^d-1 where M=2**d-1 +function mod_mersenne(n, d, M) + while n > M + n = (n & M) + (n >> d) + end + return n != M ? n : 0 +end + +# Calculates base^exp % 2^d-1 where M=2**d-1 +function power_mod_mersenne(base, exp, d, M) + ans = 1 + while exp > 0 + if exp & 1 == 1 + ans = mod_mersenne(ans * base, d, M) + end + exp >>= 1 + base = mod_mersenne(base * base, d, M) + end + return ans +end + +# Returns if M=2^prime-1 have a factor 2*k*prime+1 +# such that the largest prime factor of k is less than limit""" +function p_minus1_mersenne(M) + d = ndigits(M, base=2) + B1 = d + log_B1 = log(B1) + P = prod(p^round(BigInt, log_B1/log(p)) for p in primes(B1)) + P = power_mod_mersenne(3, 2*P*d, d, M) + g = gcd(M, P-1) + 1 < g < M && return true + g == M && return false + return false + #return p_minus1_stage_2(prime, mersenne, P, p, 10*B1, primes) +end + """ isrieselprime(k::Integer, Q::Integer) -> Bool From 8b6072d34d7f9e8931204f9762e08cff9e37a9b3 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Sun, 28 Feb 2021 11:16:28 -0600 Subject: [PATCH 2/9] fix typo --- src/Primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Primes.jl b/src/Primes.jl index 6dc980a..36fbdda 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -459,7 +459,7 @@ function ismersenneprime(M::Integer; check::Bool = true) throw(ArgumentError("The argument given is not a valid Mersenne Number (`M = 2^p - 1`).")) end M < 7 && return M == 3 - mersennehassmallfactor(M) && return false + has_small_factor_mersenne(M) && return false p_minus1(M) && return false return ll_primecheck(M) end From 33dff5e24c3cf80691976801bbffed428a4ff7ed Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Sun, 28 Feb 2021 13:01:59 -0600 Subject: [PATCH 3/9] speed up ll_primecheck by using mod_mersenne --- src/Primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Primes.jl b/src/Primes.jl index 36fbdda..4f33579 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -557,7 +557,7 @@ function ll_primecheck(X::Integer, s::Integer = 4) S, N = BigInt(s), BigInt(ndigits(X, base=2)) X < 7 && throw(ArgumentError("The condition X ≥ 7 must be met.")) for i in 1:(N - 2) - S = (S^2 - 2) % X + S = mod_mersenne((S^2 - 2), N, X) end return S == 0 end From 983b45dec619e6f82f4fd1d27b887a245027889f Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Sun, 28 Feb 2021 13:04:26 -0600 Subject: [PATCH 4/9] fix another typo --- src/Primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Primes.jl b/src/Primes.jl index 4f33579..541d10c 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -460,7 +460,7 @@ function ismersenneprime(M::Integer; check::Bool = true) end M < 7 && return M == 3 has_small_factor_mersenne(M) && return false - p_minus1(M) && return false + p_minus1_mersenne(M) && return false return ll_primecheck(M) end From 79453130593626802e0fff05f01c31424736f907 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 1 Mar 2021 22:29:27 -0600 Subject: [PATCH 5/9] use linear sieve to dramatically improve runtime of has_small_factor_mersenne --- src/Primes.jl | 40 ++++++++++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index 541d10c..3785468 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -464,22 +464,42 @@ function ismersenneprime(M::Integer; check::Bool = true) return ll_primecheck(M) end +# Prime sieve for numbers of the form ax+b +# Filters out all composites with divisors less than 2^16 +function linear_sieve(a,b,max) + arr = falses(div((max-b),a)) + sqrtmax = isqrt(max) + for p in PRIMES + p > sqrtmax && break + g,x,_ = gcdx(a, p) + if g!=1 + if mod(b,p) == 0 + return typeof(max)[] + end + continue + end + # x*a = 1 mod p + k = mod(-b*x, p) + arr[k+1:p:end] .= true + # don't set p as composite. (only applies if arr contains elements as small as p) + if a*k+b == p + arr[k+1]=false + end + end + if b == 1 + arr[1] = true + end + return (a*(i-1)+b for (i,x) in enumerate(arr) if !x) +end + # Returns whether M=2^d-1 has a small factor # For each factor, checks if 2^d = 1 \pmod factor. # Uses the fact that all factors must be of the form # 2kd+1 for integer k function has_small_factor_mersenne(M::Integer) d = ndigits(M, base=2) - #limit = round(BigInt, 2^(22.85*d^(1/16))/(2*d)) - limit = round(BigInt, 7.3e-6 * d^2.32) - step = 2 * d - wheel = Iterators.cycle(d % 4 == 1 ? - (false,false,true,true) : - (true,false,false,true)) - for (factor, w) in zip(1+step : step : step*limit, wheel) - !w && continue + for factor in linear_sieve(2 * d, 1, 2^40) factor & 7 in (3,5) && continue - !isprime(factor, 4) && continue if powermod(2, d, factor) == 1 return true end @@ -579,7 +599,7 @@ function totient(f::Factorization{T}) where T <: Integer end result end - + 0 """ totient(n::Integer) -> Integer From 5c2efe59aa90a34c6b2d586f9181627bb940f298 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 1 Mar 2021 22:53:46 -0600 Subject: [PATCH 6/9] P-1 isn't worth the effort for practical sizes --- src/Primes.jl | 29 ----------------------------- 1 file changed, 29 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index 3785468..cba98a0 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -460,7 +460,6 @@ function ismersenneprime(M::Integer; check::Bool = true) end M < 7 && return M == 3 has_small_factor_mersenne(M) && return false - p_minus1_mersenne(M) && return false return ll_primecheck(M) end @@ -515,34 +514,6 @@ function mod_mersenne(n, d, M) return n != M ? n : 0 end -# Calculates base^exp % 2^d-1 where M=2**d-1 -function power_mod_mersenne(base, exp, d, M) - ans = 1 - while exp > 0 - if exp & 1 == 1 - ans = mod_mersenne(ans * base, d, M) - end - exp >>= 1 - base = mod_mersenne(base * base, d, M) - end - return ans -end - -# Returns if M=2^prime-1 have a factor 2*k*prime+1 -# such that the largest prime factor of k is less than limit""" -function p_minus1_mersenne(M) - d = ndigits(M, base=2) - B1 = d - log_B1 = log(B1) - P = prod(p^round(BigInt, log_B1/log(p)) for p in primes(B1)) - P = power_mod_mersenne(3, 2*P*d, d, M) - g = gcd(M, P-1) - 1 < g < M && return true - g == M && return false - return false - #return p_minus1_stage_2(prime, mersenne, P, p, 10*B1, primes) -end - """ isrieselprime(k::Integer, Q::Integer) -> Bool From e78537d97bf7a1b18965efdc2e66cc93ff89f8e6 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 2 Mar 2021 17:50:19 -0600 Subject: [PATCH 7/9] make additions a non-regression for small primes --- src/Primes.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index cba98a0..f2ed4e3 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -453,13 +453,14 @@ true ``` """ function ismersenneprime(M::Integer; check::Bool = true) + d = ndigits(M, base=2) if check - d = ndigits(M, base=2) M >= 0 && isprime(d) && (M >> d == 0) || throw(ArgumentError("The argument given is not a valid Mersenne Number (`M = 2^p - 1`).")) end M < 7 && return M == 3 - has_small_factor_mersenne(M) && return false + # If d is large, it is worth it to try to wead out easy cases + d > 10007 && has_small_factor_mersenne(M) && return false return ll_primecheck(M) end @@ -471,7 +472,7 @@ function linear_sieve(a,b,max) for p in PRIMES p > sqrtmax && break g,x,_ = gcdx(a, p) - if g!=1 + if g != 1 if mod(b,p) == 0 return typeof(max)[] end From 1cbbf9e778f1715ce81d6e7d064778b916686ea0 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Wed, 3 Mar 2021 19:38:43 -0600 Subject: [PATCH 8/9] add test for improvment. Will time out CI if broken, but take no time if working --- test/runtests.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 0cf78d6..6bb47a7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -242,6 +242,8 @@ end # Lucas-Lehmer @test !ismersenneprime(2047) @test ismersenneprime(8191) +# With has_small_factor_mersenne this should only take about .01 seconds +@test !ismersenneprime(big(2)^1257953-1) @test_throws ArgumentError ismersenneprime(9) @test_throws ArgumentError ismersenneprime(9, check=true) # test the following does not throw From 09c319a49eb6ff2abea5d5323e4da4604cd47433 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 4 Mar 2021 09:24:59 -0600 Subject: [PATCH 9/9] use prime that hassmallfactor will catch and fix 32 bit build --- src/Primes.jl | 2 +- test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index f2ed4e3..9698001 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -498,7 +498,7 @@ end # 2kd+1 for integer k function has_small_factor_mersenne(M::Integer) d = ndigits(M, base=2) - for factor in linear_sieve(2 * d, 1, 2^40) + for factor in linear_sieve(2 * d, 1, Int64(2)^40) factor & 7 in (3,5) && continue if powermod(2, d, factor) == 1 return true diff --git a/test/runtests.jl b/test/runtests.jl index 6bb47a7..817e8c2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -243,7 +243,7 @@ end @test !ismersenneprime(2047) @test ismersenneprime(8191) # With has_small_factor_mersenne this should only take about .01 seconds -@test !ismersenneprime(big(2)^1257953-1) +@test !ismersenneprime(big(2)^1257869-1) @test_throws ArgumentError ismersenneprime(9) @test_throws ArgumentError ismersenneprime(9, check=true) # test the following does not throw