[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>

1 files changed, 74 insertions, 0 deletions

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
   #