[Edu-sig] More on Factoring

[Kirby Urner]

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