[ruby/prime] Optimize `Integer#prime?`

Miller Rabin algorithm can be used to test primality for integers smaller than a max value "MaxMR" (~3e24)

It can be much faster than previous implementation: ~100x faster for numbers with 13 digits, at least 5 orders of magnitude for even larger numbers (previous implementation is so slow that it requires more patience than I have for more precise estimate).

Miller Rabin test becomes faster than previous implementation at somewhere in the range 1e5-1e6. It seems that the range 62000..66000 is where Miller Rabin starts being always faster, so I picked 0xffff arbitrarily; before that, or above "MaxMR", the previous implementation remains.

I compared the `faster_prime` gem too. It is slower than previous implementation up to ~1e4. After that it becomes faster and faster compared to previous implementation, but is still slower than Miller Rabin starting at ~1e5 and up to MaxMR. Thus, after this commit, builtin `Integer#prime?` will be similar or faster than `faster_prime` up to "MaxMR".

Adapted from patch of Stephen Blackstone [Feature #16468]

Benchmark results and code: https://gist.github.com/marcandre/b263bdae488e76dabdda84daf73733b9

Co-authored-by: Stephen Blackstone <sblackstone@gmail.com>
This commit is contained in:
Marc-Andre Lafortune 2020-12-05 00:20:39 -05:00 коммит произвёл Marc-André Lafortune
Родитель dea600046a
Коммит 1866d483dc
2 изменённых файлов: 85 добавлений и 0 удалений

Просмотреть файл

@ -31,8 +31,14 @@ class Integer
end
# Returns true if +self+ is a prime number, else returns false.
# Not recommended for very big integers (> 10**23).
def prime?
return self >= 2 if self <= 3
if (bases = miller_rabin_bases)
return miller_rabin_test(bases)
end
return true if self == 5
return false unless 30.gcd(self) == 1
(7..Integer.sqrt(self)).step(30) do |p|
@ -43,6 +49,73 @@ class Integer
true
end
MILLER_RABIN_BASES = [
[2],
[2,3],
[31,73],
[2,3,5],
[2,3,5,7],
[2,7,61],
[2,13,23,1662803],
[2,3,5,7,11],
[2,3,5,7,11,13],
[2,3,5,7,11,13,17],
[2,3,5,7,11,13,17,19,23],
[2,3,5,7,11,13,17,19,23,29,31,37],
[2,3,5,7,11,13,17,19,23,29,31,37,41],
].map!(&:freeze).freeze
private_constant :MILLER_RABIN_BASES
private def miller_rabin_bases
# Miller-Rabin's complexity is O(k log^3n).
# So we can reduce the complexity by reducing the number of bases tested.
# Using values from https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
i = case
when self < 0xffff then
# For small integers, Miller Rabin can be slower
# There is no mathematical significance to 0xffff
return nil
# when self < 2_047 then 0
when self < 1_373_653 then 1
when self < 9_080_191 then 2
when self < 25_326_001 then 3
when self < 3_215_031_751 then 4
when self < 4_759_123_141 then 5
when self < 1_122_004_669_633 then 6
when self < 2_152_302_898_747 then 7
when self < 3_474_749_660_383 then 8
when self < 341_550_071_728_321 then 9
when self < 3_825_123_056_546_413_051 then 10
when self < 318_665_857_834_031_151_167_461 then 11
when self < 3_317_044_064_679_887_385_961_981 then 12
else return nil
end
MILLER_RABIN_BASES[i]
end
private def miller_rabin_test(bases)
return false if even?
r = 0
d = self >> 1
while d.even?
d >>= 1
r += 1
end
self_minus_1 = self-1
bases.each do |a|
x = a.pow(d, self)
next if x == 1 || x == self_minus_1 || a == self
return false if r.times do
x = x.pow(2, self)
break if x == self_minus_1
end
end
true
end
# Iterates the given block over all prime numbers.
#
# See +Prime+#each for more details.
@ -156,6 +229,7 @@ class Prime
end
# Returns true if +value+ is a prime number, else returns false.
# Integer#prime? is much more performant.
#
# == Parameters
#

Просмотреть файл

@ -257,7 +257,18 @@ class TestPrime < Test::Unit::TestCase
assert_not_predicate(-2, :prime?)
assert_not_predicate(-3, :prime?)
assert_not_predicate(-4, :prime?)
assert_equal 1229, (1..10_000).count(&:prime?)
assert_equal 861, (100_000..110_000).count(&:prime?)
end
=begin
# now Ractor should not use in test-all process
def test_prime_in_ractor
# Test usage of private constant...
assert_equal false, Ractor.new { ((2**13-1) * (2**17-1)).prime? }.take
end if defined?(Ractor)
=end
end
def test_eratosthenes_works_fine_after_timeout