%binary exponentiation modulo N;
function[Rs]=emn(x,r,N); %function calculates x^r modN;
rb=bitget(r,[1:ceil(log2(r+1))]); %binary representing r
l=length(rb); y=x; Rs=1; %initialization of raising;
for k=1:l
    if rb(k)==1 Rs=rem(Rs*y,N); end; %calculating current product of degrees of x modN;    
    y=rem(y^2,N); %passing on to the next step;
end;