Using the Quick Raise of Matrices to a Power to Write a Fast Interpreter
kukuruku.co
kukuruku.co
It should be quite clear that a more "powerful" language would not admit this last meaning/interpretation/denotation.
[0] http://blog.higher-order.com/blog/2014/12/21/maximally-power...
Rather, the OP is simply interested in carrying out arithmetic calculations of a sort which are trivially computable (after all, the inputs the OP is dealing with are the very programs which would compute the answer!), but notes that one can in a principled way transform naive specifications of these computations (including loops like "Do this a million times") into faster methods of attaining the same answer (by thinking of such loops as matrix exponentiations and then using quick, rather than naive, exponentiation).
If the power is zero, we will return the identity matrix. If the power is even (2N), we can recursively compute M^N, and then return the square of the obtained matrix. If the power is odd (2N+1), it’s enough to recursively compute M^2N value, and return the obtained matrix, that has been multiplied by M."
Tidbit: that likely is the best approach in practice, but perhaps surprisingly, it is not optimal in the number of matrix multiplications. See http://en.m.wikipedia.org/wiki/Addition-chain_exponentiation and https://oeis.org/A003313. Wikipedia has a link mentioning addition-subtraction chains, so it may even be beneficial to compute the inverse of the matrix in some cases (M^1023 likely is such a case, if inverting M isn't too hard)
(It's worth noting that while finding short addition chains isn't known to be NP-complete, a slight variant of the problem is.)
I think there are methods of finding short addition chains that reliably do better than the binary method, though; in TAOCP volume 2, "Evaluation of powers" section (which is a very good reference on addition chains in general!), Knuth talks about the power tree method, which from a brief skim seems like it might be better even when you include the time to find the chain. I think there might also be better methods based on continued fractions, going from memory.
For a fixed exponent, certainly; you could just spend additional time in the compiler finding short addition chains for all the exponents that are used. When exponents are allowed to vary, then we're in the case I described in my post above -- maybe using the power-tree method, or certain other methods, would improve on the binary method; but finding the shortest addition chain at runtime would not. (By which I mean, not with existing methods, and it seems unlikely that there's any fast way.)
caveat: will fail for large N due to limited precision of floats... works fine on paper though.
[update] caveat 2: assumes a model of computation where ϕ^N can be computed in constant time, which only works on paper.
ϕ_mult = [[1 5] [1 1]]
This is because ϕ(x+y√5)=(x+5y)+(x+y)√5, and
[[1 5] [1 1]] [x y] = [x+5y x+y]
Now, to compute ϕ^N, we just have to compute [x y] = ϕ_mult^N[1 0] (since [1 0] is the representation of 1=1+0√5). Then ϕ^N=x+y√5. This is all exact, no floating point arithmetic!
It's actually just an example of an algebraic number field [1] -- the idea is that "algebraic numbers" (i.e., roots of polynomials) can be equivalently thought of as vectors and matrices. The simplest example, which everyone is familiar with, is the complex numbers, where i=√(-1). In this case, [x y] represents x+y*i, and multiplication by i corresponds to the matrix i_mult=[[0 1] [-1 0]]. Then i_mult^2= -[[1 0] [0 1]], in other words, multiplication by -1!
Original article: http://habrahabr.ru/post/148901/ Post date: 2012
More from the same author: http://habrahabr.ru/users/skidanovalex/topics/
On its own, this would be great, but the issue is that these translations are not necessarily authorized in advance, and the attribution was not originally present. But there is now an 'original source' link at the top of each translation, and the translators feel they are adhering to the use policies of Habrhabr.
I was dubious initially, but at this point I'm mostly convinced that their intents are good, and I appreciate having access to these articles that I otherwise would not be able to enjoy. I think they are open to suggestions of how to improve their approach, and hope they eventually find a solution that satisfies all parties.
I'm not sure what classical algorithm of computing Fibonacci numbers OP is mentioning, but using Mathematica on my machine takes about 4 ms to compute the 1,000,001st Fibonacci number exactly…
Maybe by classical he means naive, in that case: doesn't it ever get tired to pick as an example an algorithm that will be optimized away by any half-decent compiler?
> I have recently read an article about calculating the Nth Fibonacci number in O(log N) arithmetic operations.
Either you assume fixed-precision integers, in which case you just use a lookup table and call it O(1), or you assume arbitrary-precision integers, in which case you cannot have anything less than a O(n) algorithm, as the number of digits of the Nth Fibonacci number is O(n) (as the Fibonacci sequence is exponential).
I see this all over the place, and it's frustrating every time.
A similar situation arises in, say, mergesort when we say that the algorithm takes O(n log n) comparisons. That's not the same as an O(n log n) time algorithm; for instance if you think of sorting distinct integers, each integer will take at least O(log n) bits to represent, and to compare two integers you could potentially need to read all the bits.
And as the Fibonacci sequence is exponential, any lookup table will be relatively small. (The maximum number of elements in your lookup table is linear w.r.t. the number of digits of your (fixed-precision) integers. I mean, even a 128-bit return value only requires 185 elements in the table, for a total of 2960 bytes.)
> A similar situation arises in, say, mergesort
This is why you have to specify exactly what you mean by the big-O notation you specify. Do you mean comparisons? Time? Time assuming everything but the variables specified are constant? Etc.
>It's true that you can't have an O(log N) time algorithm, as arithmetic operations cannot be performed in constant time.
Yep. You can perform some operations in sublinear time w.r.t. the number of bits, however, by exploiting hardware-level parallelization. It would be interesting to see what the best (in terms of big-O) possible implementation of fib(n) would be, given constant-time operations for a fixed-length machine word.
He did specify it, it's O(log N) arithmetic operations.
Note: I don't mean "assuming arithmetic operations take constant time, we can compute F(n) in O(log N) time", in which case you can correctly say that, since arithmetic operations don't take constant time in arbitrary precision arithmetic, we must be using fixed precision (or something essentially equivalent), in which case a lookup table takes O(1) time and O(precision) space. I mean "F(n) can be computed with O(log N) arithmetic operations". If I perform my arithmetic operations over arbitrary precision integers, and I implement simple arithmetic operations (so each arithmetic operation takes O((log n)^2) time), then F(n) takes O((log n)^3) time.
If so that's fine but it's probably a rare use case, since loops usually have conditional branches, don't they? If they don't the algorithm is poorly written in the first place as it should have been noticed that there was a linear transformation involved.
Leslie G. Valiant: General Context-Free Recognition in Less than Cubic Time. J. Comput. Syst. Sci. 10(2): 308-315 (1975)
a+=1; b+=1/a;