## Monte Carlo Factorization

### June 19, 2009

Though the explanation is lengthy, the code is agreeably short:

`(define (factor n . c)`

(define (f x c) (modulo (+ (* x x) c) n))

(let ((c (if (pair? c) (car c) 1)))

(let loop ((x 2) (y 2) (d 1))

(cond ((= d 1)

(let ((x (f x c)) (y (f (f y c) c)))

(loop x y (gcd (- x y) n))))

((= d n) (factor n (+ c 1)))

(else d)))))

`Factor`

finds a single factor of *n*. The pseudo-random sequence is generated by *f*. The main loop starts with *x* = *y* = 2; it loops if the greatest common divisor *d* is 1, restarts with the next *c* if *d* is *n*, and otherwise reports *d* as a factor of *n*.

Pollard’s method works only if *n* is composite, and the `factor`

function finds only a single factor. `Factors`

, shown below, finds all the factors of *n*; `factors`

is recursive, stopping when *n* is prime:

`(define (factors n)`

(sort < (let fact ((n n) (fs '()))

(cond ((= n 1) fs)

((even? n) (fact (/ n 2) (cons 2 fs)))

((prime? n) (cons n fs))

(else (let ((f (factor n)))

(append fs (fact f '()) (fact (/ n f) '()))))))))

`Factors`

uses the `prime?`

function, and its companion `check?`

, from the exercise on the Rabin-Miller primality checker. The factors of 2^{98} – 1 = 316912650057057350374175801343 are 3, 43, 127, 4363953127297 and 4432676798593. You can run the program at http://programmingpraxis.codepad.org/4mQExWYY.

Pages: 1 2

[…] Praxis – Monte Carlo factorization By Remco Niemeijer In today’s Programming Praxis problem we have to implement John Pollard’s factorization algorithm. Our […]

My Haskell solution (see http://bonsaicode.wordpress.com/2009/06/19/programming-praxis-monte-carlo-factorization/ for a version with comments):

Here’s my attempt in Python. A couple of issues in the code remain. The factors that it discovers aren’t guaranteed to be prime. I cribbed the Miller-Rabin test from one of the python code repositories. And, I don’t really understand exactly how this works. :-) Back to the reference books.

Okay, I fixed a couple of things, and extended the program a tiny bit. It now is a numeric calculator of sorts. It’s not industrial strength or anything, but you can basically type any python numeric expression, and it will use eval() (at least with a predefined environment) to evaluate the number. I’ve also predefined a couple of built in functions. prime(n) will return an n digit prime. rsa(n) will return an rsa key which is the combination of two n/2 digit primes. factor(n) factors n. I’ve also added code to do some trial division as well, to get rid of small factors, and it collapses multiple occurrences of a factor (instead of printing 128 copies of 2 when factoring 2^128, it outputs “2**128”).

Instead of eval, you might want to build your own calculator. See the very first Programming Praxis exercise for an RPN calculator.