; leonardo numbers
(define (make-matrix rows columns . value)
(do ((m (make-vector rows)) (i 0 (+ i 1)))
((= i rows) m)
(if (null? value)
(vector-set! m i (make-vector columns))
(vector-set! m i (make-vector columns (car value))))))
(define (matrix-rows x) (vector-length x))
(define (matrix-cols x) (vector-length (vector-ref x 0)))
(define (matrix-ref m i j) (vector-ref (vector-ref m i) j))
(define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))
(define-syntax for
(syntax-rules ()
((for (var first past step) body ...)
(let ((ge? (if (< first past) >= <=)))
(do ((var first (+ var step)))
((ge? var past))
body ...)))
((for (var first past) body ...)
(let* ((f first) (p past) (s (if (< first past) 1 -1)))
(for (var f p s) body ...)))
((for (var past) body ...)
(let* ((p past)) (for (var 0 p) body ...)))))
(define (matrix-multiply a b)
(let ((ar (matrix-rows a)) (ac (matrix-cols a))
(br (matrix-rows b)) (bc (matrix-cols b)))
(if (not (= ac br))
(error 'matrix-multiply "incompatible matrices")
(let ((c (make-matrix ar bc 0)))
(for (i ar)
(for (j bc)
(for (k ac)
(matrix-set! c i j
(+ (matrix-ref c i j)
(* (matrix-ref a i k)
(matrix-ref b k j)))))))
c))))
(define (matrix-power m n)
(cond ((= n 1) m)
((even? n)
(matrix-power
(matrix-multiply m m)
(/ n 2)))
(else (matrix-multiply m
(matrix-power
(matrix-multiply m m)
(/ (- n 1) 2))))))
(define (leo1 n)
(if (< n 2) 1
(+ (leo1 (- n 2)) (leo1 (- n 1)) 1)))
(display (leo1 30)) (newline)
(define (leo2 n)
(if (< n 2) 1
(let loop ((L-2 1) (L-1 1) (L 3) (k 2))
(if (= n k) L
(loop L-1 L (+ L-1 L 1) (+ k 1))))))
(display (leo2 30)) (newline)
(define (fib n)
(if (zero? n) 0
(matrix-ref (matrix-power #(#(1 1) #(1 0)) n) 1 0)))
(define (leo3 n) (- (* (fib (+ n 1)) 2) 1))
(display (leo3 30)) (newline)