miller-rabin.rb

require 'inline'

class MyMath

  # Thanks, Bruce Schneier!
  inline do |builder|
    builder.c <<-EOS
      unsigned long expmod(unsigned long base, unsigned long exponent, unsigned long modulus) {
        unsigned long result = 1;
        while(exponent > 0) {
          if ((exponent & 1) == 1) {
            result = (result * base) % modulus;
          }
          exponent >>= 1;
          base = (base * base) % modulus;
        }
      return result;
    }
    EOS
  end

end

class Fixnum
  
  def prime? # miller-rabin test
    n = self

    math = MyMath.new
    
    return false if n < 2
    return true  if n < 4

    if n < 1_373_653
      primes = [2, 3]
    elsif n < 9_080_191
      primes = [31, 73]
    elsif n < 4_759_123_141
      primes = [2, 7, 61]
    elsif n < 2_152_302_989_747
      primes = [2, 3, 5, 7, 11]
    elsif n < 3_474_749_660_383
      primes = [2, 3, 5, 7, 11, 13]
    elsif n < 341_550_071_728_321
      primes = [2, 3, 5, 7, 11, 13, 17]
    else
      raise ArgumentError, "Use a different primality test."
    end
    
    k = 0
    m = n - 1
    while m.even?
      m /= 2
      k += 1
    end

    primes.each do |a|
      b = math.expmod(a, m, n) # (a**m)%n
      
      next if b==1

      prime = false
      k.times do
        if b == n - 1
          prime = true
          break
        end
        b = math.expmod(b, 2, n) # (b*b)%n
      end
        
      return false unless prime
    end
   
    return true
  end
end