[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
> = 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