From 54f3ed17b542ceb03eb32a1da888a763f393e35b Mon Sep 17 00:00:00 2001
From: NagaChaitanya Vellanki <nagachaitanya.vellanki@gmail.com>
Date: Thu, 8 Oct 2015 20:42:53 -0700
Subject: [PATCH] Add Prime.probably_prime?

method to test if a given number is probably prime
or composite using Miller-Rabin Test
---
 lib/prime.rb       | 53 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 test/test_prime.rb | 25 +++++++++++++++++++++++++
 2 files changed, 78 insertions(+)

diff --git a/lib/prime.rb b/lib/prime.rb
index c9575f6..897b30e 100644
--- a/lib/prime.rb
+++ b/lib/prime.rb
@@ -145,6 +145,42 @@ class Prime
     end
   end
 
+  # Miller-Rabin Test to check if a given number is probably prime
+  #  or composite
+  # Returns true if n is probably prime, false if n is composite
+  # == Parameters
+  #
+  #  +n+:: an arbitary integer to be checked
+  #  +k+:: optional, A parameter that determines the accuracy
+  #                  of the test
+  #  References: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
+  #              http://rosettacode.org/wiki/Miller-Rabin_primality_test#Ruby
+  def probably_prime?(n, k = 30)
+    raise ArgumentError, 'n should be >= 3' if n < 3
+    d = n - 1
+    r = 0
+    while d & 1 == 0 do
+      d = d / 2
+      r += 1
+    end
+
+    k.times do
+      a = 2 + rand(n - 4)
+      x = mod_exp(a, d, n)
+      next if x == 1 || x == (n - 1)
+
+      (r - 1).times do
+        x = mod_exp(x, 2, n)
+        return false if x == 1
+        break if x == (n - 1)
+      end
+      return false if x != (n - 1)
+    end
+
+    true
+  end
+
+
   # Re-composes a prime factorization and returns the product.
   #
   # == Parameters
@@ -453,4 +489,21 @@ class Prime
       @max_checked = segment_max
     end
   end
+
+  private
+
+  # returns (base ^ exp) % n
+  # modular exponentiation by squaring
+  # https://en.wikipedia.org/wiki/Modular_exponentiation
+  def mod_exp(base, exp, n)
+    product = 1
+    base = base % n
+    while exp != 0 do
+      product = (product * base) % n if exp & 1 == 1
+      exp = exp >> 1
+      base = (base * base) % n
+    end
+
+    product
+  end
 end
diff --git a/test/test_prime.rb b/test/test_prime.rb
index 1788d71..9cc3b41 100644
--- a/test/test_prime.rb
+++ b/test/test_prime.rb
@@ -190,4 +190,29 @@ class TestPrime < Test::Unit::TestCase
 
     assert_not_include Prime.each(7*37).to_a, 7*37, "[ruby-dev:39465]"
   end
+
+  def test_probably_prime
+    assert_equal Prime.probably_prime?(961_748_941), true
+    assert_equal Prime.probably_prime?(4547337172376300111955330758342147474062293202868155909489), true
+    assert_raise ArgumentError do
+      Prime.probably_prime?(2)
+    end
+    PRIMES[1, PRIMES.length - 1].each do |prime|
+      assert_equal Prime.probably_prime?(prime), true
+    end
+
+    assert_equal Prime.probably_prime?(314159 * (10 ** 765) + 951413), true
+    assert_equal Prime.probably_prime?(1749343240116807117649823543576480140353607282475249), true
+    assert_equal Prime.probably_prime?(999999999999999999999999999999999841), true
+    assert_equal Prime.probably_prime?(25262728293031323334353637383940414243444546474849), true
+    assert_equal Prime.probably_prime?(33452526613163807108170062053440751665152000000001), true
+    assert_equal Prime.probably_prime?(1808422353177349564546512035512530001279481259854248860454348989451026887), true
+    assert_equal Prime.probably_prime?(13579111_3151719313_3353739515_3555759717_3757779919_3959799111_1131151171_1913113313_5137139151_1531551571_5917117317_5177179191_1931951971_9931131331_5317319331_3333353373_3935135335_5357359371_3733753773_7939139339_5397399511_5135155175_1953153353_5537539551_5535555575_5957157357_5577579591_5935955975_9971171371_5717719731_7337357377_3975175375_5757759771, 1000), true
+    assert_equal Prime.probably_prime?(98765432101456789), true
+    assert_equal Prime.probably_prime?(99194853094755497), true
+    assert_equal Prime.probably_prime?((10 ** 100) - 1), false
+    assert_equal Prime.probably_prime?(2 ** 256), false
+    assert_equal Prime.probably_prime?(9 ** 512), false
+    assert_equal Prime.probably_prime?(13579111_3151719313_3353739515_3555759717_3757779919_3959799111_1131151171_1913113313_5137139151_1531551571_5917117317_5177179191_1931951971_9931131331_5317319331_3333353373_3935135335_5357359371_3733753773_7939139339_5397399511_5135155175_1953153353_5537539551_5535555575_5957157357_5577579591_5935955975_9971171371_5717719731_7337357377_3975175375_5757759770, 1000), false
+  end
 end
-- 
2.3.8 (Apple Git-58)