[Edu-sig] Miller-Rabin (first draft, slow)

Thu, 14 Dec 2000 10:51:26 -0800

>not) -- I don't have the tests for that written.  A Fermat test 
>seems impractical.  Has anyone done Miller-Rabin in pure Python, 
>for pedagogic purposes?  Still looking into it.  Yes, I know, 
>I know, I need to join the Python cryptography SIG.
>
>Kirby

I haven't chatted with the Python crypto heads yet, but using
http://www.middlebury.edu/~pemerson/thesis/node5.html and the
related source code, I was able to come up with a first draft
of Miller-Rabin's primality test in pure Python.

Although I believe it implements the algorithm "by the book",
it's quite slow -- no way I'd want to sic this on some randomly
generated 100-digit possible-prime.  Like, it took a minute and 
a half just to verify that 399173 was very likey prime -- 
something lowly trial-by-division can handle (with certainty)
in the blink of an eye.

The paper's author, Paul Emerson, was getting his program to 
verify 57*2**502 +1 in under 3 minutes, checking against 20 bases.  
There's gotta be some drastic difference in our hardware and/or 
software in that case.  His is written in C and uses mpz utilities.  
I hypothesize that I should likewise be using Python's optional 
mpz module -- but that one doesn't come with the default binary 
for Windows.  Maybe in my Linux RPM?

Note, in the code below, some of the functions are actually defined 
in primes.py, which is at my website:  basically gcd(a,b) and 
get2nb(n) -- find the greatest common divisor (Euclid's Algorithm) 
and get n successive primes from 2 (trial-by-division).

See: http://www.python.org/pipermail/edu-sig/2000-August/000599.html
for more on Euclid's Algorithm (includes Python and DrScheme source), 
plus a discussion of precision issues.  I think GCD, LCM and 
Euclid's Algorithm in particular belong at the grade 6-8 level.

Miller-Rabin, along with other primality tests, should be mentioned
in high school, but possibly only in a story-telling context (along
with some play with the Euler-Fermat test -- see related reading
below), i.e. the point might be to give the flavor of number theory 
applied to cryptography -- a topic with many historical ramifications.

import primes, time 

def ppt(n,i=6):
    """
    Probable Primality Test (Miller-Rabin)
    """
    start = time.clock()
    # get i successive primes from 2 to < n
    bases  = filter(lambda x,n=n: x<n, primes.get2nb(i))
    isprob = 0
    # if any of the primes is a factor, we're done
    if 0 in map(lambda i,n=n,b=bases:
                primes.gcd(n,b[i])==1, range(len(bases))):
       print "Definitely not prime"
       return
    tests  = 0
    m = n-1
    r=0
    # turning (n-1) into form (2**r) * m
    while m%2==0:
        m = m/2
        r += 1
    print "(n-1,r,m) = (%s,%s,%s)" % (2L**r*m,r,m)        
    for b in bases:
        tests += 1
        print "Testing base %s" % b
        isprob = miller(m,r,b,n)
        if not isprob: break

    if not isprob:  print "Definitely not prime"
    else: print "%s chance of being prime" % (1-(1./(4**tests)))
    print "Elapsed time: %s seconds" % (time.clock()-start)

def miller(m,r,b,n):
    result = 0    
    bignum = b**m
    # Miller test: if b**m mod n = -1 or 1, pass, or...
    if bignum%n ==1 or bignum%n == n-1:
       result = 1
    else: # if b**(m*(2**k)) mod n = -1 for some 1<k<r, pass
       for k in range(1,r+1):
          expo = 2L**k
          if (bignum**expo)%n==n-1:
              result = 1
              break
    if result:  print "Passed for base %s" % b
    else: print "Failed for base %s" % b
    return result

Usage:

 >>> rsa.ppt(399173L)
 (n-1,r,m) = (399172,2,99793)
 Testing base 2
 Passed for base 2
 Testing base 3
 Passed for base 3
 Testing base 5
 Passed for base 5
 Testing base 7
 Passed for base 7
 Testing base 11
 Passed for base 11
 Testing base 13
 Passed for base 13
 0.999755859375 chance of being prime
 Elapsed time: 89.9671910357 seconds

Kirby