# millerRabin.py
# Miller-Rabin primality test
# https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

'''
>>> isPrime(13)
True

>>> isPrime(36557580306963713)
True

>>> isPrime(5166523324060302275610181)
False
'''

import random 


def isPrime(n, numTries=4):
  if (n <= 1) or (n == 4):
    return False;
  if n <= 3:
    return True;

  # Find r such that n =  2^d * r + 1 
  # for some r >= 1
  d = n - 1;
  while (d % 2) == 0:
    d //= 2;

  for i in range(numTries):
    if not mrTest(d, n):
      return False;
  return True;


def mrTest(d, n):
  '''
    returns false if n is composite and 
    returns true if n is probably prime. 
    d is an odd number such that d*2^r = n-1
    for some r >= 1
  '''

  a = random.randint(2, n-2);
  x = modPow(a, d, n);
  if (x == 1) or (x == (n-1)):
    return True;

  # Keep squaring x while 
  # (i) d does not reach n-1
  # (ii) (x^2) % n is not 1
  # (iii) (x^2) % n is not n-1
  while d != (n-1):
    x = (x * x) % n;
    d *= 2;
    if x == 1:
      return False;
    if x == (n - 1):
      return True;

  return False;


def modPow(x, n, mod):
  # returns (x^n) % mod
  res = 1; 
  x = x % mod; 
  while (n > 0):
    # If n is odd, multiply x with result
    if (n & 1):
      res = (res * x) % mod;
    # n must be even now
    n = n >> 1;   # n = n/2
    x = (x * x) % mod;
  return res;


if __name__ == "__main__":
  print("All primes smaller than 100: ");
  for n in range(1,100):
    if (isPrime(n, 4)):
      print(n , end=" ");