The Fibonacci sequence is a sequence
Fnof natural numbers defined recursively:
F0 = 0 F1 = 1 Fn = Fn-1 + Fn-2, if n>1
Write a function to generate the
Solutions can be iterative or recursive (though recursive solutions are generally considered too slow and are mostly used as an exercise in recursion).
In the following, we will discuss iterative, recursive and numeric implementations and compare them in terms of speed. We will use several functions that we showed this week in the "PicoLisp Explored" series too.
The Iterative Approach
Let's start with the iterative approach that we would probably also use if we calculated it by hand. The implementation is straightforward, the only interesting point is the
swap function. From the documentation:
(swap 'var 'any) -> any
Set the value of
any, and return the previous value.
Let's start with
A+B will be our new
B, and the "old"
B will be
A. This is equivalent to
A=A+B and then simply swapping the values of
B with each other, right?
So this is what we will do:
(swap 'A B) sets
B, but returns the old value, which is the previous
A. This means that those two expressions are equivalent:
: (+ (swap 'A B) B ) # is same as: : (+ A B)
Now we only need to swap
(+ A B)and we're done:
(de fib (N) (let (A 0 B 1) (do N (swap 'B (+ (swap 'A B) B)) ) ) )
Each time the
do-loop is finished,
B have both iterated to the next higher Fibonacci-Number.
The Recursive Approach
We have discussed the recursive solution extensively in the Binary Tree series, so let's just take the best parts and quickly go through it. The "basic" implementation is easy, but awfully slow since the number of calls grow exponentially.
(de fibo (N) (if (>= 2 N) 1 (+ (fibo (dec N)) (fibo (- N 2))) ) )
A susbstantial improvement can be reached by introducing caching. Caching means that each calculation that has already been done is stored in a binary tree. This means that every number is only calculated once instead of multiple times like in the standard version.
(de fiboCache (N) (cache '(NIL) N (if (>= 2 N) 1 (+ (fiboCache (dec N)) (fiboCache (- N 2))) ) ) )
There is even room for further improvement: Instead of the
cache function, we could take the
enum function, which also builds up a binary tree of already reached values. However, instead of calculating a key to store the result,
enum takes a number as argument for the storing position.
(de fiboEnum (N) (let E (enum '(NIL) N) (or (val E) (set E (if (>= 2 N) 1 (+ (fiboEnum (dec N)) (fiboEnum (- N 2))) ) ) ) ) )
The Numeric Approach
Interestingly, there is also a third approach which is entirely numeric. It builds up on Binet's formula to calculate the famous golden ratio. The golden ratio is famous for appearing in some patterns in nature, such as the spiral arrangement of leaves and other plant parts.
An example is the sunflower in the cover picture: "florets in spirals of 34 and 55 around the outside", according to Wikipedia. (I didn't count it, so let's just believe it).
Anyway, this is Binet's formula to calculate the n'th Fibonacci number:
We can implement it in PicoLisp, but we need to take care of a few things:
- The square root of 5 is an irregular number, therefore we need a high precision in order to get correct results. This we can do using the
sclfunction (read this article if you don't know about
- we need to load the libary
math.l(included in the standard pil21-distribution) in order to be able to use the function
We call the first term inside the brackets
P and the second one
(scl 9) (load "@lib/math.l") (de analytic_fibonacci (N) (setq N (* N 1.0)) (let (Sqrt5 (sqrt 5.0 1.0) P (/ (+ 1.0 Sqrt5) 2) Q (*/ `(* 1.0 1.0) P) ) (/ (+ (*/ (+ (pow P N) (pow Q N)) 1.0 Sqrt5 ) 0.5 ) 1.0 ) ) )
The ` in the 6th line of the
analytic_fibonacci function is a so-called Read-Macro that causes the reader to evaluate the expression.
However, when we run the code, we still see a couple of problems. Things start to get wrong at
N=43: instead of 433494437 we get 433494428. The error is further increasing with higher numbers.
A second problem is that for all numbers
N>47, the results are 0!!
What is happening? The
@lib/math.l is calling a native C-library to calculate the exponent Q^N (the
pow function). However, C is calculating in floating point numbers. Analysis shows that the Binet's equation is not solvable with floating point numbers due to several issues:
Tracing reveals that
N=48, and since
Q<1, the whole term
Q^N is getting closer and closer to
0 until we get an underflow. On the other hand, we know that
P>1 which means that
P^N is growing larger and larger, until the term
P - Q is smaller than the mantisse of the floating number.
So we either get an underflow error, or we get 0 - in either way, it's a problem. This leaves only one solution: We must get rid of the
pow function that forces us to use floating point numbers.
Let's try it with integers!
Now we're getting quite deep into maths! The first problem we're facing is that the results are wrong since our representation of
(sqrt 5) was too unprecise. So let's increase it to
Secondly, we need to replace our power function
pow by an integer one, because we know now that floating numbers won't work. So, let's re-define the
pow function by a new exponential function using only integers:
(de fixPow (X N) # N th power of X (if (ge0 N) (let Y 1.0 (loop (when (bit? 1 N) (setq Y (*/ Y X 1.0)) ) (T (=0 (setq N (>> 1 N))) Y ) (setq X (*/ X X 1.0)) ) ) 0 ) )
(credits to Alexander Burger)
Now let's try again! The test reveals two things:
- The numbers are correct until
N=470. In order to further increase the precision, we simply would need to further increase the scale
- The upper limit has vanished: We can even calculate the 1.000.000 Fibonacci number, which has an incredible number of 208989 digits!
Let's compare the six approaches that we have seen in this post. Which is the highest
N we can get within 0.1 sec (depending on the hardware, obviously)?
- Recursive: 30
- Numeric, Float: 42 (only correct results)
- Numeric, Integer: 470 (only correct results)
- Recursive with
- Recursive with
- Iterative: 30.000
- Numeric, Integer: 400.000 (error unknown)
As you can see, the numeric solution is by far the fastest one (though the execution speed is also decreasing with higher
N due to the loop in the
Nevertheless, since we probably want the results to be exact, the iterative solution wins.
The source code to all tasks can be downloaded here.