Miller Rabin primality test

miller_rabin.py

import random

def is_probably_prime(n, num_iterations):
    s = bin(n - 1)[::-1].index('1')
    d = n // 2**s

    for k in range(num_iterations):
        a = random.randint(2, n - 2)
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
            
        for r in xrange(1, s):
            x = pow(x, 2, n)
            if x == 1:
                return False
            if x == n - 1:
                break
        else:
            return False
    return True

There are some caveats/corner cases. E.g. n has to be > 3.