brainwagon

Sometimes, your interests converge. Over on Programming Praxis, he had a coding challenge to implement Monte Carlo factorization. The last couple of days, I have been thinking a bit more about factorization, so it seemed like a good idea to give it a whack. I cribbed a version of the Rabin Miller primality tester off the web, and then coded the rest up in python.

#!/usr/bin/env python

# 
# a basic implementation of the Pollard rho factorization
# Written by Mark VandeWettering <markv@pixar.com>
#

import sys
import locale
import random

class FactorError(Exception):
    def __init__(self, value):
	self.value = value 
    def __str__(self):
	return repr(self.value)

def miller_rabin_pass(a, n):
    d = n - 1
    s = 0
    while d % 2 == 0:
        d >>= 1
        s += 1

    a_to_power = pow(a, d, n)
    if a_to_power == 1:
        return True
    for i in xrange(s-1):
        if a_to_power == n - 1:
            return True
        a_to_power = (a_to_power * a_to_power) % n
    return a_to_power == n - 1

def isprime(n):
    for repeat in xrange(20):
        a = 0
        while a == 0:
            a = random.randrange(n)
        if not miller_rabin_pass(a, n):
            return False
    return True



def gcd(a, b):
    while b != 0:
	a, b = b, a%b
    return a 

def findfactor(n):
    for c in range(1, 50):
	x = y = random.randint(1, n-1)
	x = (x * x + c) % n
	y = (y * y + c) % n
	y = (y * y + c) % n
	while True:
	    t = gcd(n, abs(x-y))
	    if t == 1:
		x = (x * x + c) % n
		y = (y * y + c) % n
		y = (y * y + c) % n
	    elif t == n:
		break
	    else:
		return t
    raise FactorError("couldn't find a factor.")
	
def factor(n):
    r = []
    while True:
	if isprime(n):
	    r.append(n)
	    break
	try:
	    f = findfactor(n)
	    r.append(f)
	    n = n / f
	except FactorError, msg:
	    r.append(n)
	    break
    r.sort()
    return r 

def doit(n):
	flist = factor(n)
	print locale.format("%d", n, True), "="
	for f in flist:
		print "\t%s" % locale.format("%d", f, True)

locale.setlocale(locale.LC_ALL, "")

doit(2**98-1)

This is a fairly challenging factorization, and after a few seconds, it outputs:

316,912,650,057,057,350,374,175,801,343 =
	3
	43
	127
	4,363,953,127,297
	4,432,676,798,593