## Williams’ *P*+1 Factorization Algorithm

### June 4, 2010

The code is far simpler to write than to describe. Here is the basic arithmetic over the Lucas sequences:

`(define (add a b a-b n) (modulo (- (* a b) a-b) n))`

`(define (double a n) (modulo (- (* a a) 2) n))`

`(define (mult k v n)`

(let loop ((ks (cdr (digits k 2))) (x v) (y (double v n)))

(cond ((null? ks) x)

((odd? (car ks)) (loop (cdr ks) (add x y v n) (double y n)))

(else (loop (cdr ks) (double x n) (add x y v n))))))

Then the first stage of the *p*+1 algorithm is the same as the *p*-1 algorithm, with *v* the random starting point, but the second stage changes somewhat:

`(define (pplus1-factor n b1 b2 v)`

(let stage1 ((p 2) (v v))

(let ((a (ilog p b1)))

(if (< p b1) (stage1 (next-prime p) (mult (expt p a) v n))

(let ((g (gcd (- v 2) n))) (if (< 1 g n) g

(let* ((vr v) (c (+ (* (quotient b1 6) 6) 12))

(v6 (mult 6 vr n)) (v12 (mult 12 vr n))

(z (modulo (* (- v6 vr) (- v12 vr)) n)))

(let stage2 ((c c) (x v6) (y v12) (z z))

(if (< c b2)

(let ((x+y (add y v6 x n)))

(stage2 (+ c 6) y x+y (modulo (* z (- x+y vr)) n)))

(let ((g (gcd z n))) (if (< 1 g n) g #f)))))))))))

We used `ilog`

and `digits`

from the Standard Prelude and `next-prime`

and `prime?`

from previous exercises. You can run the program at http://programmingpraxis.codepad.org/gAMB9T4D. Here’s an example:

`> (pplus1-factor 451889 10 50 7)`

139

If you wish, you can add Williams’ *p*+1 method to the integer factorization program of a previous exercise. A one-stage version of the algorithm, parameters for the `factors`

function, and the calling code are shown below. The *p*+1 method fits best between the *p*-1 method and the elliptic curve method:

`(define (pplus1-factor n b v)`

(define (m x) (modulo x n))

(define (add a b a-b) (m (- (* a b) a-b)))

(define (double a) (m (- (* a a) 2)))

(define (mult k v)

(let loop ((ks (cdr (digits k 2))) (x v) (y (double v)))

(cond ((null? ks) x)

((odd? (car ks)) (loop (cdr ks) (add x y v) (double y)))

(else (loop (cdr ks) (double x) (add x y v))))))

(let loop ((p 2) (v v) (k 0))

(cond ((< b p) (let ((g (gcd (- v 2) n))) (if (< 1 g n) g #f)))

((zero? (modulo k 10 0)) (let ((g (gcd (- v 2) n))) (if (< 1 g n) g

(loop (next-prime p) (mult (expt p (ilog p b)) v) (+ k 1)))))

(else (loop (next-prime p) (mult (expt p (ilog p b)) v) (+ k 1))))))

` (define pplus1-limit 100000) ; iteration limit`

(define pplus1-trials 7) ; number of p+1 constants to try

` ; williams pplus1`

(let loop ((k pplus1-trials) (v (randint 3 n)))

(when (positive? k)

(msg "Williams p+1: bound=" pplus1-limit ", constant=" v)

(if (factor? 'p+1 (pplus1-factor n pplus1-limit v))

(loop k (randint n))

(loop (- k 1) (randint n)))))

Pages: 1 2