primes.py

import random


random = random.SystemRandom()


def naive_is_prime(n):
    if n <= 1:
        return False
    if n == 2:
        return True
    for i in range(2, int(n ** 0.5) + 1):
        if n % i == 0:
            return False
    return True


def is_prime(n, k=64):
    """
    Test whether n is prime probabilisticly.
    This uses the Miller-Rabin primality test. If n is composite,
    then this test will declare it to be probably prime with a
    probability of at most 4**-k.
    To be on the safe side, a value of k=64 for integers up to
    3072 bits is recommended (error probability = 2**-128). If
    the function is used for RSA or DSA, NIST recommends some
    values in FIPS PUB 186-3:
    <http://csrc.nist.gov/publications/fips/fips186-3/fips_186-3.pdf>
    """
    def check_candidate(a, d, n, s):
        if pow(a, d, n) == 1:
            return False
        for i in range(s):
            if pow(a, 2 ** i * d, n) == n - 1:
                return False
        return True
    if n < 100000000:
        return naive_is_prime(n)
    for i in range(3, 2048):  # performace optimisation
        if n % i == 0:
            return False
    s = 0
    d = n - 1
    while True:
        q, r = divmod(d, 2)
        if r == 1:
            break
        s += 1
        d = q
    for i in range(k):
        a = random.randint(2, n - 1)
        if check_candidate(a, d, n, s):
            return False
    return True