Miller-Rabin code - can't find any mistakes?
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: 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](http://legacy.python.org/dev/peps/pep-3136/) and
[rejected](https://mail.python.org/pipermail/python-3000/2007-July/008663.html)
by Guido.

