## Monday, 3 December 2012

My previous post raised some interest (i.e., one person asked me about it), so I decided to follow up a little bit on that. In a nutshell, I suggested to use logarithms and exponents to raise large matrices to some non-integer power because it's faster. That is, `exp(log(A)*x)` is faster than `A.^x`, subject to the constraint that `A>0` everywhere. Referring to this, I can confirm that it works on Octave as well as in Matlab. I wouldn't know about other programming languages though, if you feel like trying and letting me know in the comments, I'll create a new post about it.

To that, I wanted to add something else. If we want to raise a large matrix to some integer power, it's worth using the exponentiation by squaring. The following Matlab code is significantly faster than the `.^` operator already when the exponent is 4.

```function y = pow(x,p)

y = ones(size(x));

while (p > 0)
if (mod(p,2) == 1)
y = y.*x;
end

p = floor(p/2);
x = x.*x;
end
```

Interestingly, this does not work for Octave, which probably implements a similar algorithm when it realises that the exponent is an integer. So, if you use Octave keep going with your `.^` operator.

In Matlab the code doesn't look great, probably because most of its commands are native for double precision floating-point arithmetic, while here we're dividing the exponent which is an integer. Somehow, the same function written in C looks a little bit more neat.

```float powN(float x, unsigned int p) {
float y = 1.0;

while (p > 0) {
if (p & 1)
y *= x;

x *= x;
p >>= 1;
}

return y;
}
```

Here we can use directly bitwise operators in the condition and in place of the division by 2. That said, this is a very low-level operation which is likely to be implemented more efficiently in whichever API you're using (except for Matlab, it seems). So I'd recommend you to double-check first, because it might very well be that you can use a vectorised version that outperforms this one.