Programming Praxis

We begin with a function that can compute a Lucas sequence, using the recursion X_n = P X_n-1 − Q X_n-2:

(define (lseq p q x0 x1 n)
  (let loop ((x0 x0) (x1 x1) (k 1) (xs (list x1 x0)))
    (if (= k n) (reverse xs)
      (let ((x2 (- (* p x1) (* q x0))))
        (loop x1 x2 (+ k 1) (cons x2 xs))))))

Here are three familiar examples:

> (lseq 1 -1 0 1 20) ; fibonacci numbers
(0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765)
> (lseq 1 -1 2 1 20) ; lucas numbers
(2 1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364 2207 3571 5778 9349 15127)
> (lseq 3 2 0 1 20) ; mersenne numbers 2^n - 1
(0 1 3 7 15 31 63 127 255 511 1023 2047 4095 8191 16383 32767 65535 131071 262143 524287 1048575)

We look next at the function that computes the nth element of a Lucas sequence left-to-right. Since the U sequence depends on the V sequence, it is easiest to compute both at the same time, and return both values. We add two features that we will frequently need; we return Qⁿ in addition to U_n and V_n, and we take a modulus m, which can be disabled by setting m = 0. The comments show how the chain is built depending on the binary digits of n:

(define (lucas p q m n) ; left-to-right
  (define (mod n) (if (zero? m) n (modulo n m)))
  (let bits ((n n) (bs (list)))
    (if (positive? n) (bits (quotient n 2) (cons (modulo n 2) bs))
      (let loop ((un 0) (uu 1) (vn 2) (vv p) (qn 1) (bs bs))
        (cond ((null? bs) (values un vn qn))
              ((odd? (car bs)) ; 1-bit
                (loop (mod (- (* uu vn) qn)) ; u(2n+1)
                      (mod (* uu vv)) ; u(2n+2)
                      (mod (- (* vv vn) (* p qn))) ; v(2n+1)
                      (mod (- (* vv vv) (* 2 q qn))) ; v(2n+2)
                      (mod (* qn qn q)) (cdr bs)))
              (else ; 0-bit
                (loop (mod (* un vn)) ; u(2n)
                      (mod (- (* uu vn) qn)) ; u(2n+1)
                      (mod (- (* vn vn) (* 2 qn))) ; v(2n)
                      (mod (- (* vv vn) (* p qn))) ; v(2n+1)
                      (mod (* qn qn)) (cdr bs))))))))

We will find it useful for testing to define utility functions that compute the nth elements of the U and V sequences:

(define (u p q n)
  (call-with-values
    (lambda () (lucas p q 0 n))
    (lambda (un vn qn) un)))

(define (v p q n)
  (call-with-values
    (lambda () (lucas p q 0 n))
    (lambda (un vn qn) vn)))

We have two ways to test our functions. The first is to compare the two functions against each other: compute a variety of Lucas sequences, and for each compute every element of the sequence, comparing the results. The second is to compare our computations to an outside oracle.

To compare our two functions against each other, we compute and store the first 101 elements of each sequence for -10 ≤ P ≤ 10 and -10 ≤ Q ≤ 10:

(define (test1)
  (do ((p -10 (+ p 1))) ((< 10 p))
    (do ((q -10 (+ q 1))) ((< 10 q))
      (display "P=") (display p)
      (display ", Q=") (display q) (newline)
      (assert (lseq p q 0 1 100) (map (lambda (n) (u p q n)) (range 101)))
      (assert (lseq p q 2 p 100) (map (lambda (n) (v p q n)) (range 101))))))

Running the test shows the 21 × 21 = 441 P and Q parameters being tested for both the U and V sequences:

> (test1)
P=-10, Q=-10
P=-10, Q=-9
P=-10, Q=-8
P=-10, Q=-7
P=-10, Q=-6
… 431 lines elided …
P=10, Q=6
P=10, Q=7
P=10, Q=8
P=10, Q=9
P=10, Q=10

For an outside oracle we scrape Lucas sequences from the supremely valuable On-Line Encyclopedia of Integer Sequences. There are 45 Lucas sequences in the OEIS database, which we enumerate in an easy-to-parse format:

