Miller-Rabin primality in Python

gistfile1.py

import random

_rng = random.SystemRandom()

def modexp(base, exp, mod):
    '''
    Classic modular exponentiation by repeated squaring.
    
    '''
    acc = 1
    while exp != 0:
        if (exp%2) == 1:
            acc *= base
            acc %= mod
        base *= base
        base %= mod
        exp >>= 1
    return acc

def _factor_binary_powers(n):
    '''
    Computes n = 2^s * d as a helper for Miller-Rabin.

    '''
    s = 0
    d = n
    while d % 2 == 0:
        s += 1
        d /= 2
    return s, d


def is_prime(n, iterations = 16):
    '''
    Tests n for primality using the classic Miller-Rabin algorithm.

    '''
    s, d = _factor_binary_powers(n - 1)
    for _ in xrange(iterations):
        a = _rng.randint(2, n - 2)
        x = modexp(a, d, n)
        if x == 1 or x == (n-1):
            continue
        test_pass = False
        for _r in xrange(1, s+1):
            x = x*x % n
            if x == 1:
                return False
            if x == (n-1):
                test_pass = True
                break
        if not test_pass:
            return False
    return True