
FYI (fer ya'll's information), I did a simple implementation of Pollard's rho method for factoring into primes at: http://www.mathforum.org/epigone/math-teach/blerlplerdghi (the 2nd post in that thread contains important bug fixes and enhancements). This includes a reimplementation of Miller-Rabin without eye balling my source of Dec 2000 which Tim Peters helped me speed up tremendously back then. Except just now I reread those old posts and decided to bring back another of Tim's suggestions: substitute 1L<<n for 2**n (included in the version below). Jumping to the end of that thread: http://mail.python.org/pipermail/edu-sig/2000-December/000827.html For continuity, I'm appending the newer source code here (the Miller-Rabin part). I'm pretty sure I'm using the same motherboard as in 2000 (1.2 Mhz) but now have more RAM (1 gig). Also, this is Python 2.3 and that was what? 2.0? These factors would likely contribute to speed gain. But I bet it's partly the way it's written that makes it somewhat faster. Note I'm not using any prime bases in a pre-test as per the thesis I was following on a web page long gone. Most implementations of Miller-Rabin just pick random integers for each trial (every time you pass a test, similar to a Fermat test, you're more likely to be prime) -- which is all I did here, using the library's random.randint(). from random import randint import time def millrab(n, max=30, timer=False): """ Miller-Rabin primality test as per the following source: http://www.wikipedia.org/wiki/Miller-Rabin_primality_test """ # returns probability p is prime: either p = 0 or ~1 if timer: start = time.clock() if not n%2: return 0 k = 0 z = n - 1 # compute m,k such that (2**k)*m = n-1 while not z % 2: k += 1 z //= 2 m = z ok = 1 trials = 0 p = 1 # try tests with max random integers between 2,n-1 while trials < max and ok: a = randint(2,n-1) trials += 1 test = pow(a,m,n) if (not test == 1) and not (test == n-1): # if 1st test fails, fall through ok = 0 for r in range(1,k): test = pow(a, (1L<<r)*m, n) if test == (n-1): ok = 1 # 2nd test ok break else: ok = 1 # 1st test ok if ok==1: p *= 0.25 if timer: print "Elapsed time: %s seconds" % (time.clock()-start) if ok: return 1 - p else: return 0 Testing:
from mathstuff.pollard2 import millrab bignum = 57*(1L<<502)+1 millrab(bignum,20,timer=True) Elapsed time: 105.602022883 seconds 0.99999999999909051
That's down from 508 seconds. Kirby