How Fermat Factored Integers
July 18, 2014
We begin by calculating the table of squares. Instead of printing it as we did on the previous page, we store it in a vector, plus one additional square to act as a sentinal:
(define limit 10000)
(define squares
(let ((squares (make-vector (+ limit 1))))
(do ((i 0 (+ i 1))) ((< limit i) squares)
(vector-set! squares i (* i i)))))
Then isqrt, which calculates the largest number that, when multiplied by itself, doesn’t exceed n, and in-table, which returns the square root of n if n is a square or #f otherwise, are simple:
(define (isqrt n)
(do ((i 1 (+ i 1))) ((< n (vector-ref squares i)) (- i 1)))))
(define (in-table n)
(let loop ((i 0))
(cond ((< n (vector-ref squares i)) #f)
((= n (vector-ref squares i)) i)
(else (loop (+ i 1))))))
We calculate the digital root algorithmically; Fermat would certainly have calculated it visually by casting out nines:
(define (digital-root n)
(let loop ((n n) (r 0))
(cond ((zero? n) (if (< r 10) r (digital-root r)))
((< n 10) (loop 0 (+ r n)))
(else (let ((d (modulo n 10)))
(loop (/ (- n d) 10) (+ r d)))))))
Here we encode the rules for identifying squares. The square? function returns the square root if it exists or #f otherwise:
(define (square? n)
(let* ((ones (modulo n 10))
(tens (modulo (/ (- n ones) 10) 10)))
(and (or (and (= 0 ones) (zero? tens))
(and (= 1 ones) (even? tens))
(and (= 4 ones) (even? tens))
(and (= 5 ones) (= 2 tens))
(and (= 6 ones) (odd? tens))
(and (= 9 ones) (even? tens)))
(member (digital-root n) '(1 4 7 9))
(in-table n))))
We finish with the function that simulates Fermat’s calculation. First we ensure n fits in our table and remove factors of 2. Then we initialize r and t, iterate until we find a square, and calculate the factors:
(define (fermat n)
(if (not (< -1 n (* limit limit))) #f
(if (even? n) (list 2 (/ n 2))
(let ((x (isqrt n)))
(if (= (* x x) n) (list x x)
(let loop ((r (- (* x x) n)) (t (+ x x 1)))
(if (not (square? r))
(loop (+ r t) (+ t 2))
(let ((x (/ (- t 1) 2)) (y (isqrt r)))
(list (- x y) (+ x y))))))))))
You can run the program at http://programmingpraxis.codepad.org/NtvkBMaG:
