Here's a fairly efficient Python (2.5) and well-documented implementation of the Rabin-Miller primality test, based on section 33.8 in CLR's Introduction to Algorithms. Due to Python's built-in arbitrary precision arithmetic, this works for numbers of any size.

from random import randint

def _bits_of_n(n):
    """ Return the list of the bits in the binary
        representation of n, from LSB to MSB
    bits = []

    while n:
        bits.append(n % 2)
        n /= 2

    return bits

def _MR_composite_witness(a, n):
    """ Witness functions for the Miller-Rabin
        test. If 'a' can be used to prove that
        'n' is composite, return True. If False
        is returned, there's high (though < 1)
        probability that 'n' is prime.
    rem = 1

    # Computes a^(n-1) mod n, using modular
    # exponentation by repeative squaring.
    for b in reversed(_bits_of_n(n - 1)):
        x = rem
        rem = (rem * rem) % n

        if rem == 1 and x != 1 and x != n - 1:
            return True

        if b == 1:
            rem = (rem * a) % n

    if rem != 1:
        return True
    return False

def isprime_MR(n, trials=6):
    """ Determine whether n is prime using the
        probabilistic Miller-Rabin test. Follows
        the procedure described in section 33.8
        in CLR's Introduction to Algorithms

            The amount of trials of the test.
            A larger amount of trials increases
            the chances of a correct answer.
            6 is safe enough for all practical
    if n < 2:
        return False

    for ntrial in xrange(trials):
        if _MR_composite_witness(randint(1, n - 1), n):
            return False

    return True

The function you should call is isprime_MR.

Although this test is probabilistic, the chances of it erring are extremely low. According to Bruce Schneier in "Applied Cryptography", the chances of error for a 256-bit number with 6 trials are less than one in - this is very low.

Therefore, you should always use this method instead of the naive one (trying do divide by all primes up to ), because it's much faster.


comments powered by Disqus