A linear algebra trick for computing Fibonacci numbers fast

What is this matrix doing? When you left multiply it, you're summing the elements of the left column and putting them in the top left (next fib number). Then you're moving the current number to bottom left (and top right). Incidentally, you move the previous to bottom right.

That's interesting!! It's very simple

But, all these method, in practicality, leads to integer overflows.

Works very well for smaller Fibonacci numbers, but not for, say, the 1500th Fibonacci number.

That is very easy and practical to do instead with the naive formula using an array based approach.

    fibs = [0, 1]
    for i in range(2,n):
        fibs.append(fibs[-2] + fibs[-1])

I have been meaning to write it up, but dealing with a bad flu.

The closed form leads to approximation errors.

I tried Julia, Numpy, Numba, Jax, Rust, but integer overflows are there in any parallelized approach.

The naive implementation with LRU caching is the fastest and most reliable so far.

I am surely missing something, you can actually use some things to overcome integer overflows, right?

I hoped to see those in this article, but they were simply not there.

Can anyone else help?

Even UInt128 is not enough for 1500th Fibonacci number. The upside of Rust, is that it doesn't fail silently. The compiler is really helpful.

Jax, Julia, Numpy confidently give wrong answers. Only way to get the feel is to eyeball a bunch of entries and check for negative numbers. Very easy to spot overflows when you are at, say, 1600th fibonacci, but one can feel confident on wrong answers on smaller numbers.

You can carry out the exponentiations in the field Q[sqrt(5)], which is two-dimentional over Q. The interesting thing here is that diagonalization is trading one 2d vector space for another -- one with a very well-behaved multiplication operation (it's distributive, commutative, associative, and invertible).

There's no need to do this to quickly compute fibonacci numbers, but it's pretty neat.

For all other things that I tried, I either get integer overflows (for linear algebra based implementations) or floating point errors (for the closed form).