3 ms·

The comparison with the closed form is weird to me: since it uses sqrt(5), I suppose the author intended to use floating point numbers to compute it. But when d

The comparison with the closed form is weird to me: since it uses sqrt(5), I suppose the author intended to use floating point numbers to compute it. But when doing so: - You get a constant time algorithm, since the power function is constant time on most mathematical librarie. Without entering into details on how the pow function is actually computed, on real number you can compute it with: pow(x,y) = exp(y * log(x)). It is not visible on the performance graph probably because the matrix-based algo and the closed formula are dwarfed by the linear time algo. - You gen an incorrect result very quickly: not even considering the fact that sqrt(5) can not be represented exactly by a floating point number, the result will be out of the range of floating point numbers with n proportional to the number of bits your floating point number holds, so probably between n=50 and n=100.

abhi9u 3y ago
Author here: At the time of writing, I was not aware that the matrix form and the closed form are same, because you can derive the golden ratio based formula from the matrix formula. Full disclaimer: I am not a math major, and I can make such mistakes :) Although, I am curious to learn about exponentiation of real numbers. While writing this article, I was investigating how different libraries/languages implement exponentiation. While for matrix exponentiation, numpy used the algorithm I showed in the article. I also looked at the pow() function in libm of FreeBSD and it used the form you provided. I've been trying to look at an explanation behind this. Please share a pointer if you have any references.
- SuchAnonMuchWow 3y ago
  The main issue with floating point number is that it doesn't act like a ring (https://en.wikipedia.org/wiki/Ring_(mathematics) https://en.wikipedia.org/wiki/Ring_(mathematics)) because of rounding errors. For example the equality x*(y*z) = (x*y)*z doesn't hold because each multiplication introduce small rounding errors. And this kind of equations are required for the exponentiation-by-squaring method (consider for example how x**5 is computed: x**5 = (((x*x)*x)*x)*x) using a naive exponentiation approach, but also x**5 = x*(x*x)*(x*x) = x*(x*x)**2 using the exponentiation-by-squaring method). Integer operations and matrix operations are both rings, so the exponentiation-by-squaring approach works. For floating points, it doesn't give an exact result, but could still be used in practice to get an approximate result. The pow function doesn't use this approach however, because both arguments of the pow function are floating points numbers, and can thus compute things like pow(5, 0.5) = sqrt(5), or even 534.2**12.7. The exponentiation-by-squaring approach works only when the exponent is an integer, so we need other solutions. How exactly the pow function (or other math function like exp, log, cos, ...) are computed is a vast research subject in itself, but if you are interested, the Handbook of Floating-Point Arithmetic (https://perso.ens-lyon.fr/jean-michel.muller/Handbook.html https://perso.ens-lyon.fr/jean-michel.muller/Handbook.html) is both a good introduction to the domain (and the problems associated with the finite precision of floating point number), and provides a lot of techniques to work around or minimize the rounding errors that happens when manipulating floating point numbers.
  - abhi9u 3y ago
    Thank you so much, both for providing more background on exponentiation, and for the reference. :-)
contravariant 3y ago
You can also use exact arithmetic. This is not as bad as it seems, all values you'll encounter are of the form a + sqrt(5) b, with a and b rational (and not even all rationals, but I can't be bothered). You can even make this somewhat more concrete by using the matrix representation* of a + sqrt(5) b: [[a, 5b], [b, a]] Of course this reduces the whole problem to a matrix exponential again, so it's somewhat pointless. *: The matrix representation is what you get if you view each value a + sqrt(5)b as a vector, and notice that multiplication by a constant is a linear map. The matrices are these linear maps, and it can be shown that matrix addition and multiplication preserve the ring structure. You can forget about the vectors at this point.
- SuchAnonMuchWow 3y ago
  That's true. Integers of this form a+b*sqrt(5) are called quadratic integers (https://en.wikipedia.org/wiki/Quadratic_integer https://en.wikipedia.org/wiki/Quadratic_integer), and define a ring for which you can define addition and multiplication rules: (a+b*sqrt(5)) + (c+d*sqrt(5)) = (a+c) + (b+d)*sqrt(5) (a+b*sqrt(5)) * (c+d*sqrt(5)) = (ac+5bd) + (ad+bc)*sqrt(5) Note that seeing the operations this way reduce a bit the number of operations compared to the matrix representation you provided: a multiplication using the formula above require 5 integer multiplications instead of 8 with the matrix [[a, 5b], [b, a]]. But this is only a constant factor improvement: the algorithm complexity is still log(n). And the same result can be achieved by using the special structure of the matrix: for your representation, the matrix are all toeplitz matrices, and multiplying those kind of matrices require less operations than full generic matrix multiplication. For the representation in the original article, If I remember correctly, the matrices are triangular, and again you can multiply them more efficiently using this specificity.

