function rem = powMod(base,exponent,modulo)
% POWMOD(BASE,EXPONENT,MODULO) computes BASE^EXPONENT mod MODULO using the
% right-to-left binary method.
if (modulo^2 > bitmax)
error('Modulos is too large!')
end
rem = 1;
base = mod(base,modulo);
e = fliplr(str2num(dec2bin(exponent)')'); % this is the best/worst line ever; lsb in e(1).
n = length(e);
for i = 1:n
if e(i) == 1
rem = mod(rem*base,modulo);
end
base = mod(base^2,modulo);
end
Notice that even here we could have double precision overflow if the square of the modulus is too large. I've been thinking about either implementing a variable precision method, or switching out of Matlab entirely the next time I run the course. I think a natural alternative would be Python. If anyone has any thoughts, feel free to post or email.