r8_agm <- function ( a, b )

#*****************************************************************************80
#
## r8_agm() computes the arithmetic-geometric mean of A and B.
#
#  Discussion:
#
#    The AGM is defined for nonnegative A and B.
#
#    The AGM of numbers A and B is defined by setting
#
#      A(0) = A,
#      B(0) = B
#
#      A(N+1) = ( A(N) + B(N) ) / 2
#      B(N+1) = sqrt ( A(N) * B(N) )
#
#    The two sequences both converge to AGM(A,B).
#
#    In Mathematica, the AGM can be evaluated by
#
#      ArithmeticGeometricMean [ a, b ]
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    07 September 2026
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Stephen Wolfram,
#    The Mathematica Book,
#    Fourth Edition,
#    Wolfram Media / Cambridge University Press, 1999.
#
#  Input:
#
#    real A, B, the arguments whose AGM is to be computed.
#    0 <= A, 0 <= B.
#
#  Output:
#
#    real VALUE, the arithmetic-geometric mean of A and B.
#
{
  it_max <- 1000

  if ( a < 0.0 )
  {
    cat ( '\n' )
    cat ( 'r8_agm(): Fatal error!\n' )
    cat ( '  Argument A < 0.\n' )
    stop ( 'r8_agm(): Fatal error!' )
  }

  if ( b < 0.0 )
  {
    cat ( '\n' )
    cat ( 'r8_agm(): Fatal error!\n' )
    cat ( '  Argument B < 0.\n' )
    stop ( 'r8_agm(): Fatal error!' )
  }

  if ( a == 0.0 | b == 0.0 )
  {
    value <- 0.0
    return ( value )
  }

  if ( a == b )
  {
    value <- a
    return ( value )
  }

  it <- 0
  tol <- 100.0 * .Machine$double.eps
  a1 <- a
  b1 <- b

  while ( TRUE )
  {
    it <- it + 1

    a2 <- ( a1 + b1 ) / 2.0
    b2 <- sqrt ( a1 * b1 )

    if ( abs ( a2 - b2 ) <= tol * ( a2 + b2 ) )
    {
      break
    }

    if ( it_max < it )
    {
      cat ( '\n' )
      cat ( 'r8_agm(): Warning!\n' )
      cat ( '  No convergence.\n' )
      break
    }

    a1 <- a2
    b1 <- b2
  }

  value <- a2

  return ( value )
}