Rabin-Miller primality test implementation

February 21st, 2009 at 12:19 pm

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

        trials:
            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
            purposes.
    """
    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 $$2^{51}$$ – this is very low.

Therefore, you should always use this method instead of the naive one (trying do divide by all primes up to $$\sqrt{N}$$), because it’s much faster.

Related posts:

  1. hard IQ test
  2. Python unit testing: parametrized test cases
  3. SICP section 1.2.6
  4. Book review: “Test driven development by example”, Kent Beck
  5. did I finally “grok” Test:: ?

11 Responses to “Rabin-Miller primality test implementation”

  1. TheranNo Gravatar Says:

    You can also use the three-argument form of python’s builtin pow() function to compute a^(n-1) mod n. Just make sure you don’t hide it by importing the two argument pow from the math module with “from math import *”.

  2. elibenNo Gravatar Says:

    @Theran,

    pow is faster than a custom built modular exponentation, but for Rabin-Miller I need to perform the non-trivial-root test after each squaring, so I don’t think I can use it.

  3. MartinNo Gravatar Says:

    An alternative implementation of the bit extraction function:
    def _bits_of_n(n):
    return [int(b) for b in bin(n)[2:]]

  4. elibenNo Gravatar Says:

    @Martin,

    Yes, but this is Python 2.6 only. I’m on 2.5 for now :-)

  5. Guilherme BoaviagemNo Gravatar Says:

    I didn’t understand how the function “_MR_composite_witness” applies Miller Rabin test: where is the representation of n as (2^s)*d, d odd? Have you done each congruence of the test (“a^d = 1 mod n”, and “a^(d*2^i) = -1 mod n”, with i between 0 and s-1)?
    I need to test primality of numbers with more than 20 decimal digits, and for that I tried to implement Miller Rabin test in Python exactly as it is, but my algorithm is too much inefficient.

  6. elibenNo Gravatar Says:

    Guilherme,

    It implements modular exponentiation of a^(n-1) mod n by an iterative algorithm. See here: http://eli.thegreenplace.net/2009/03/28/efficient-modular-exponentiation-algorithms/

  7. amnNo Gravatar Says:

    is there a particular reason why it seems to fail on checking 97?
    came across this problem as i am trying to implement Miller-Rabin in SML and it seems to work for most of the primes i try (returns true) but not for 97.

    Any ideas as to why? increasing the number of trials (to 90) in your implementation apparently does not help :/

    i thought Miller-Rabin had false positives (on Pseudoprimes, not false negatives on primes?) or am I confusing something and 97 is not a real prime? o.O

  8. elibenNo Gravatar Says:

    amn,

    97 is prime and the Python implementation presented here correctly returns True for 97.

  9. amnNo Gravatar Says:

    not on my system o.O
    =(

  10. amnNo Gravatar Says:

    oh, now it does o.O

  11. John McNallyNo Gravatar Says:

    Thanks Eli,

    I had been looking for an efficiently implemented statistical primality test for Python.

    Your code works great in Python 2.7.1 — very brisk.