#! /usr/bin/env python
#
def 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:
#
#    05 April 2013
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer A, the base of the expression to be tested.
#    A should be nonnegative.
#
#    Input, integer N, the power to which the base is raised.
#    N should be nonnegative.
#
#    Input, 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 x

  if ( int ( a ) != a ):
    x = -1
    return x

  if ( n < 0 ):
    x = -1
    return x

  if ( int ( n ) != n ):
    x = -1
    return x

  if ( m <= 0 ):
    x = -1
    return x

  if ( int ( m ) != m ):
    x = -1
    return x
#
#  A contains the successive squares of A.
#
  x = 1

  while ( 0 < n ):

    d = ( n % 2 )

    if ( d == 1 ):
      x = ( ( x * a ) % m )

    a = ( ( a * a ) %  m )
    n = ( ( n - d ) // 2 )
#
#  Ensure that 0 <= X.
#
  while ( x < 0 ):
    x = x + m

  return x

def power_mod_test ( ):

#*****************************************************************************80
#
## POWER_MOD_TEST tests POWER_MOD.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2014
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'POWER_MOD_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  POWER_MOD computes the remainder of a power' )
  print ( '  of an integer modulo another integer.' )

  a = 7
  n = 50
  m = 11

  x = power_mod ( a, n, m )

  print ( '' )
  print ( '  A = %d' % ( a ) )
  print ( '  N = %d' % ( n ) )
  print ( '  M = %d' % ( m ) )
  print ( '  mod ( A^N, M ) = %d' % ( x ) )

  a = 3
  n = 118
  m = 119

  x = power_mod ( a, n, m )

  print ( '' )
  print ( '  A = %d' % ( a ) )
  print ( '  N = %d' % ( n ) )
  print ( '  M = %d' % ( m ) )
  print ( '  mod ( A^N, M ) = %d' % ( x ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'POWER_MOD_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  power_mod_test ( )
  timestamp ( )