Prime number algo... what's wrong?

[A.M. Kuchling]

On Sat, 22 Mar 2003 13:08:37 GMT, 
	L. B. <lorenzo at mysurname.net> wrote:
> Also are you aware of more sophisticated and efficient algos for prime
> testing implemented in Python?

Here's a straightforward implementation of the Rabin-Miller probabilistic
primality test.  The Python cryptography toolkit contains a faster
implementation of it that's less readable.  mathworld.wolfram.com has an
explanation of the Rabin-Miller test; textbooks on algorithms or number
theory may also discuss it.

(Has anyone implemented the polynomial-time primality test discovered a few
months back?  Wonder how complicated it is...)

--amk                                                    (www.amk.ca)
HAMLET: Angels and ministers of grace defend us!
      -- _Hamlet_, I, iv


def rm_test(n):
    # Say that 1 and 2 are prime.
    if n==1 or n==2:
        return 1

    # If even, it's obviously composite.
    if (n % 2) == 0:
        return 0

    # Straightforward Rabin-Miller test...
    # Compute r,s, such that N = s*2**r + 1.
    n1 = n-1
    r = 1
    while 1:
        div = n1/(2**r)
        if div*(2**r) != n1:
            break
        r += 1

    r = r-1 ; s = n1/(2**r)

    # Try a bunch of integers.  If a**s == 1, or
    # a**(s*2**j) == -1, for j in 0...r, the
    # number is prime.  (Trying more values of 'a' 
    # will decrease the risk that n isn't really prime.)
    for a in range(2, min(10,n)):
        t1 = pow(a,s,n)
        if t1 == 1:
            continue
        prime = 0
        for j in range(0, r):
            t2 = pow(a, s*(2**j), n)
            if t2 == (n-1):
                prime=1
                break
        else:
            return 0
    else:
        return 1

# Test small numbers
for i in range(1, 255, 2):
    p = rm_test(i)
    if p:
        print i,
print

# Look for a 256-bit prime
i = 2**256+1
while not rm_test(i):
    i += 2
print 'Prime:', i