[Edu-sig] Miller-Rabin (second draft, much faster)

[Kirby Urner]

Chuckle.  Thanks to Tim Peters, my 2nd draft of Miller-Rabin
computes the same result in 0.3 seconds instead of 90 seconds.  
And if I take out those "I'm still doin' stuff" print statements, 
it drops to 0.06 or 0.07 seconds.  Heh heh heh... (rubbing hands).

So let's try that biggie prime from the paper I cited, 
57**(2**502)+1 (a decimal with 153 digits):

>>> bignum = 57*(1L<<502)+1
>>> bignum
74633305860032034636300725087669260670539438649781
87720021904319259185055802657985133855810503061478
30402163981083696190101534492461869616158904290377
729L
>>> rsa.ppt(bignum)
0.999755859375 chance of being prime
Elapsed time: 149.893572638 seconds

Under 3 minutes.  But Pete Emerson's thesis program ran 
through 20 iterations w/ successive prime bases vs. my 
default 6.  

I'll try that too:

>>> rsa.ppt(bignum,20)
0.999999999999 chance of being prime 
Elapsed time: 508.681804925 seconds

(I see I need to fix that print statement to show more 
significant digits)

OK, almost 9 minutes -- more than 3 times slower than 
the thesis program, but still, I'm in the ball park, within 
range of being able to check those 100-randomly-generated-
digit numbers for probable primehood.  And there may be 
some more easy optimizations I'm missing here too.

So here's the code for the 2nd draft (more comments after):

import primes, time

def ppt(n,i=6):
    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
    for b in bases:
        if n%b==0:
            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 >>= 1
        r += 1
    print "(n-1,r,m) = (%s,%s,%s)" % (2L**r*m,r,m)        
    for b in bases:
        tests += 1
        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./(4L**tests)))
    print "Elapsed time: %s seconds" % (time.clock()-start)

def miller(m,r,b,n):
    result = 0    
    testval = pow(b,m,n)
    # Miller test: if b**m mod n = -1 or 1, pass, or...
    if testval==1 or testval==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 = (1L<<k)*m
          if pow(b,expo,n)==n-1:
              result = 1
              break
    #if result:  print "Passed for base %s" % b
    #else: print "Failed for base %s" % b
    return result

>> = me
>  = Tim's comments

>>     # 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
>
> *Try* writing a nice, readable Python loop -- it might grow on you <wink>.
> Note that since you're using primes as bases, gcd(n,b[i]) == 1 is equivalent
> to the simpler n % b[i] != 0.  So:
>
>    for base in bases:
>        if n % base == 0:
>            print "Definitely not prime"
>            return
>
> from-the-convoluted-to-the-obvious-ly y'rs  - tim

Hah hah!  Of COURSE (duh).

looking-wonderingly-at-tortured-code-the-next-morning-ly y'rs -- Kirby