## Daniel Shanks’ Square Form Factorization Algorithm

### August 24, 2010

We begin with the definition of ρ; since the original definition is rather complicated, we use a computationally-simpler derivation by James Pate Williams Jr.

`(define (rho f)`

(let* ((a (car f)) (b (cadr f)) (c (caddr f))

(d (- (* b b) (* 4 a c))) (s (isqrt d)) (l (abs c))

(r (if (< l (quotient s 2))

(- (* 2 l (quotient (+ b s) (* 2 l))) b)

(- (* 2 l) b))))

(list c r (quotient (- (* r r) d) (* 4 c)))))

Then `reduce`

just cycles until the input form is reduced:

`(define (reduce f)`

(let* ((a (car f)) (b (cadr f)) (c (caddr f))

(d (- (* b b) (* 4 a c))) (s (isqrt d)))

(if (< (abs (- s (* 2 (abs a)))) b s) f (reduce (rho f)))))

Our version of `squfof`

gives the forward loop in `loop`

and the backward loop in `pool`

:

`(define (squfof n)`

; assumes n is odd, composite, and not a perfect square

(define (enqueue f m q l)

(let* ((x (abs (caddr f))) (g (quotient x (gcd x m))))

(if (and (<= g l) (not (member g q))) (cons g q) q)))

(define (inv-sqrt f c)

(list (- c) (cadr f) (* -1 c (car f))))

(let* ((n4 (= (modulo n 4) 1)) (m (if n4 1 2))

(d (if n4 n (* 4 n))) (s (isqrt d))

(b (if n4 (+ (* 2 (quotient (- s 1) 2)) 1) (* 2 (quotient s 2))))

(delta (quotient (- (* b b) d) 4))

(f (list 1 b delta)) (l (isqrt s)) (bound (* 4 l)))

(let loop ((i 2) (f (rho f)) (q (enqueue f m '() l)))

(cond ((or (< bound i) (= (caddr f) 1)) #f)

((square? (caddr f)) => (lambda (c)

(if (member c q)

(loop (+ i 1) (rho f) (enqueue f m q l))

(let pool ((g (reduce (inv-sqrt f c))))

(let ((new-g (rho g)))

(if (= (cadr g) (cadr new-g))

(let ((c (abs (caddr g))))

(if (odd? c) c (/ c 2)))

(pool new-g)))))))

(else (loop (+ i 1) (rho f) (enqueue f m q l)))))))

We changed some of the variable names from the article by Jason Gower and Samuel Wagstaff that we are following; names like *d* and *D* are just too similar for comfort, and too easy to confuse. Note too that the article errs in the computation of the inverse square root, which was maddening.

Since the reduction of the square forms in the forward loop eventually cycles, Shanks’ method can fail if the cycle ends before a proper square form is found. If that happens, the easiest course is to multiply *N* by a small prime number, or a few prime numbers, and try again. For instance, consider the factorization of the twentieth repunit, R_{20} = (10^{20}−1)/9. `Squfof`

initially fails, but finds the factor 10000000001 of 11·R_{20}. Then `squfof`

fails to factor 10000000001, but finds a factor 101 of 11·10000000001, giving our first prime factor of R_{20}. And the cycle repeats: `squfof`

fails to factor 10000000001/101, but multiplying by 11 `squfof`

finds the factor 27961. That leaves a co-factor of 3541, and 101·3541·27961 = 10000000001, so we are finished with that branch of the factorization. Returning to the top, we have R_{20}/10000000001 = R_{10}, which `squfof`

fails to factor. Multiplying by 7, `squfof`

finds a factor 20867, which is divisible by 7, so the factor is 2981. Again, `squfof`

fails to factor 2981, but multiplying by 5, factoring with `squfof`

, and dividing by 5 gives us the factors 11 and 271. We’re almost finished. The remaining co-factor of R_{20} is 372731, which `squfof`

quickly factors as 41·9091. Thus, the complete factorization is R_{20} = 11 · 41 · 101 · 271 · 3541 · 9091 · 27961. Most factorizations aren’t this difficult; it is normal to remove small factors by trial division, leaving `squfof`

to do what it does best, which is to split medium-size numbers into two roughly equal factors. For instance, the factorization of the sixth fermat number, 2^{64}+1, is easily done by removing the small factors 3, 5, 17, 257 and 641, leaving a co-factor 439125228929. `Squfof`

fails to factor that number, but multiplying by 17 finds the factor 65537, and the remaining co-factor 6700417 is prime, completing the factorization.

We used isqrt from the Standard Prelude and `square?`

from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/jc1aCpdB.

Here is a version of

`squfof`

that builds in the calculations for the multiplier, given as`k`

:`(define (squfof n k)`

; assumes n is odd, composite, and not a perfect square

(define (enqueue f m q l)

(let* ((x (abs (caddr f))) (g (quotient x (gcd x m))))

(if (and (<= g l) (not (member g q))) (cons g q) q)))

(define (inv-sqrt f c)

(list (- c) (cadr f) (* -1 c (car f))))

(let* ((kn (* k n)) (n4 (= (modulo n 4) 1)) (m (if n4 k (* 2 k)))

(d (if n4 kn (* 4 k n))) (s (isqrt d))

(b (if n4 (+ (* 2 (quotient (- s 1) 2)) 1) (* 2 (quotient s 2))))

(delta (quotient (- (* b b) d) 4))

(f (list 1 b delta)) (l (isqrt s)) (bound (* 4 l)))

(let loop ((i 2) (f (rho f)) (q (enqueue f m '() l)))

(cond ((or (< bound i) (= (caddr f) 1)) #f)

((square? (caddr f)) => (lambda (c)

(if (member c q)

(loop (+ i 1) (rho f) (enqueue f m q l))

(let pool ((g (reduce (inv-sqrt f c))))

(let ((new-g (rho g)))

(if (= (cadr g) (cadr new-g))

(let* ((c (abs (caddr g)))

[…] factoring algorithms. We studied the “binary quadratic forms” variant of SQUFOF in a previous exercise; today we study the “continued fraction” variant of […]