# Efficient implementations of Fibonacci function with interesting results

**URL:** <https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123>\
**Category:** Numerics\
**Created:** [November 29, 2018, 1:34am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123 "2018-11-29T01:34:52Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![mhw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mhw/32/8469_2.png) [@mhw](https://discourse.julialang.org/u/mhw)\
**Post date:** [November 29, 2018, 1:34am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/1 "2018-11-29T01:34:53Z")

</div>

It is well known that the standard recursive implementation of the Fibonacci function defined by

```julia
F(0) = 0
F(1) = 1
F(n) = F(n-1) + F(n-2) for n > 1

```

runs in exponential time. Therefore most often the Fibonacci function is implemented either iteratively

```julia
function itfib(n)
    x,y = BigInt(0), BigInt(1)
    for i in 1:n
        x,y = y, x+y
    end
    return x
end

```

or by matrix exponentation

```julia
const b = [BigInt(1) 1; 1 0]

function mfib(n)
    if n == 0
        return 0
    elseif n == 1
        return 1
    else
        return (b^(n+1))[2,2]
    end
end

```

The underlying algorithms run in linear time. Using a divide-and-conquer strategy for matrix exponentation we get an algorithm running in time `O(log n)`:

```julia
function pow(a,n) # n > 0 !
    if n == 1
        return a
    else
        c = pow(a, n >> 1)
        if iseven(n)
            return c*c
        else
            return c*c * a
        end
    end
end  

function afib(n)    
    if n == 0
        return 0
    elseif n == 1
        return 1
    else
        return (pow(b,n+1))[2,2]
    end
end

```

Another implementation running in time `O(log n)` uses for `n>1`Dijkstra’s recurrence equation

```julia
F(2n) = F(n) ( 2F(n-1) + F(n) )

F(2n+1) = F(n)^2 + F(n+1)^2

```

and memoization:

```julia
const fibmem = Dict{Int,BigInt}()
function fib(n)
    get!(fibmem, n) do
        if n <= 0 
            return BigInt(0)
        elseif n == 1 
            return BigInt(1)
        else
            m = n >> 1
            if iseven(n)
                return fib(m)*(2*fib(m-1) + fib(m))
            else
                return fib(m+1)^2 + fib(m)^2
            end
        end
    end
end

```

I tested the iterative version itfib for `10^7` and all the other version for `10^9` and got the following results

```julia
@elapsed(fib(10^9)) = 12.006108423

@elapsed(afib(10^9)) = 52.248328254

@elapsed(mfib(10^9)) = 78.922591838

@elapsed(itfib(10^7)) = 882.27293546

```

I was a little bit surprised that Dijkstra’s recurrence relation together with memoization was so much faster than the efficient matrix exponentation method (afib). I was also surprised that the standard matrix exponenation method (mfib) was only a “bit slower” than afib, but much faster than the iterative implementation. What are the reasons for these facts?

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [November 29, 2018, 1:47am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/2 "2018-11-29T01:47:22Z")

</div>

The fastest method to compute these is the Binet formula F\_n = \frac1{\sqrt5}(\Phi^n-\Phi^{-n}), where \Phi = \frac{1+\sqrt5}2.

```nohighlight
julia> F(n) = ((1+sqrt(big(5)))^n-(1-sqrt(big(5)))^n)/(sqrt(big(5))*big(2)^n)
F (generic function with 1 method)

julia> @btime F(10^7);
  167.953 μs (25 allocations: 1.19 MiB)

```

It’s slightly slower if you convert it to `BigInt` afterwards.

---

<div class="post-metadata">

**Author:** ![xiaodai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xiaodai/32/15937_2.png) [@xiaodai](https://discourse.julialang.org/u/xiaodai)\
**Post date:** [November 29, 2018, 2:05am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/3 "2018-11-29T02:05:10Z")

</div>

Well. You can express Fibonacci sequence as powering a matrix, but you can eigen-decompose the matrix into A^{-1}DA where D is diagonal, so powering that is dead easy! that’s one way to derive the Binet formula above.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 3:57am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/4 "2018-11-29T03:57:12Z")

</div>

> [@chakravala](#):
>
> The fastest method to compute these is the Binet formula

The Binet formula is useless if you want an exact integer result for large n.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 4:02am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/5 "2018-11-29T04:02:36Z")

</div>

The fastest method in Julia is probably to call the [optimized C implementation](https://github.com/alisw/GMP/blob/cc09e59afe9a23de9e7d0ef10598a395d7e8c850/mpz/fib_ui.c) from the underlying `libgmp` library:

```julia
function fastfib(n)
    z = BigInt()
    ccall((:__gmpz_fib_ui, :libgmp), Cvoid, (Ref{BigInt}, Culong), z, n)
    return z
end

```

This is about 50% faster than your memo-ized `fib` implementation above (which means your `fib` function is quite good)!

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [November 29, 2018, 11:40am UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/6 "2018-11-29T11:40:13Z")

</div>

> [@stevengj](#):
>
> The Binet formula is useless if you want an exact integer result for large n.

Actually, the formula produces an EXACT integer, if you have trouble getting one, it is due to the limitations of your floating point arithmetic.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 29, 2018, 12:30pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/7 "2018-11-29T12:30:41Z")

</div>

> [@chakravala](#):
>
> it is due to the limitations of your floating point arithmetic

I was under the impression that this was @stevengj’s point, as all floating point arithmetic has limitations of some kind, starting from the fact that the golden ratio is irrational, and it only gets worse from there.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 1:14pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/8 "2018-11-29T13:14:01Z")

</div>

> [@Tamas\_Papp](#):
>
> I was under the impression that this was @stevengj’s point, as all floating point arithmetic has limitations of some kind, starting from the fact that the golden ratio is irrational, and it only gets worse from there.

Right. You need rapidly increasing precision with `n` in order to guarantee a correct rounded result. Of course, this is possible to do with `BigFloat` via an appropriate call to `setprecision`, but that probably negates any performance advantage of this approach. Just using the default `BigFloat` precision of 256 bits (about 77 decimal digits), as @chakravala did above for `n=10^7`, gives the wrong results.

For example, I tried repeatedly doubling the `BigFloat` precision until the Binet formula gave an exact result (when rounded to `BigInt`) for `n=10^7`, and it required a precision of 8388608 bits. At this precision, the Binet `F(n)` function posted by @chakravala takes 3.2 seconds on my machine, whereas `fastfib(n)` takes 0.05 seconds.

For `n=10^9` as in the original post (for which `fastfib` takes about 8 seconds), I tried repeatedly doubling the precision to get an exact answer from `F(n)`, but I ran out of patience after waiting a few minutes.

---

<div class="post-metadata">

**Author:** ![bpr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bpr/32/2710_2.png) [@bpr](https://discourse.julialang.org/u/bpr)\
**Post date:** [November 29, 2018, 2:26pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/9 "2018-11-29T14:26:53Z")

</div>

It’s good that calling out to C is so easy, but part of Julia’s elevator pitch is that it solves the two language problem in data science. Can that 50% difference be reduced much more by a pure Julia function?

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [November 29, 2018, 2:29pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/10 "2018-11-29T14:29:43Z")

</div>

@stevengj the performance is bad because you are using naive implementations. It would probably be much quicker using `arb` with Nemo package and using the Chakravala algorithm to approximate the square root. The chakravala algorithm combined with the bhavana principle will significantly outperform the rate of convergence of regular continued fractions, and it has not even been tried to apply it to the square root of 5 here to obtain the accuracy of the roots to a higher precision in shorter time. Perhaps, I will give it a try when I have some free time, I will try to see if the chakravala algorithm and bhavana principle can be used to converge quicker than these by using them to compute square root of 5 for Binet’s formula, also using the arbitrary precision facilities of Nemo.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 2:40pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/11 "2018-11-29T14:40:36Z")

</div>

> [@bpr](#):
>
> Can that 50% difference be reduced much more by a pure Julia function?

Probably, but with some caveats. First, of course, you’d want to look at the GMP implementation to see if they are using any high-level algorithmic tricks that can easily be copied.

However, realize that the `BigInt` type in Julia is based on an external C library (GMP), not pure Julia code, so to get the maximum possible optimizations you may need to call low-level GMP functions. In particular low-level GMP code can eke out some additional performance gains by careful use of in-place operations on pre-allocated `BigInt` objects (for example, by re-using the same object for more than one calculation).

Automating those kind of gains seems like it will require lower-level support for memory management of bignum objects in the Julia language implementation, which has not been a high priority of most of the developers so far. See, for example, [#4176](https://github.com/JuliaLang/julia/issues/4176) and [#11202](https://github.com/JuliaLang/julia/issues/11202).

(If you compare the code length and clarity for the GMP C code and the Julia `fib(n)` implementation, in my opinion a 50% performance penalty is quite good for such short and readable code. In general, getting the “last factor of 2” for a highly-optimized implementation of a nontrivial algorithm often takes an unreasonable amount of code in _any_ language, and is rarely worth it except for critical library functions…you’re usually better off focusing on extensibility, generality, and clarity.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 2:47pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/12 "2018-11-29T14:47:15Z")

</div>

> [@chakravala](#):
>
> the performance is bad because you are using naive implementations. It would probably be much quicker using `arb` with Nemo package and using the Chakravala algorithm to approximate the square root.

The `sqrt` function for `BigFloat` is calling `mpfr_sqrt` from the MPFR library; if MPFR is using a “naive” algorithm, you should file an issue with MPFR, but I’m skeptical of that assertion without evidence.

---

<div class="post-metadata">

**Author:** ![bpr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bpr/32/2710_2.png) [@bpr](https://discourse.julialang.org/u/bpr)\
**Post date:** [November 29, 2018, 4:17pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/13 "2018-11-29T16:17:34Z")

</div>

> [@stevengj](#):
>
> (If you compare the code length and clarity for the GMP C code and the Julia `fib(n)` implementation, in my opinion a 50% performance penalty is quite good for such short and readable code. In general, getting the “last factor of 2” for a highly-optimized implementation of a nontrivial algorithm often takes an unreasonable amount of code in _any_ language, and is rarely worth it except for critical library functions…you’re usually better off focusing on extensibility, generality, and clarity.)

I agree, but I’m curious about whether the same “unreasonable” code can be written in Julia as it can in C or C++. Could BigInt have been written in pure Julia, given the manpower and time, or will two languages always be required to achieve the highest performance?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 4:35pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/14 "2018-11-29T16:35:57Z")

</div>

> [@bpr](#):
>
> Could BigInt have been written in pure Julia, given the manpower and time.

Yes, probably, but it would require a huge amount of effort to replicate GMP.

One tricky thing to replicate in “pure” Julia is code that uses inline assembly to access specialized CPU instructions, e.g. the `crc32c` function in Julia’s `CRC32c` standard library calls a [highly optimized C implementation](https://stackoverflow.com/questions/17645167/implementing-sse-4-2s-crc32c-in-software) that utilizes inline assembly to access the crc32c instruction of modern x86 processors. It is [possible to access this in Julia](https://github.com/andrewcooke/CRC.jl/issues/5) via the `Base.llvmcall` “escape hatch” that lets you insert LLVM instructions directly, but there was not much point in rewriting the C in this case.

Another difficulty is to replicate code that uses low-level thread management, since you can’t call these directly from Julia without worrying about the thread-safety of the Julia runtime. However, the native Julia multithreading support has been getting better by leaps and bounds, so there’s not much need for lower level access anymore.

In general, excluding libraries like GMP or OpenBLAS or FFTW that involve tens of thousands of lines of low-level code and require a huge engineering effort to replicate, I can’t recall any case where careful tuning of Julia code was not able to replicate the performance of C or Fortran within a few percent. I’ve sometimes found that one can [even exceed the speed of optimized C or Fortran](https://github.com/JuliaLang/julia/pull/2987), because Julia features like metaprogramming allow us to code in a style that would be difficult to replicate in a lower-level language.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [November 29, 2018, 4:50pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/15 "2018-11-29T16:50:41Z")

</div>

> [@stevengj](#):
>
> It is [possible to access this in Julia](https://github.com/andrewcooke/CRC.jl/issues/5) via the `Base.llvmcall` “escape hatch” that lets you insert LLVM instructions directly, but there was not much point in rewriting the C in this case.

I wrote a bit about calling intrinsics in [SIMD and SIMD-intrinsics in Julia | Kristoffer Carlsson](http://kristofferc.github.io/post/intrinsics/) (under the header “Using intrinsics”). It can be done slightly simpler than in the linked version:

```julia
julia> crc32(crc::Int32, x::Int32) = ccall("llvm.x86.sse42.crc32.32.32", llvmcall, Int32, (Int32, Int32), crc, x);

julia> crc32(Int32(2), Int32(3))
-582636872

```

---

<div class="post-metadata">

**Author:** ![jwmerrill](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jwmerrill/32/6261_2.png) [@jwmerrill](https://discourse.julialang.org/u/jwmerrill)\
**Post date:** [November 29, 2018, 7:07pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/16 "2018-11-29T19:07:23Z")

</div>

> [@mhw](#):
>
> I was also surprised that the standard matrix exponenation method (mfib) was only a “bit slower” than afib

Julia already uses the efficient power-by-squaring method for matrix exponentiation. You can use the Base implementation to compute integer powers of any type that implements \* by calling Base.power\_by\_squaring.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [November 29, 2018, 7:44pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/17 "2018-11-29T19:44:31Z")

</div>

> [@stevengj](#):
>
> The `sqrt` function for `BigFloat` is calling `mpfr_sqrt` from the MPFR library; if MPFR is using a “naive” algorithm, you should file an issue with MPFR, but I’m skeptical of that assertion without evidence.

Most of the time for the Binet formula calculation actually goes into the exponentiation:

```Julia
julia> using Nemo; ff = RealField(8388608);

julia> @time five = sqrt(ff(5));
  0.129010 seconds (1.77 k allocations: 41.310 MiB, 0.66% gc time)

julia> @time (1+five)^(10^7);
  1.668563 seconds (109 allocations: 289.760 MiB, 3.17% gc time)

```

So that is where the problem is with the Binet formula, the problem is actually not with the square root.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 29, 2018, 8:13pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/18 "2018-11-29T20:13:46Z")

</div>

> [@chakravala](#):
>
> So that is where the problem is with the Binet formula, the problem is actually not with the square root.

The `BigFloat^Integer` operation is also done by calling an optimized MPFR routine. In short, I don’t think there is any problem with the implementation being especially naive here—although you could save a few operations by being more careful, there’s no way to rewrite the Binet calculations in floating-point arithmetic to save orders of magnitude in time. The problem, as I said, is that the Binet formula requires such high precision that it is effectively useless (compared to `BigInt` methods) if you want an exact result. (It is very nice if you just want an approximate answer, however!)

Another way to see a problem with the Binet formula is to look at the complexity. Exponentiation via power-by-squaring requires O(log n) multiplications, so at first glance it may seem comparable to the O(log n) `BigInt` calculations for matrix exponentiation or for Dijkstra’s recurrence with memoization. This is deceptive, however, because not all bignum multiplications have the same cost. For the `BigInt` calculations, the numbers start small (hence cheap) and get bigger (more expensive) towards the end. For Binet’s formula via `BigFloat` arithmetic, however, you have to use the full precision from the very beginning (to get the square root accurately). If you work it out, I think this leads to roughly an extra log(n) factor in the complexity.

---

<div class="post-metadata">

**Author:** ![jwmerrill](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jwmerrill/32/6261_2.png) [@jwmerrill](https://discourse.julialang.org/u/jwmerrill)\
**Post date:** [November 29, 2018, 11:52pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/20 "2018-11-29T23:52:32Z")

</div>

Here’s one more variation that’s somewhere in between Matrix exponentiation and the memoized Dijkstra recurrence.

Using the recurrence

```julia
F(n*m) = F(m+1)*F(n+1) - F(n-1)*F(m-1)
F(n*m-1) = F(n)*F(m) + F(n-1)*F(m-1)

```

we can implement a multiplication on Fibonacci triples:

```julia
import Base.*

struct FibTrip{T}
  a::T
  b::T
  c::T
end

function *(x::FibTrip, y::FibTrip)
    q = x.c*y.c
    b = x.a*y.a - q
    c = x.b*y.b + q
    FibTrip(b+c, b, c)
end

```

Then using your `pow` function for power by squaring,

```julia
function pow(a,n) # n > 0 !
    if n == 1
        return a
    else
        c = pow(a, n >> 1)
        if iseven(n)
            return c*c
        else
            return c*c * a
        end
    end
end

```

we can get any Fibonacci number by exponentiating these triples:

```julia
function fibt(n)
    n < 0 && return iseven(n) ? -fibt(-n) : fibt(-n)
    n == 0 && return BigInt(0)
    n == 1 && return BigInt(1)
    return pow(FibTrip(BigInt(1), BigInt(1), BigInt(0)), n-1).a
end

```

```julia
> @elapsed fibt(10^9)
13.762949728

```

I’m finding that this is 30-40% slower than the memoized algorithm, and a factor of 3-5 times faster than the matrix exponentiation algorithms.

One thing that surprised me a lot is that using `Base.power_by_squaring` in place of `pow` slows this algorithm down by a factor of 2. I spent a little bit of time trying to figure out why this might be and came up with nothing. Maybe the multiplications are accumulated in a slightly different order?

One thing that I think could improve this algorithm is if Julia’s BigInts had operations that could write to a pre-allocated destination instead of always allocating new memory to store the result. Similar to `mul!` for matrices.

---

<div class="post-metadata">

**Author:** ![jakewilliami](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakewilliami/32/18945_2.png) [@jakewilliami](https://discourse.julialang.org/u/jakewilliami)\
**Post date:** [December 28, 2022, 12:12pm UTC](https://discourse.julialang.org/t/efficient-implementations-of-fibonacci-function-with-interesting-results/18123/21 "2022-12-28T12:12:12Z")

</div>

Note that this is implemented as `fibonaccinum` in Combinatorics.jl, so you don’t have to implement it yourself.