primetests_mr.rb

require 'openssl'
require 'parallel'

class Integer
  # Returns true if +self+ passes Miller-Rabin Test on +b+
  def miller_rabin_test(b)             # b is a witness to test with
    return self >= 2 if self <= 3 
    return false unless 6.gcd(self % 6) == 1
    n = d = self - 1
    d >>= 1 while d.even?
    y = b.to_bn.mod_exp(d, self)       # y = (b**d) mod self
    while d != n && y != n && y != 1
      y = y.to_bn.mod_exp(2, self)     # y = (y**2) mod self
      d <<= 1
    end
    y == n || d.odd?
  end

  # Returns true if +self+ is a prime number, else returns false.
  def prime?
    return true if [2, 3, 5, 7].include? self
    return false unless self > 1 and 210.gcd(self % 210) == 1
    Prime::A014233.each do |pair|
      ix, mx = pair
      return false unless miller_rabin_test(ix)
      break if self < mx
    end
    self == 41 || miller_rabin_test(41)
  end
  
  def prime1?
    return true if [2, 3, 5, 7].include? self
    return false unless self > 1 and 210.gcd(self % 210) == 1
    witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    val = Prime::A014233.values.sort.detect {|v| v > self}
    witnesses.select! {|p| p <= Prime::A014233.key(val)} if val
    witnesses.each { |p| return false unless miller_rabin_test(p) }
    true
  end

  def prime2?    # perform miller-rabin tests in parallel
    return true if [2, 3, 5, 7].include? self
    return false unless self > 1 and 210.gcd(self % 210) == 1
    witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    val = Prime::A014233.values.sort.detect {|v| v > self}
    witnesses.select! {|p| p <= Prime::A014233.key(val)} if val
    Parallel.each(witnesses) { |p| return false unless miller_rabin_test(p) }
    true
  end

  def prime3?    # perform miller-rabin tests in parallel
    return true if [2, 3, 5, 7].include? self
    return false unless self > 1 and 210.gcd(self % 210) == 1
    witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    wits = Prime::A014233NEW.sort.detect {|range, wits| range > self} # [k, [wits]] or nil
    Parallel.each(wits && wits[1] || witnesses) { |p| return false unless miller_rabin_test(p) }
    true
  end

  def prime4?    # perform miller-rabin tests in parallel
    witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    return true if witnesses.include? self
    m = witnesses.reduce(:*)
    return false unless self > 1 and m.gcd(self % m) == 1
    wits = Prime::A014233NEW.sort.detect {|range, wits| range > self} # [k, [wit_prms]] or nil
    Parallel.each(wits && wits[1] || witnesses) { |p| return false unless miller_rabin_test(p) }
    true
  end

  def prime5?    # perform miller-rabin tests in parallel
    witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    return true if witnesses.include? self
    m = witnesses.reduce(:*)
    return false unless self > 1 and m.gcd(self % m) == 1
    wits = Prime::A014233NEW.sort.detect {|range, wits| range > self} # [k, [wit_prms]] or nil
    Parallel.each(wits && wits[1] || witnesses) do |p| 
      raise Parallel::Break unless miller_rabin_test(p)
    end
    true
  end
end

class Prime

  # https://oeis.org/A014233
  #
  # Smallest odd number for which Miller-Rabin primality test
  # on bases <= n-th prime does not reveal compositeness.
  #
  A014233 = {
    2  => 2_047,
    3  => 1_373_653,
    5  => 25_326_001,
    7  => 3_215_031_751,
    11 => 2_152_302_898_747,
    13 => 3_474_749_660_383,
    17 => 341_550_071_728_321,
    19 => 341_550_071_728_321,
    23 => 3_825_123_056_546_413_051,
    29 => 3_825_123_056_546_413_051,
    31 => 3_825_123_056_546_413_051,
    37 => 318_665_857_834_031_151_167_461
  }

  A014233NEW = {
    2_047 => [2],
    1_373_653 => [2, 3],
    25_326_001 => [2, 3, 5],
    3_215_031_751 => [2, 3, 5, 7],
    2_152_302_898_747 => [2, 3, 5, 7, 11],
    3_474_749_660_383 => [2, 3, 5, 7, 11, 13],
    341_550_071_728_321 => [2, 3, 5, 7, 11, 13, 17],
    3_825_123_056_546_413_051 => [2, 3, 5, 7, 11, 13, 17, 19, 23],
    318_665_857_834_031_151_167_461 => [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
  }
end

def tm; s=Time.now; yield; Time.now-s end