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

Kirby Urner pdx4d@teleport.com
Thu, 14 Dec 2000 21:03:42 -0800

```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

>>     # 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

```