Xref: utzoo sci.math:16709 comp.lang.scheme:2317
Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!pacific.mps.ohio-state.edu!linac!att!news.cs.indiana.edu!ariel.unm.edu!cie!scavo
From: scavo@cie.uoregon.edu (Tom Scavo)
Newsgroups: sci.math,comp.lang.scheme
Subject: Re: fibonacci Numbers
Summary: an O(log n) scheme algorithm is given
Keywords: Fibonacci Numbers - routines?
Message-ID: <1991Apr14.030812.887@ariel.unm.edu>
Date: 14 Apr 91 03:08:12 GMT
References: <1991Apr12.051844.15063@milton.u.washington.edu>
Reply-To: scavo@cie.uoregon.edu (Tom Scavo)
Followup-To: comp.lang.scheme
Organization: University of Oregon Campus Information Exchange
Lines: 129

In article <1991Apr12.051844.15063@milton.u.washington.edu> amigo@milton.u.washington.edu (The Friend) writes:
>
>     Can someone send me source code for a routine to calculate 
>Fibonacci numbers to _x_ (x being input from user)? Additionally it'd
>be great if this found the PRIME numbers in the set. It could print its
>results as it went through (if that makes it any easier).


;*************************************************************************
; The recursive definition of the Fibonacci sequence yields an O(k^n)
; algorithm, where k is the golden ratio = (/ (+ 1 (sqrt 5)) 2).

(define (fib1 n)
  (if (< n 2)
      n
      (+ (fib1 (- n 1))
         (fib1 (- n 2)))))


;*************************************************************************
; An auxiliary procedure that carries along the previous Fibonacci
; number helps to eliminate redundant calculations and produce O(n)
; behavior.

(define (fib2 n)
  (define (loop i fi fi-1)
    (if (= i n)
        fi
        (loop (+ i 1) (+ fi fi-1) fi)))
  (cond ((< n 2) n)
        ((exact? n) (loop 1 1 0))
        (else (loop 1 1.0 0.0))))


;*************************************************************************
; This procedure implements the closed form solution to the Fibonacci
; difference equation.  Its efficiency hinges on the unknown efficiency
; of some primitive numerical routines.

(define (fib3 n)
  (let ((sqrt5 (sqrt 5.0)))
    (round (/ (expt (/ ( + 1.0 sqrt5) 2.0) n)
              sqrt5))))


;*************************************************************************
; This O(log n) algorithm is an extension of the corresponding algorithm
; for exponentiation of numbers found in just about any introductory
; textbook.  Here we take advantage of the fact (given to me by Will
; Clinger) that

;
;                                n
;                      | 1   1 |        | fn+1  fn   |
;                      |       |    =   |            |
;                      | 1   0 |        | fn    fn-1 |


; and exponentiate the 2x2 matrix by successive squaring.


(define realIdentity '((1.0 0.0)
                       (0.0 1.0)))

(define integerIdentity '((1 0)
                          (0 1)))

(define (extract M)
  (caar M))


; This top level procedure determines the type of its lone argument.

(define (fib4 n)
  (extract (if (exact? n)
               (power (- n 1) integeridentity)
               (power (- n 1) realidentity))))


; The actual calculations take place here.  If the exponent is even,
; we square the matrix and halve the exponent since x^n = (x^2)^(n/2)
; when n is even.  Otherwise, we factor out an x before squaring.

(define power
  (lambda (exponent identity)
    (define (exponential i)
      (cond ((= i 0) identity)
            ((odd? i) (product (square (exponential (quotient (- i 1)
                                                              2)))))
            (else (square (exponential (quotient i 2))))))
    (exponential exponent)))


; This special purpose routine multiplies the 2x2 matrix given at the 
; outset of this discussion by its argument.  It also takes advantage
; of the fact that x is symmetric.

(define (product x)
  (let* ((a11 (caar x))
         (a12 (cadar x))
         (a21 a12)                             ;in general, (caadr x)
         (a22 (cadadr x)))
    (list (list (+ a11 a21)
                (+ a12 a22))
          (list a11
                a12))))


; This procedure is a bit more general in that it will square and return
; any 2x2 symmetric matrix.  All potentially redundant calculations are
; being performed once and for all.

(define (square x)
  (let* ((a11 (caar x))
         (a12 (cadar x))
         (a21 a12)                             ;in general, (caadr x)
         (a22 (cadadr x))
         (trace (+ a11 a22))
         (crossproduct (* a12 a21))
         (aij (* a12 trace)))                  ;or (* a21 trace)
    (list (list (+ crossproduct
                   (* a11 a11))
                aij)
          (list aij
                (+ crossproduct
                   (* a22 a22))))))

Tom Scavo
scavo@cie.uoregon.edu