Feature #16468
Updated by steveb3210 (Stephen Blackstone) over 4 years ago
The miller-rabin algorithm is a non-deterministic primality test, however it is known that below 2**64, you can always get a deterministic answer by only checking a=[2,3,5,7,11,13,17,19,23, 29, 31, 37]
Given that Prime.prime? would never respond in a reasonable amount of time for larger numbers, we can gain much more utility and performance by switching..
```
user system total real
miller_rabin: random set 0.150000 0.000000 0.150000 ( 0.152212)
Prime.prime?: random set 0.270000 0.000000 0.270000 ( 0.281257)
user system total real
miller_rabin: 16 digits 0.010000 0.000000 0.010000 ( 0.000300)
Prime.prime? 16 digits 2.200000 0.020000 2.220000 ( 2.368247)
user system total real
miller_rabin: 2-10000 0.030000 0.000000 0.030000 ( 0.035752)
Prime.prime? 2-10000 0.020000 0.000000 0.020000 ( 0.022948)
```
```
require 'benchmark'
require 'prime'
def modpow(base, power, mod)
result = 1
while power > 0
result = (result * base) % mod if power & 1 == 1
base = (base * base) % mod
power >>= 1;
end
result
end
def miller_rabin(n)
return false if n < 2 # 0, 1
return true if n < 4 # 2, 3
return false if n % 2 == 0 # 4, 6, 8, 10
d = (n-1)
s = 0
while (d % 2 == 0)
s +=1
d /= 2
end
# https://arxiv.org/pdf/1509.00864.pdf
#
# if n < 18,446,744,073,709,551,616 = 2**64, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37.
#
[2,3,5,7,11,13,17,19,23, 29, 31, 37].each do |a|
x = modpow(a,d,n)
next if x == 1 or x == (n - 1) or n == a
skip = false
1.upto(s-1) do |k|
x = x*x % n
if x == 1
return false
elsif x == (n - 1)
skip = true
break
end
end
next if skip
return false
end
return true
end
test_set = []
while test_set.length < 50_000
test_set << rand(1_000_000)
end
c = 1
Benchmark.bm(40) do |bm|
bm.report("miller_rabin: random set") do
test_set.each do |x|
miller_rabin(x)
end
end
bm.report("Prime.prime?: random set") do
test_set.each do |x|
Prime.prime?(x)
end
end
end
puts
Benchmark.bm(40) do |bm|
bm.report("miller_rabin: 16 digits") do
miller_rabin(1000000000100011)
end
bm.report("Prime.prime? 16 digits") do
Prime.prime?(1000000000100011)
end
end
puts
Benchmark.bm(40) do |bm|
bm.report("miller_rabin: 2-10000") do
(2..10000).each do |x|
miller_rabin(x)
end
end
bm.report("Prime.prime? 2-10000") do
(2..10000).each do |x|
Prime.prime?(x)
end
end
end
```