Miller-Rabin code - can't find any mistakes? [closed]

Question

Closed. This question needs details or clarity. It is not currently accepting answers.

---

Want to improve this question? Add details and clarify the problem by editing this post.

Closed 9 years ago.

Improve this question

I used some code taken from Rosetta Code. I renamed a few things but I didn't really change anything.

import random

def is_probable_prime(n, num_trials = 5):
    assert n >= 2
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    s = 0
    d = n-1
    while True:
        quotient, remainder = divmod(d, 2)
        if remainder == 1:
            break
        s += 1
        d = quotient
    assert(2**s * d == n-1)

    def try_composite(a):
        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

    for i in range(num_trials):
        a = random.randrange(2, n)
        if try_composite(a):
            return False
    return True

It pretty closely matches some psuedocode. However, when I test the number

123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901

it returns False. Other (python and java) implementations of Miller-Rabin return True for probable prime. After some testing, try_composite returns True after only 2 rounds! I would really like to know any error, I'm guessing a mis-indent or some feature I don't know about.

Answer 1

In your try_composite function, the for loop should be for i in range(1,s). Do not test the case where i is zero.

EDIT: Also, you are missing a test in your try_composite function. Here is my version of the pseudocode:

def isPrime(n, k=5):
    def isComposite(s, d):
        x = pow(randrange(2,n-1), d, n)
        if x == 1 or x == n-1: return False
        for r in range(1, s):
            x = pow(x, 2, n)
            if x == 1: return True
            if x == n-1: return False
        return True
    if n < 2: return False
    for p in [2, 3, 5, 7, 11, 13, 17]:
        if n % p == 0: return n == p
    s, d = 0, n-1
    while d % 2 == 0: s, d = s+1, d/2
    for i in range(k):
        if isComposite(s, d): return False
    return True

It's too bad that Python doesn't allow labels on break or continue statements. Here's pseudocode for a much prettier version of the function:

function isPrime(n, k=5)
    if n < 2 then return False
    for p in [2,3,5,7,11,13,17,19,23,29]
        if n % p == 0 then return n == p
    s, d = 0, n-1
    while d % 2 == 0
        s, d = s+1, d/2
    for i from 0 to k
        x = powerMod(randint(2, n-1), d, n)
        if x == 1 or x == n-1 then next i
        for r from 1 to s
            x = (x * x) % n
            if x == 1 then return False
            if x == n-1 then next i
        return False
    return True

Notice the two places where the control flow goes to next i. There is no good way to write that in Python. One choice uses an extra boolean variable that can be set and tested to determine when to bypass the rest of the code. The other choice, which I took above, is to write a local function to perform the task. This "loop-and-a-half" idiom is convenient and useful; it was proposed in PEP 3136 and rejected by Guido.