Programming Praxis

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:

Pages: 1 2 3