The Fibonacci Sequence
Tags: picolisp, lisp, algorithms
Published: 2021-12-30 10:25:12

### 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](https://picolisp-blog.hashnode.dev/series/intermediate)" 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](https://software-lab.de/doc/index.html):

> (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)
```
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)) ) ) )
```
Each time the ``do``-loop is finished, ``A`` and ``B`` have both iterated to the next higher Fibonacci-Number.

---------------------------

### The Recursive Approach

![ozymandias.png](https://cdn.hashnode.com/res/hashnode/image/upload/v1631200419634/RNUOzlpfx.png)

We have discussed the recursive solution extensively in the [Binary Tree](https://picolisp-blog.hashnode.dev/picolisp-explored-the-cache-function) 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](https://picolisp-blog.hashnode.dev/picolisp-explored-the-cache-function)**. 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](https://picolisp-blog.hashnode.dev/picolisp-explored-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](https://en.wikipedia.org/wiki/Fibonacci_number#Binet's_formula) to calculate the famous [golden ratio](https://en.wikipedia.org/wiki/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](https://simple.wikipedia.org/wiki/Fibonacci_number). (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](https://cdn.hashnode.com/res/hashnode/image/upload/v1631202892019/1kIGoFZIX.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](https://picolisp-blog.hashnode.dev/floating-point-vs-fixed-point-arithmetics) 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 ) ) )
```

*The ` in the 6th line of the ``analytic_fibonacci`` function is a so-called [Read-Macro](https://software-lab.de/doc/ref.html#macro-io) 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](https://cdn.hashnode.com/res/hashnode/image/upload/v1632920160527/EfoL4neKQ.png)

---------------------------

The source code to all tasks can be downloaded [here](https://gitlab.com/picolisp-blog/single-plage-scripts/-/tree/main/rosetta/fibonacci).

--------------------------

### Sources
- https://xkcd.com/1557/  
- https://xkcd.com/2205/  
- https://rosettacode.org/wiki/Fibonacci_sequence  
- https://en.wikipedia.org/wiki/Fibonacci_number#Binet's_formula