scavo@cie.uoregon.edu (Tom Scavo) (04/14/91)
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
haltraet@gondle.idt.unit.no (Hallvard Traetteberg) (04/14/91)
In a class on Knowledge Based Software Design we had an assignment on
fibonacci-numbers. We were supposed to demonstrate how one could logically
deduce the standard linear algorithm from the exponential one. When trying to
solve the problem I found the following formula:
fib(n) = fib(m) * fib(n-m) + fib(m-1) * fib (n-m-1)
It's easy to prove correct by induction and equally easy to "prove"
intuitively by expanding the standard formula fib(n) = fib(n-1) + fib(n-2) a
couple of times.
By choosing m = (n+ n mod 2)/2 one is able to divide the problem into two
approximately equal problems. Writing the equations for n mod 2 = 0 and
n mod 2 = 1 gives:
fib(n) = (fib(n/2))^2 + (fib(n/2 - 1))^2 ; n mod 2 = 0
fib(n) = fib((n-1)/2)(fib((n-1)/2) + 2*fib((n-3)/2)) ; n mod 2 = 1
Each of the equations gives two fib-calculations around n/2, so the complexity
is linear. Since the n's are divided by two and thus reach 1 after log(n)
steps, the computation tree is bound to calculate som fib-numbers more than
once. The same kind of reasoning tells us that the standard fib-formula makes
redundant calculations since the complexity is exponential and n decreases by
at least 1 in each call.
By expanding the formulas I was able to find where the redundant calculations
were made. It wasn't difficult to reformulate the algorithm so to calculate
each fib-number only once. Since we still reach 1 after log(n) steps the
algorithm is O(log(n)). Unfortunately it is written in Common Lisp and uses
multiple values, but it shouldn't be hard to translate to Scheme:
;;; Auxillary function
(defun fib2-aux (n f1 f2)
(cond ((< n 2) (values 1))
((evenp n) (+ (* f1 f1) (* f2 f2)))
(t (* f1 (+ f1 (* 2 f2))))
))
;;; The fib-function itself (fib1 was the linear version, hence the name).
(defun fib2 (n)
(case n
((0 1) (values 1 n 0))
(2 (values 2 1 1))
(t (multiple-value-bind (f1 f2 f3) (fib2 (truncate n 2))
(values (fib2-aux n f1 f2)
(if (evenp n)
(fib2-aux (1- n) f2 f3)
(fib2-aux (1- n) f1 f2))
(fib2-aux (- n 2) f2 f3))
))
))
It's a bit hard to explain, but it works! It actually calculates fib(n),
fib(n-1) and fib(n-2) at the same time, so it returns three numbers:
>(trace fib2)
>(fib2 100)
1> (FIB2 100)
2> (FIB2 50)
3> (FIB2 25)
4> (FIB2 12)
5> (FIB2 6)
6> (FIB2 3)
7> (FIB2 1)
<7 (FIB2 1 1 0)
<6 (FIB2 3 2 1)
<5 (FIB2 13 8 5)
<4 (FIB2 233 144 89)
<3 (FIB2 121393 75025 46368)
<2 (FIB2 20365011074 12586269025 7778742049)
<1 (FIB2 573147844013817084101 354224848179261915075 218922995834555169026)
573147844013817084101
354224848179261915075
218922995834555169026
The trace shows that it is logarithmic, it calculates the three fib-numbers
fib(n), fib(n-1) and fib(n-2) for each n. These are the only fib-numbers that
are calculated and additionally they are calculated only once.
The use of multiple values make the code rather nice. I guess an equivalent
scheme implementation will have to return a list and destructure it with car
and cdr. Common Lisp isn't that bad is it? :-)
- hal (Student at the Norwegian Institute of Technology)
--
- haltraet (@idt.unit.no)