Power algorithm for finding rem(a^n,m)

% This implements the power algorithm for finding rem(a^n,m) (the remainder on dividing a^n by m) reliably for large numbers (up to about 8 digits).
%
% Usage: Type pow(a,n,m) where a,n,m are already defined, or are explicit numbers as in pow(2,1000,25).


function x=pow(a,n,m)
b = a;
x = 1;
while n > 0
d=rem(n,2);
if d==1
x=rem(x*b,m);
end
b=rem(b * b,m);
n=(n-d)/2;
end