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))
> 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).
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.
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.
On an unbounded chessboard it's not too difficult to represent the movement of a knight by a polynomial (or really a 2D convolution), but when your chessboard is bounded then the boundary effects become annoying. Maybe there is some way to mitigate those but I'm not aware of any.
… 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;
}
}
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.
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.
To do f(10^8) the OP says his/her code would take the matrix and square it about 26 times (2^26 = 67 million and change), then do some more in other ways. That's where the speedup comes from. Are you also doing that?
(The mechanism for the speedup is described in passing, in a little note about the Haskell library function he's using to calculate a power of the matrix. For whatever it's worth this trick is excellent, but it dates from before computers were invented.)
Yes, as you can see in the code, the loop executes once for each bit to the right of the leftmost set bit or equivalently ⌊log₂(n)⌋ times. Naively this suggest a runtime of O(log(n)) but the numbers become huge pretty quickly and then all the arithmetic operations are far from O(1). The actual runtime is probably more like O(n²) or O(n² log(n)) after one accounts for the runtime of the big integer operations.
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.
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.
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} φ.
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.
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.
There is not much to learn, a monoid is a very simple abstraction.
Think for example about learning about the "iterator" abstraction. It's just something with "hasNext" and "next" operations, nothing more to understand. But when you know it you start seeing iterators everywhere, and you get frustrated when something that should be an iterator isn't so you cannot write a simple foreach loop. And you'll start implementing iterators in your own code, because you know that provides a lot of convenience for free to consumers.
Same thing for monoids. They're everywhere, they offer some convenience for free, and you'll be wondering why things that should have monoid instances don't.
A monoid is an algebraic data structure that has a binary operation (+, *, ++, etc), and an empty element (0, 1, [], etc). In terms of relationship to other algebraic data structure, it "inherits" the binary operation from its Semigroup buddy, which is a very nice data structure to represent non-empty containers (and statically enforce that property).
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.
In abstract algebra there really is a hierarchy of group-like-structures with different levels of symmetry (see the table in [1]), and a monoid is one of such structures.
I've had some luck with the first chapter of this (free) [2] book, although they explain groups and not monoids but for me it was then rather quick to gain intuition about monoids.
Side note: the study of monoids isn't a very big field in mathematics, because they're too simple to say too much about them. But if you modify them just a bit, either in the direction of groups or in the direction of categories, you get massive fields of study.
I think the existence of an identity isn't very important. You can adjoin an identity to any semigroup to turn it into a monoid.
In any case, doesn't adding more properties make a theory simpler? For example the classification of finite simple groups is much harder than the classification of finite simple commutative groups.
(It's not focused on Monoids. It also covers groups, rings, and fields from a mathematical point of view. It's written for an undergraduate math major with no prior experience in these areas.)
> 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....
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
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) (another co-author is Zeilberger, whose name anyone interested in computational discrete math will recognise).
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.
In the case of Fibonacci matrices, they can be represented with their first column only since they have the form
so these are just