Prev: Blockwise Matrix Expansion
Next: help with matlab code
From: Bruno Luong on 20 Jul 2010 02:49
I dig out the old thread where the binary decomposition power algorithm is provided (need to change '*' to '.*' to adapt to the context):

http://www.mathworks.com/matlabcentral/newsreader/view_thread/264887#692226

Bruno
From: Bruno Luong on 20 Jul 2010 07:10
Here is a benchmark with the binary method for integer power (Matlab 2010A 64, Vista):

function TestIntOpt(NumEl)
if nargin<1; NumEl = 1e7; end

tt=rand(1,NumEl);
k = 10;

tic
x = tt.^k; %takes 1159msec
T1=toc;fprintf('Tenth power: %.0f ms\n',1000*T1);

tic
x = exp(k*log(tt)); %takes 560 msec
T1=toc;fprintf('Tenth power via exp 10 log: %.0f ms\n',1000*T1);

tic
x = intpow(tt,k); %takes 265 msec
T1=toc;fprintf('Tenth power via intpow: %.0f ms\n',1000*T1);

end % TestIntOpt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = intpow(a, n)
% d = intpow(a, n)
% power with integer exponent
% return the array d = a.^n

% binary decomposition
nb = dec2bin(n)-'0';
ak = a;
if nb(end)
d = a;
else
d = 1;
end
for nbk = nb(end-1:-1:1)
ak = ak.*ak;
if nbk
d = d.*ak;
end
end

end % intpow

% Bruno
From: dpb on 20 Jul 2010 08:38
Steve Amphlett wrote:
....

> Related...
>
> If you are writing MEX files (or just normal programs/functions) that
> use powers, consider switching to C++ and use std::pow() in place of C's
> very restrictive pow() and powf() functions. It has overloads for
> integral exponents, which can choose the optimum method (multiple
> multiplications vs logarithms vs ???).
>
> I cannot comment on what FORTRAN and its modern incarnations do for
> intergral exponents.
....

All combinations of intrinsic numeric types are permitted with the **
operator but you cannot raise a negative real value to a real power.

However, the Standard does not prescribe the implementation so the
actual algorithm/optimization is "processor dependent" (where
"processor" in the Standard means in generic terms the compiler).

For example, some compilers may recognize the case of (-1.0)**2.0 and
replace the real exponent w/ the integer; others may not. Similarly for
the subject here of what may be done in speed optimization will be a
"quality of implementation" issue not mandated by the Standard.

--
From: Bruno Luong on 20 Jul 2010 10:35
This is a more compact version of the function allpow() in post #4

function powers = allpow(x, pmax)
% function powers = allpow(x, pmax)
%
% compute all powers of exponent P = (0,1,...,PMAX) of vector X.
% return the result in 2d array POWERS having size [length(X) x PMAX+1)]

powers= zeros(1, pmax+1);

if pmax>0
powers(:,2) = x;
p2 = 1;
r = 1;
for k=2:pmax
powers(:,k+1) = powers(:,p2+1).*powers(:,r+1);
if (p2==r)
p2 = 2*p2;
r = 1;
else
r = r+1;
end
end % for-loop
end

end %% allpow
From: Irl on 20 Jul 2010 11:13
Dear Bruno, Steve, Walter, et al. --
Thanks for the interesting and useful info. My main question was whether there was any way to control how the internal vectorized power function deals with small integer powers, and so far it appears that there is not. For my application, the powers are no higher than ten or so, and my little test program suggests that for that regime the repeated-multiply method, possibly augmented with a binary decomposition, will be the fastest. I can do the binary decomposition manually for this case, although the possibility of doing it automatically is a useful trick as well.
Again, thanks for the help.
Yours,
Irl
"Irl " <irl_smith(a)raytheonRemoveTheseWords.com> wrote in message <i2236g$oco$1(a)fred.mathworks.com>...
> MATLAB seems to implement x^2 as a multiply but x^3, or higher power, via something much slower (exp of number times log?). For large enough n, exp n log will be faster than multiplying n times, but the internal code generator seems to err wildly on the side of using the slow method, whatever it is, rather than just repeated multiplies. But it does do x^2, apparently, as x*x. Is there a way to make it use repeated multiplies for higher exponents? I wrote a Zernike polynomial program and it's very tiresome to go change all the powers to multiply strings.
> Example program, including comments giving the execution times:
> function TestIntOpt(NumEl)
> if nargin<1; NumEl = 1e7; end
> % Commented results are for 1e7 as input argument
> tt=rand(1,NumEl);
> tic
> x = tt; %zero ms
> T1=toc;fprintf('Simple assign, %.0f elements: %.0f ms\n',NumEl,1000*T1);
> tic
> x = tt.^2; %takes 70 msec
> T1=toc;fprintf('Second power (square): %.0f ms\n',1000*T1);
> tic
> x = tt.^3; %takes 1200 msec
> T1=toc;fprintf('Third power (cube): %.0f ms\n',1000*T1);
> tic
> x = tt.^2.*tt; %takes 100 msec
> T1=toc;fprintf('Cube via square times first power: %.0f ms\n',1000*T1);
> tic
> x = tt.^10; %takes 1200 msec
> T1=toc;fprintf('Tenth power: %.0f ms\n',1000*T1);
> tic
> x = tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt.*tt; %takes 170 msec
> T1=toc;fprintf('Tenth power via nine multiplies: %.0f ms\n',1000*T1);
> tic
> x = exp(10*log(tt)); %takes 660 msec
> T1=toc;fprintf('Tenth power via exp 10 log: %.0f ms\n',1000*T1);
First  |  Prev  | 
Pages: 1 2
Prev: Blockwise Matrix Expansion
Next: help with matlab code