function b = powermod(a,z,n)
% This function calculates y = a^z mod n
% If z is a matrix, it calculates a^z(j,k) mod for every element in a
[ax,ay]=size(z);


% If a is negative, put it back to between 0 and n-1
g=mod(a,n);
binaryExpansion = dec2bin(z)
a = g;
b = ones(ax,ay);
l = length(binaryExpansion(1,:));
%COMPUTE THE POWERS OF a:
powers = [1 a*ones(1,l)];
for I = [2:l]
    powers(I+1) = mod(powers(I)^2,n);
end
powers
bitmasks = zeros(ax*ay,l);
comparison = dec2bin(ones(ax*ay,1));
for I = [1:l]
    bitmasks = binaryExpansion(:,l-I+1) == comparison;
    indices= find(bitmasks);
    b(indices) = mod(b(indices)*powers(I+1),n);
end