number_powermod — Modular exponentiation.
x = number_powermod ( a, n, m )
a 1x1 matrix of floating point integers, must be positive
a 1x1 matrix of floating point integers, must be positive
a 1x1 matrix of floating point integers, must be positive
a 1x1 matrix of floating point integers
Returns x = mod ( a^n, m ).
Uses a repeated squaring algorithm.
The algorithm may fail if the number is so large that the
square a^2
is out of range.
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 then A**13 = A * A**4 * A**8. This leads to 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.