abhi9u 3y ago

Author here: At the time of writing, I was not aware that the matrix form and the closed form are same, because you can derive the golden ratio based formula from the matrix formula. Full disclaimer: I am not a math major, and I can make such mistakes :) Although, I am curious to learn about exponentiation of real numbers. While writing this article, I was investigating how different libraries/languages implement exponentiation. While for matrix exponentiation, numpy used the algorithm I showed in the article. I also looked at the pow() function in libm of FreeBSD and it used the form you provided. I've been trying to look at an explanation behind this. Please share a pointer if you have any references.

SuchAnonMuchWow 3y ago
The main issue with floating point number is that it doesn't act like a ring (https://en.wikipedia.org/wiki/Ring_(mathematics) https://en.wikipedia.org/wiki/Ring_(mathematics)) because of rounding errors. For example the equality x*(y*z) = (x*y)*z doesn't hold because each multiplication introduce small rounding errors. And this kind of equations are required for the exponentiation-by-squaring method (consider for example how x**5 is computed: x**5 = (((x*x)*x)*x)*x) using a naive exponentiation approach, but also x**5 = x*(x*x)*(x*x) = x*(x*x)**2 using the exponentiation-by-squaring method). Integer operations and matrix operations are both rings, so the exponentiation-by-squaring approach works. For floating points, it doesn't give an exact result, but could still be used in practice to get an approximate result. The pow function doesn't use this approach however, because both arguments of the pow function are floating points numbers, and can thus compute things like pow(5, 0.5) = sqrt(5), or even 534.2**12.7. The exponentiation-by-squaring approach works only when the exponent is an integer, so we need other solutions. How exactly the pow function (or other math function like exp, log, cos, ...) are computed is a vast research subject in itself, but if you are interested, the Handbook of Floating-Point Arithmetic (https://perso.ens-lyon.fr/jean-michel.muller/Handbook.html https://perso.ens-lyon.fr/jean-michel.muller/Handbook.html) is both a good introduction to the domain (and the problems associated with the finite precision of floating point number), and provides a lot of techniques to work around or minimize the rounding errors that happens when manipulating floating point numbers.
- abhi9u 3y ago
  Thank you so much, both for providing more background on exponentiation, and for the reference. :-)

contravariant 3y ago

You can also use exact arithmetic. This is not as bad as it seems, all values you'll encounter are of the form a + sqrt(5) b, with a and b rational (and not even all rationals, but I can't be bothered). You can even make this somewhat more concrete by using the matrix representation* of a + sqrt(5) b: [[a, 5b], [b, a]] Of course this reduces the whole problem to a matrix exponential again, so it's somewhat pointless. *: The matrix representation is what you get if you view each value a + sqrt(5)b as a vector, and notice that multiplication by a constant is a linear map. The matrices are these linear maps, and it can be shown that matrix addition and multiplication preserve the ring structure. You can forget about the vectors at this point.

SuchAnonMuchWow 3y ago
That's true. Integers of this form a+b*sqrt(5) are called quadratic integers (https://en.wikipedia.org/wiki/Quadratic_integer https://en.wikipedia.org/wiki/Quadratic_integer), and define a ring for which you can define addition and multiplication rules: (a+b*sqrt(5)) + (c+d*sqrt(5)) = (a+c) + (b+d)*sqrt(5) (a+b*sqrt(5)) * (c+d*sqrt(5)) = (ac+5bd) + (ad+bc)*sqrt(5) Note that seeing the operations this way reduce a bit the number of operations compared to the matrix representation you provided: a multiplication using the formula above require 5 integer multiplications instead of 8 with the matrix [[a, 5b], [b, a]]. But this is only a constant factor improvement: the algorithm complexity is still log(n). And the same result can be achieved by using the special structure of the matrix: for your representation, the matrix are all toeplitz matrices, and multiplying those kind of matrices require less operations than full generic matrix multiplication. For the representation in the original article, If I remember correctly, the matrices are triangular, and again you can multiply them more efficiently using this specificity.