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, R20 = (1020−1)/9. Squfof initially fails, but finds the factor 10000000001 of 11·R20. Then squfof fails to factor 10000000001, but finds a factor 101 of 11·10000000001, giving our first prime factor of R20. 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 R20/10000000001 = R10, 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 R20 is 372731, which squfof quickly factors as 41·9091. Thus, the complete factorization is R20 = 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, 264+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.

About these ads

Pages: 1 2

2 Responses to “Daniel Shanks’ Square Form Factorization Algorithm”

  1. programmingpraxis said

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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 600 other followers

%d bloggers like this: