
Fast Fibonacci numbers using Monoids - lelf
http://www.haskellforall.com/2020/04/blazing-fast-fibonacci-numbers-using.html
======
edflsafoiewq
Exponentiation-by-squaring is usually phrased as the problem of computing X^n,
but its actually usually more natural to ask about X^n Y. Unlike finding X^n
directly, this reduces directly to another instance of the same problem.

    
    
      X^n Y = (X^2)^(n/2) Y if n is even
      X^n Y = X^(n-1) (XY) if n is odd
    

So you don't really need to take a full matrix-matrix product; you only need
to square a matrix, X^2, and take a matrix-vector product, XY.

In the case of Fibonacci matrices, they can be represented with their first
column only since they have the form

    
    
      A  B
      B A+B
    

so these are just

    
    
      square (A, B) = (AA+BB, 2AB+BB)
      mulVec (A, B) (u, v) = (Au+Bv, Av+B(u+v))

------
chillee
> You can more easily generalize this solution to other arithmetic sequences.

Fun trick, if you ever find yourself in a situation like this - instead of
calculating the matrix by hand (a tedious chore), you can use an algorithm
called Berlekamp Massey on the terms of your sequence and then find the k-th
term in O(N log N log K) time instead of O(N^3 log K) time (like matrix
exponentiation would).

~~~
superpermutat0r
Here's a nice blog post about the algorithm
[https://codeforces.com/blog/entry/61306](https://codeforces.com/blog/entry/61306)
.

It definitely makes sense, recurrence relation does have its own polynomial.
Even the naive polynomial multiplication is very fast.

Now, what I'm wondering is if this is possible to do with adjacency matrices
of graphs.

Recently I was calculating the number of paths of a knight going from one
position to another in K moves. Ended up doing O(N^6 log K) (N is the
chessboard dimension) which was fast enough for a small chessboard but maybe
there's a way to turn that adjacency matrix to a polynomial and do the same.

~~~
chillee
Yes you can.

Cayley Hamilton theorem gives a linear relation between A^N and the terms less
than that.

So for your problem, you can calculate the number of paths from A->B of length
k in k*N^2 time, you need O(N^3) to compute enough terms for Berlekamp Massey,
O(N^2) to compute Berlekamp Massey, and then O(n log n log k) to compute the
k-th linear recurrence term with FFT.

I've actually considered writing a blog post about this :) It's a cool example
of how much further you can go with a problem where many people only know the
O(NK)DP solution.

~~~
mrgriffin
Please write a blog post!

------
danbruc
_… in fact, you can easily compute up to f(10^8) in a couple of seconds using
this code (not shown, because the result takes far longer to print than to
compute)._

Is this actually this fast or did the author get fooled by lazy evaluation or
something along that line? For the fun of it I tried to come up with an
algorithm in C# that works more or less like the matrix exponentiation one but
without explicitly using matrix exponentiation. When finished it failed to
terminate for large numbers and I eventually estimated that it would take
about six hours to find the 100 millionth Fibonacci number with this code on
my laptop. I don't see why this would be so much slower than a straight
forward implementation of the matrix exponentiation algorithm, where is the
big difference if any? Or is the .NET big integer implementation that slow? Or
did the author get fooled after all?

    
    
      private static BigInteger Fibonacci(UInt32 n)
      {
        if (n == 0)
        {
          return BigInteger.Zero;
        }
        else
        {
          var (a, b) = (BigInteger.Zero, BigInteger.One);
    
          var mask = 1u << (Int32)Math.Floor(Math.Log(n, 2));
    
          while ((mask >>= 1) != 0)
          {
              var s = a + b;
    
              var (a2, b2, s2) = (a * a, b * b, s * s);
    
              (a, b) = ((n & mask) == 0) ? (a2 + b2, s2 - a2) : (s2 - a2, s2 + b2);
          }
    
          return b;
        }
      }

~~~
edflsafoiewq
I transcribed your code into Python and it took about 90s to do
fib(10[star][star]8) on my T420. Language doesn't matter of course, this is
just an exercise of the bignum library.

Nice idea to use s.

~~~
danbruc
Okay, so the .NET big integer implementation is like two orders of magnitude
slower than it could be? I will have a look into that, I always assumed that
all big integer libraries are more or less comparable in performance. Thanks
for the hint.

EDIT: I found an issue [1] concerning .NET big integer performance, it seems
it is only using pretty basic and relatively slow algorithms, especially for
multiplication which is the relevant operation here.

[1]
[https://github.com/dotnet/runtime/issues/14407](https://github.com/dotnet/runtime/issues/14407)

------
srean
The blog post author abandons the closed form solution of (φ^n - ψ^n) / (φ -
ψ) citing loss of accuracy to demonstrate the pretty neat matrix solution.

The closed form solution is not quite a dead end. The floating point errors
can be avoided if one works in the field of number of the form a + b √5 where
a and b are rational.

All one needs to is to maintain the numbers a and b separately. The same
doubling trick applies for exponentiation.

~~~
joppy
This is equivalent to the matrix solution. The matrix has eigenvalues φ and ψ,
so working with linear combinations of powers of the matrix is the same as
working over Q(φ, ψ) = Q(√5). This is not just an equivalence abstractly
mathematically, but computationally the solutions are identical and run within
a constant factor of each other.

~~~
nimish
I believe it's equivalent to working in a basis that diagonalizes the matrix
rather than working in the standard one. Much easier to work with.

~~~
contravariant
Indeed (φ^n - ψ^n) / (φ - ψ) is what you'd get if you diagonalize the matrix.

However the _field_ you'd be working over is Q[√5], although really we're
interested in the subring Z[φ] here. And since φ is (trivially) an eigenvalue
of the matrix, let's call it M, this is equivalent to calculating in Z[M] by
the Cayley-Hamilton theorem.

Of course some of these constructions might be easier to work with, and
certainly something like Z[x]/(x^2-x-1) allows for fairly straightforward
calculation, but you're still doing the same multiplications essentially.

And if you're working in a convenient ring you don't really need the ψ^n term,
since you can just use that φ^n = F_{n-1} + F_{n} φ.

------
olooney
Very neat. Picking up the exponentiation by squaring algorithm for free by
implementing Semigroup is a really nice touch.

Since you're defining your own custom Matrix2x2 class and <> operation anyway,
you can save a little space and computation by relying on the fact that every
Fibonacci matrix is of the form: [ a+b, a; a, b ]. Thus only two member
variables are needed, and <> can be implemented with fewer multiplications.

------
phonebucket
The essence of using Exponentiation by Squaring to calculate Fibonacci
sequences in log(N) time is also covered in SICP exercise 1.19.

The elegance of the Haskell solution is that you only need to define a Monoid,
and then you get exponentiation by squaring for free. This is a great example
of abstraction done right.

------
s5ma6n
I see monoids more and more everyday. A bit outside the topic but does anyone
have a good introduction-level resource to learning monoids?

~~~
colonwqbang
It depends on what you mean by introduction. I think the Typeclassopedia entry
on Semigroups and Monoids is pretty good, if you are familiar with basic
Haskell syntax.

[https://wiki.haskell.org/Typeclassopedia#Semigroup](https://wiki.haskell.org/Typeclassopedia#Semigroup)

~~~
s5ma6n
Thanks for the link.

I am not good with Haskell but I meant in general mathematical / algebraic
introduction, explanation or a lecture.

~~~
JCoder58
Try this online book by Bartosz Milewski:

[https://bartoszmilewski.com/2014/10/28/category-theory-
for-p...](https://bartoszmilewski.com/2014/10/28/category-theory-for-
programmers-the-preface/)

~~~
colonwqbang
He also has lectures on YouTube which are great.

------
spekcular
> You can more easily generalize this solution to other arithmetic sequences.

The closed form as sums of exponents also easily generalizes. You can find the
roots of an associated polynomial, and the those roots are the numbers you
take powers of. This is explained on Wikipedia here:
[https://en.wikipedia.org/wiki/Recurrence_relation#Solving_ho...](https://en.wikipedia.org/wiki/Recurrence_relation#Solving_homogeneous_linear_recurrence_relations_with_constant_coefficients).

~~~
phonebucket
As an aside, a lovely (freely available) book on using generating functions to
calculate closed form solutions to such problems by the late Herbert Wilf:
[https://www.math.upenn.edu/~wilf/DownldGF.html](https://www.math.upenn.edu/~wilf/DownldGF.html)

One of the more enjoyable maths texts I've read.

~~~
JadeNB
Nothing to add to this except to say that Generatingfunctionology is an
awesome book. Although I haven't read it, another in the same space is A=B,
also co-authored by Wilf
([https://www.math.upenn.edu/~wilf/AeqB.html](https://www.math.upenn.edu/~wilf/AeqB.html))
(another co-author is Zeilberger, whose name anyone interested in
computational discrete math will recognise).

~~~
phonebucket
I've spent some time on A = B as well. Very pleased to find someone else
interested in these things on here!

------
dwilding
Several years ago I hacked together a little tool for exploring integer
sequences that can be generated via matrix-vector products:

[http://dpw.me/clear-the-deck/](http://dpw.me/clear-the-deck/)

Sadly it's not as functional as it once was (sequences are supposed to
continue generating as you scroll down the page). If anyone's interested in
lending a hand to get it working properly again, I'd love to hear from you! My
email is in my profile.

~~~
dwilding
I wrote a blog post about it here at the time:

[https://aperiodical.com/2014/06/discovering-integer-
sequence...](https://aperiodical.com/2014/06/discovering-integer-sequences-by-
dealing-cards/)

