From 1866d483dce614a02c5186bd0588b48a5041e55e Mon Sep 17 00:00:00 2001 From: Marc-Andre Lafortune Date: Sat, 5 Dec 2020 00:20:39 -0500 Subject: [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 --- lib/prime.rb | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) (limited to 'lib/prime.rb') diff --git a/lib/prime.rb b/lib/prime.rb index d2de8dc017..43c6f9e943 100644 --- a/lib/prime.rb +++ b/lib/prime.rb @@ -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 # -- cgit v1.2.3