#! /usr/bin/env python3
#
def agm ( x, y ):

#*****************************************************************************80
#
## agm() estimates the arithmetic geometric mean of x and y.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    24 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  n = 0
  eps = np.finfo ( np.float ).eps
  tol = 100.0 * eps * max ( abs ( x ), abs ( y ) )

  while ( True ):
    if ( n == 0 ):
      a, g = max ( x, y ), min ( x, y )
    else:
      a, g = 0.5 * ( a + g ), np.sqrt ( a * g )
#   print ( n, a, g )
    if ( abs ( a - g ) <= tol ):
      break
    n = n + 1

  return a, n

def agm_test ( ):

#*****************************************************************************80
#
## agm_test() tests agm().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    24 June 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  print ( '' )
  print ( 'agm_test():' )
  print ( '  Test agm().' )

  x = 4.0
  y = 12.0
  z, n = agm ( x, y )
  print ( x, y, z, n )

  for i in range ( 0, 5 ):
    x = 100.0 * np.random.rand ( )
    y = 100.0 * np.random.rand ( )
    z, n = agm ( x, y )
    print ( x, y, z, n )

  print ( '' )
  print ( 'agm_test():' )
  print ( '  Normal end of execution.' )

  return

if ( __name__ == "__main__" ):
  agm_test ( )