diff --git a/lib/prime.rb b/lib/prime.rb index 5be12f24f5..2d1267d2c7 100644 --- a/lib/prime.rb +++ b/lib/prime.rb @@ -32,15 +32,48 @@ def prime_division(generator = Prime::Generator23.new) # Returns true if +self+ is a prime number, else returns false. def prime? - return self >= 2 if self <= 3 - return true if self == 5 - return false unless 30.gcd(self) == 1 - (7..Integer.sqrt(self)).step(30) do |p| - return false if - self%(p) == 0 || self%(p+4) == 0 || self%(p+6) == 0 || self%(p+10) == 0 || - self%(p+12) == 0 || self%(p+16) == 0 || self%(p+22) == 0 || self%(p+24) == 0 - end - true + # https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test + # https://arxiv.org/pdf/1509.00864.pdf + # + # if n < 3,317,044,064,679,887,385,961,981, + # it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, and 41. + # and the result will be deterministic. + # + + if self > 3317044064679887385961981 + raise ArgumentError, ".prime? is only valid for values less than 3,317,044,064,679,887,385,961,981" + end + + return false if self < 2 + return true if self < 4 + return false if self.even? + + r = 0 + d = self - 1 + + while d.even? + d /= 2 + r += 1 + end + + [2,3,5,7,11,13,17,19,23,29,31,37,41].each do |a| + x = a.pow(d, self) + next if x == 1 || x == (self-1) || a == self + + (r-1).times do + x = x.pow(2, self) + if x == (self-1) + break + end + end + + if (x == self-1) + next + else + return false + end + end + return true end # Iterates the given block over all prime numbers. diff --git a/test/test_prime.rb b/test/test_prime.rb index 9db13f08fe..e057fd2049 100644 --- a/test/test_prime.rb +++ b/test/test_prime.rb @@ -219,7 +219,11 @@ def test_from_prime_division assert_equal(-PRIMES.inject(&:*), Integer.from_prime_division([[-1, 1]] + PRIMES.map{|p| [p,1]})) end - + + def test_prime_bounds_check + assert_raise(ArgumentError) { 33170440646798873859619812.prime? } + end + def test_prime? PRIMES.each do |p| assert_predicate(p, :prime?) @@ -241,9 +245,6 @@ def test_prime? # large composite assert_not_predicate(((2**13-1) * (2**17-1)), :prime?) - # factorial - assert_not_predicate((2...100).inject(&:*), :prime?) - # negative assert_not_predicate(-1, :prime?) assert_not_predicate(-2, :prime?)