DEV Community

Cover image for The Fibonacci Sequence
Mia
Mia

Posted on • Originally published at picolisp-explored.com

The Fibonacci Sequence

Task description

The Fibonacci sequence is a sequence Fn of natural numbers defined recursively:

 F0 = 0 
 F1 = 1 
 Fn = Fn-1 + Fn-2, if n>1 

Write a function to generate the nth Fibonacci number.

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 var to any, and return the previous value.

Let's start with A=0, B=1. 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 A and B with each other, right?

So this is what we will do: (swap 'A B) sets A to 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)
Enter fullscreen mode Exit fullscreen mode

Now we only need to swap B with (+ A B)and we're done:

(de fib (N)
   (let (A 0  B 1)
      (do N
         (swap 'B (+ (swap 'A B) B)) ) ) )
Enter fullscreen mode Exit fullscreen mode

Each time the do-loop is finished, A and B have both iterated to the next higher Fibonacci-Number.


The Recursive Approach

ozymandias.png

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))) ) )
Enter fullscreen mode Exit fullscreen mode

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))) ) ) )
Enter fullscreen mode Exit fullscreen mode

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))) ) ) ) ) )
Enter fullscreen mode Exit fullscreen mode

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:

binetsformula.png


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 scl function (read this article if you don't know about scl).
  • we need to load the libary math.l (included in the standard pil21-distribution) in order to be able to use the function pow.

We call the first term inside the brackets P and the second one Q:

(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 ) ) )
Enter fullscreen mode Exit fullscreen mode

The in the 6th line of theanalytic_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 43349443*7* we get 43349442*8*. 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 Q=0.618 and 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 scl 100.

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 scl.
  • The upper limit has vanished: We can even calculate the 1.000.000 Fibonacci number, which has an incredible number of 208989 digits!

A Comparison

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 cache: 10.000
  • Recursive with enum: 15.000
  • 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 fixPow function).

Nevertheless, since we probably want the results to be exact, the iterative solution wins.

pi.png


The source code to all tasks can be downloaded here.


Sources

Top comments (0)