(define oeis-lucas '((-1  3 u "A214733") (1 -1 u "A000045") (1 -1 v "A000032")
  (1  1 u "A128834") ( 1  1 v "A087204") (1  2 u "A107920") (2 -1 u "A000129")
  (2 -1 v "A002203") ( 2  1 u "A001477") (2  2 u "A009545") (2  2 v "A007395")
  (2  3 u "A088137") ( 2  4 u "A088138") (2  5 u "A045873") (3 -5 u "A015523")
  (3 -5 v "A072263") ( 3 -4 u "A015521") (3 -4 v "A201455") (3 -3 u "A030195")
  (3 -3 v "A172012") ( 3 -2 v "A206776") (3 -1 u "A006190") (3 -1 v "A006497")
  (3  1 u "A001906") ( 3  1 v "A005248") (3  2 u "A000225") (3  2 v "A000051")
  (3  5 u "A190959") ( 4 -3 u "A015530") (4 -3 v "A080042") (4 -2 u "A090017")
  (4 -1 u "A001076") ( 4 -1 v "A014448") (4  1 u "A001353") (4  1 v "A003500")
  (4  2 v "A056236") ( 4  3 u "A003462") (4  3 v "A034472") (4  4 u "A001787")
  (5 -3 u "A015536") ( 5 -2 u "A015535") (5 -1 v "A087130") (5  1 v "A003501")
  (5  4 u "A002450") ( 5  4 v "A052539")))

The sequences of the OEIS are available in an internal format that we can process easily; the only problem is that most sequences use the %S, %T and %U fields from the database, but if there are negative numbers within the sequence, the %V, %W and %X fields must be used:

(define (get-oeis seq)
  (with-input-from-url
    (string-append "http://oeis.org/search?q=id:" seq "&fmt=text")
    (lambda ()
      (let ((seq-str ""))
        (do ((line (read-line) (read-line)))
            ((eof-object? line) (map string->number (string-split #\, seq-str)))
          (when (< 2 (string-length line))
            (case (substring line 0 2)
              (("%S" "%T" "%U")
                (let ((fields (string-split #\space line)))
                  (set! seq-str (string-append seq-str (caddr fields)))))
              (("%V")
                (let ((fields (string-split #\space line)))
                  (set! seq-str (caddr fields))))
              (("%W" "%X")
                (let ((fields (string-split #\space line)))
                  (set! seq-str (string-append seq-str (caddr fields))))))))))))

The check-lucas function compares a single sequence to the OEIS, and the test2 function compares all the available sequences:

(define (check-lucas p q u-or-v seq)
  (display seq) (newline)
  (let loop ((ss (get-oeis seq)) (n 0))
    (when (pair? ss)
      (if (eq? u-or-v 'u)
          (assert (u p q n) (car ss))
          (assert (v p q n) (car ss)))
      (loop (cdr ss) (+ n 1)))))

(define (test2)
  (for-each (lambda (t) (apply check-lucas t)) oeis-lucas))

And here we perform the test:

>(test2)
A214733
A000045
A000032
A128834
A087204
… 35 lines elided …
A015535
A087130
A003501
A002450
A052539

We can do all that same testing on this version of the function that works from right to left:

(define (lucas p q m n) ; right-to-left
  (define (even e o) (if (even? n) e o))
  (define (mod n) (if (zero? m) n (modulo n m)))
  (let ((d (- (* p p) (* 4 q))))
    (let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
               (u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
      (if (zero? n) (values u v k)
        (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
              (q2 (mod (* qn qn))) (n2 (quotient n 2)))
          (if (even? n) (loop u2 v2 q2 n2 u v k)
            (let* ((uu (+ (* u v2) (* u2 v)))
                   (vv (+ (* v v2) (* d u u2)))
                   (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
                   (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
                   (uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
              (loop u2 v2 q2 n2 uu vv (* k q2)))))))))

I’ll leave it to you to run that function through the same tests as the other function.

I’m convinced now that the functions for computing Lucas sequences work properly, and hope that I’ve convinced you, my readers. Although we didn’t mention it above, we stole much code from the Standard Prelude and from previous exercises; you can see the entire code assembled at http://programmingpraxis.codepad.org/WfXn48J3. I’m in the process of fixing all of the previous exercises that computed Lucas sequences that I know were incorrect, fixing more bugs than I will ever admit. And, in case you haven’t guessed, we will soon have a use for our new Lucas functions.

Pages: 1 2