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