## 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

### 4 Responses to “Williams’ P+1 Factorization Algorithm”

1. Lucas A. Brown said
```def primegen():
"""
Generates primes lazily via the sieve of Eratosthenes
Input: none
Output:
Sequence of integers
Examples:
>>> list(takewhile(lambda x: x < 100, primegen()))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
"""
yield 2; yield 3; yield 5; yield 7; yield 11; yield 13
ps = primegen() # yay recursion
p = ps.next() and ps.next()
q, sieve, n = p**2, {}, 13
while True:
if n not in sieve:
if n < q: yield n
else:
next, step = q + 2*p, 2*p
while next in sieve: next += step
sieve[next] = step
p = ps.next()
q = p**2
else:
step = sieve.pop(n)
next = n + step
while next in sieve: next += step
sieve[next] = step
n += 2

def mlucas(v, a, n):
""" Helper function for williams_pp1().  Multiplies along a Lucas sequence modulo n. """
v1, v2 = v, (v**2 - 2) % n
for bit in bin(a)[3:]: v1, v2 = ((v1**2 - 2) % n, (v1*v2 - v) % n) if bit == "0" else ((v1*v2 - v) % n, (v2**2 - 2) % n)
return v1

def williams_pp1(n):
"""
Williams' p+1 integer factoring algorithm
Input:
n -- integer to factor
Output:
Integer.  A nontrivial factor of n.
Example:
>>> williams_pp1(315951348188966255352482641444979927)
12403590655726899403L
"""
for v in count(1):
for p in primegen():
e = ilog(sqrt(n), p)
if e == 0: break
for _ in xrange(e): v = mlucas(v, p, n)
g = gcd(v - 2, n)
if 1 < g < n: return g
if g == n: break
```
2. Ruslan Khamidullin said

Hello, what is the ilog function?

3. programmingpraxis said

@Ruslan: The integer logarithm (ilog) of n to the base b is the largest integer e such that ben. You can see an implementation by clicking to run the code.

4. Pierre said

Hi, when I run this code Python says” name ‘count’ is not defined”, where is the mistake?