The Miller-Rabin primality test.

MillerRabin(n, s = 1000) -> bool Checks whether n is prime or not. This is an extremley fast algorithm designed to test very large numbers.

s is the number of tests to perform. The chance that Rabin-Miller is mistaken about a number (i.e. thinks it's prime, but it's not) is 2^(-s). So, a value of 50 for s is more than enough for any imaginable goal (2^(-50) is 8.8817841970012523e-16).

   1  
   2  import sys
   3  import random
   4  
   5  def toBinary(n):
   6    r = []
   7    while (n > 0):
   8      r.append(n % 2)
   9      n = n / 2
  10    return r
  11  
  12  def test(a, n):
  13    """
  14    test(a, n) -> bool Tests whether n is complex.
  15  
  16    Returns:
  17      - True, if n is complex.
  18      - False, if n is probably prime.
  19    """
  20    b = toBinary(n - 1)
  21    d = 1
  22    for i in xrange(len(b) - 1, -1, -1):
  23      x = d
  24      d = (d * d) % n
  25      if d == 1 and x != 1 and x != n - 1:
  26        return True # Complex
  27      if b[i] == 1:
  28        d = (d * a) % n
  29    if d != 1:
  30      return True # Complex
  31    return False # Prime
  32  
  33  def MillerRabin(n, s = 50):
  34    """
  35      MillerRabin(n, s = 1000) -> bool Checks whether n is prime or not
  36  
  37      Returns:
  38        - True, if n is probably prime.
  39        - False, if n is complex.
  40    """
  41    for j in xrange(1, s + 1):
  42      a = random.randint(1, n - 1)
  43      if (test(a, n)):
  44        return False # n is complex
  45    return True # n is prime
  46  
  47  def main(argv):
  48    print MillerRabin(int(argv[0]), int(argv[1]))
  49  
  50  if __name__ == "__main__":
  51    main(sys.argv[1:])