function x = power_mod ( a, n, m )
%*****************************************************************************80
%
%% POWER_MOD computes mod ( A^N, M ).
%
% Discussion:
%
% Some programming tricks are used to speed up the computation, and to
% allow computations in which A^N is much too large to store in a
% real word.
%
% First, for efficiency, the power A^N is computed by determining
% the binary expansion of N, then computing A, A^2, A^4, and so on
% by repeated squaring, and multiplying only those factors that
% contribute to A^N.
%
% Secondly, the intermediate products are immediately "mod'ed", which
% keeps them small.
%
% For instance, to compute mod ( A^13, 11 ), we essentially compute
%
% 13 = 1 + 4 + 8
%
% A^13 = A * A^4 * A^8
%
% mod ( A^13, 11 ) = mod ( A, 11 ) * mod ( A^4, 11 ) * mod ( A^8, 11 ).
%
% Fermat's little theorem says that if P is prime, and A is not divisible
% by P, then ( A^(P-1) - 1 ) is divisible by P.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 15 November 2004
%
% Author:
%
% John Burkardt
%
% Input:
%
% integer A, the base of the expression to be tested.
% A should be nonnegative.
%
% integer N, the power to which the base is raised.
% N should be nonnegative.
%
% integer M, the divisor against which the expression is tested.
% M should be positive.
%
% Output:
%
% integer X, the remainder when A^N is divided by M.
%
if ( a < 0 )
x = -1;
return
end
if ( floor ( a ) ~= a )
x = -1;
return
end
if ( n < 0 )
x = -1;
return
end
if ( floor ( n ) ~= n )
x = -1;
return
end
if ( m <= 0 )
x = -1;
return
end
if ( floor ( m ) ~= m )
x = -1;
return
end
%
% A contains the successive squares of A.
%
x = 1;
while ( 0 < n )
d = mod ( n, 2 );
if ( d == 1 )
x = mod ( x * a, m );
end
a = mod ( a * a, m );
n = floor ( ( n - d ) / 2 );
end
%
% Ensure that 0 <= X.
%
while ( x < 0 )
x = x + m;
end
return
end