From aac32f1966c68bce2208f632d8104da6309a6ea1 Mon Sep 17 00:00:00 2001 From: pabloferz Date: Fri, 27 May 2016 12:15:56 +0200 Subject: [PATCH 1/4] Improve isprime performance for 32 bit integers --- src/Primes.jl | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index 8625b08..648de7b 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -133,9 +133,13 @@ true function isprime(n::Integer) # Small precomputed primes + Miller-Rabin for primality testing: # https://en.wikipedia.org/wiki/Miller–Rabin_primality_test - # - (n < 3 || iseven(n)) && return n == 2 - n <= 2^16 && return PRIMES[searchsortedlast(PRIMES,n)] == n + # https://github.com/JuliaLang/julia/issues/11594 + iseven(n) && return n == 2 + n in (3,5,7,11,13,17,19,23) && return true + for m in (3,5,7,11,13,17,19,23) + n % m == 0 && return false + end + n < 841 && return n > 1 s = trailing_zeros(n-1) d = (n-1) >>> s for a in witnesses(n) @@ -171,12 +175,14 @@ isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr # http://mathoverflow.net/questions/101922/smallest-collection-of-bases-for-prime-testing-of-64-bit-numbers # http://primes.utm.edu/prove/merged.html # http://miller-rabin.appspot.com -# -witnesses(n::Union{UInt8,Int8,UInt16,Int16}) = (2,3) -witnesses(n::Union{UInt32,Int32}) = n < 1373653 ? (2,3) : (2,7,61) +# https://github.com/JuliaLang/julia/issues/11594 +function witnesses(n::Union{UInt8,Int8,UInt16,Int16,UInt32,Int32}) + i = (((n >> 16) $ n) * 0x45d9f3b) & 0xffffffff + i = ((i >> 16) $ i) & 15 + 1 + (2,(2249,483,194,199,15,369,499,945,419,735,33,471,946,615,497,702)[i]) +end witnesses(n::Union{UInt64,Int64}) = - n < 1373653 ? (2,3) : - n < 4759123141 ? (2,7,61) : + n < 4294967296 ? witnesses(UInt32(n)) : n < 2152302898747 ? (2,3,5,7,11) : n < 3474749660383 ? (2,3,5,7,11,13) : (2,325,9375,28178,450775,9780504,1795265022) From 0f9c62f10540d5f91330ca9ec52277aca1c19a78 Mon Sep 17 00:00:00 2001 From: pabloferz Date: Fri, 27 May 2016 15:48:13 +0200 Subject: [PATCH 2/4] Forisek and Jancina algorithm --- src/Primes.jl | 35 +++++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/src/Primes.jl b/src/Primes.jl index 648de7b..965dc32 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -142,7 +142,7 @@ function isprime(n::Integer) n < 841 && return n > 1 s = trailing_zeros(n-1) d = (n-1) >>> s - for a in witnesses(n) + for a in witnesses(n)::Tuple{Vararg{Int}} x = powermod(a,d,n) x == 1 && continue t = s @@ -170,19 +170,38 @@ true isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr{BigInt}, Cint), &x, reps) > 0 - # Miller-Rabin witness choices based on: # http://mathoverflow.net/questions/101922/smallest-collection-of-bases-for-prime-testing-of-64-bit-numbers # http://primes.utm.edu/prove/merged.html # http://miller-rabin.appspot.com # https://github.com/JuliaLang/julia/issues/11594 -function witnesses(n::Union{UInt8,Int8,UInt16,Int16,UInt32,Int32}) - i = (((n >> 16) $ n) * 0x45d9f3b) & 0xffffffff - i = ((i >> 16) $ i) & 15 + 1 - (2,(2249,483,194,199,15,369,499,945,419,735,33,471,946,615,497,702)[i]) +# Forišek and Jančina, "Fast Primality Testing for Integers That Fit into a Machine Word", 2015 +const bases = UInt16[15591,2018,166,7429,8064,16045,10503,4399,1949,1295,2776, +3620,560,3128,5212,2657,2300,2021,4652,1471,9336,4018,2398,20462,10277,8028, +2213,6219,620,3763,4852,5012,3185,1333,6227,5298,1074,2391,5113,7061,803,1269, +3875,422,751,580,4729,10239,746,2951,556,2206,3778,481,1522,3476,481,2487,3266, +5633,488,3373,6441,3344,17,15105,1490,4154,2036,1882,1813,467,3307,14042,6371, +658,1005,903,737,1887,7447,1888,2848,1784,7559,3400,951,13969,4304,177,41, +19875,3110,13221,8726,571,7043,6943,1199,352,6435,165,1169,3315,978,233,3003, +2562,2994,10587,10030,2377,1902,5354,4447,1555,263,27027,2283,305,669,1912,601, +6186,429,1930,14873,1784,1661,524,3577,236,2360,6146,2850,55637,1753,4178,8466, +222,2579,2743,2031,2226,2276,374,2132,813,23788,1610,4422,5159,1725,3597,3366, +14336,579,165,1375,10018,12616,9816,1371,536,1867,10864,857,2206,5788,434,8085, +17618,727,3639,1595,4944,2129,2029,8195,8344,6232,9183,8126,1870,3296,7455, +8947,25017,541,19115,368,566,5674,411,522,1027,8215,2050,6544,10049,614,774, +2333,3007,35201,4706,1152,1785,1028,1540,3743,493,4474,2521,26845,8354,864, +18915,5465,2447,42,4511,1660,166,1249,6259,2553,304,272,7286,73,6554,899,2816, +5197,13330,7054,2818,3199,811,922,350,7514,4452,3449,2663,4708,418,1621,1171, +3471,88,11345,412,1559,194] + +function _witnesses(n::UInt64) + i = ((n >> 16) $ n) * 0x45d9f3b + i = ((i >> 16) $ i) * 0x45d9f3b + i = ((i >> 16) $ i) & 255 + 1 + @inbounds (Int(bases[i]),) end -witnesses(n::Union{UInt64,Int64}) = - n < 4294967296 ? witnesses(UInt32(n)) : +witnesses(n::Integer) = + n < 4294967296 ? _witnesses(UInt64(n)) : n < 2152302898747 ? (2,3,5,7,11) : n < 3474749660383 ? (2,3,5,7,11,13) : (2,325,9375,28178,450775,9780504,1795265022) From 77debe1175b7895ffb79464565e68cb9127320f9 Mon Sep 17 00:00:00 2001 From: pabloferz Date: Fri, 27 May 2016 15:53:47 +0200 Subject: [PATCH 3/4] Add missing return --- src/Primes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Primes.jl b/src/Primes.jl index 965dc32..2228af7 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -198,7 +198,7 @@ function _witnesses(n::UInt64) i = ((n >> 16) $ n) * 0x45d9f3b i = ((i >> 16) $ i) * 0x45d9f3b i = ((i >> 16) $ i) & 255 + 1 - @inbounds (Int(bases[i]),) + @inbounds return (Int(bases[i]),) end witnesses(n::Integer) = n < 4294967296 ? _witnesses(UInt64(n)) : From 4dc4fae980967c32902e3b2b51ba902f08e65f3b Mon Sep 17 00:00:00 2001 From: pabloferz Date: Wed, 1 Jun 2016 14:12:40 +0200 Subject: [PATCH 4/4] Indicate algorithm used --- src/Primes.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Primes.jl b/src/Primes.jl index 2228af7..1a21069 100644 --- a/src/Primes.jl +++ b/src/Primes.jl @@ -176,6 +176,7 @@ isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr # http://miller-rabin.appspot.com # https://github.com/JuliaLang/julia/issues/11594 # Forišek and Jančina, "Fast Primality Testing for Integers That Fit into a Machine Word", 2015 +# (in particular, see function FJ32_256, from which the hash and bases were taken) const bases = UInt16[15591,2018,166,7429,8064,16045,10503,4399,1949,1295,2776, 3620,560,3128,5212,2657,2300,2021,4652,1471,9336,4018,2398,20462,10277,8028, 2213,6219,620,3763,4852,5012,3185,1333,6227,5298,1074,2391,5113,7061,803,1269,