function c = FastExp(a,b,n)
% c = FastExp(a,b,n)
% Input:  a, b, < n positive integers, n = the modulus
% Output: c = a^b (mod n) 
% Program uses the fast exponentiation method to calculate c = a^b (mod n).  
% Beginning with factor = a^1: 
% on each pass, it is determined whether this is a desired factor.  If so, then 
% it is multiplied into a variable containing a cumulative product of desired 
% factors.  Then the factor is squared, and the process repeated.

expon = b;
c = ones(size(a));
factor = a;

while(expon > 0)
    if (mod(expon,2) == 1)
        c = mod(c .* factor,n);
    end
    
    expon = floor(expon/2);
    factor = mod(factor.*factor,n);
end