Programming Praxis

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₂₀ = (10²⁰−1)/9. Squfof initially fails, but finds the factor 10000000001 of 11·R₂₀. Then squfof fails to factor 10000000001, but finds a factor 101 of 11·10000000001, giving our first prime factor of R₂₀. 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₂₀/10000000001 = R₁₀, 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₂₀ is 372731, which squfof quickly factors as 41·9091. Thus, the complete factorization is R₂₀ = 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⁶⁴+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.

Pages: 1 2