#! /usr/bin/env python
#
def monte_carlo_log_test ( ):

#*****************************************************************************80
#
## MONTE_CARLO_LOG_TEST tests the Monte Carlo rule on the log singular integrand.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 November 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from test_functions import integral_log
  from test_functions import integrand_log
  from r8vec_uniform_01 import r8vec_uniform_01

  print ( '' )
  print ( 'MONTE_CARLO_LOG_TEST' )
  print ( '  Test the Monte Carlo rule on the log singular integrand.' )
  print ( '' )
  print ( '        N        Estimate           Error' )
  print ( '' )
  seed = 123456789
  n = 17
  v2 = integral_log ( )
  for nlog in range ( 5, 21 ):
    n = ( n - 1 ) * 2 + 1
    h = 1.0 / float ( n )
    x, seed = r8vec_uniform_01 ( n, seed )
    f = integrand_log ( n, x )
    v1 = h * np.sum ( f )
    print ( '  %7d  %14.6g  %14.6g' % ( n, v1, abs ( v1 - v2 ) ) )
  print ( '' )
  print ( '    Exact: %14.6g' % ( v2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'MONTE_CARLO_LOG_TEST' )
  print ( '  Normal end of execution.' )
  return

def monte_carlo_power_test ( ):

#*****************************************************************************80
#
## MONTE_CARLO_POWER_TEST tests the Monte Carlo rule on the power singular integrand.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 November 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from test_functions import integral_power
  from test_functions import integrand_power
  from r8vec_uniform_01 import r8vec_uniform_01

  print ( '' )
  print ( 'MONTE_CARLO_POWER_TEST' )
  print ( '  Test the Monte Carlo rule on the power singular integrand.' )
  print ( '' )
  print ( '        N        Estimate           Error' )
  print ( '' )
  seed = 123456789
  n = 17
  v2 = integral_power ( )
  for nlog in range ( 5, 21 ):
    n = ( n - 1 ) * 2 + 1
    h = 1.0 / float ( n )
    x, seed = r8vec_uniform_01 ( n, seed )
    f = integrand_power ( n, x )
    v1 = h * np.sum ( f )
    print ( '  %7d  %14.6g  %14.6g' % ( n, v1, abs ( v1 - v2 ) ) )
  print ( '' )
  print ( '    Exact: %14.6g' % ( v2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'MONTE_CARLO_POWER_TEST' )
  print ( '  Normal end of execution.' )
  return

def monte_carlo_regular_test ( ):

#*****************************************************************************80
#
## MONTE_CARLO_REGULAR_TEST tests the Monte Carlo rule on the regular integrand.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 November 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from test_functions import integral_regular
  from test_functions import integrand_regular
  from r8vec_uniform_01 import r8vec_uniform_01

  print ( '' )
  print ( 'MONTE_CARLO_REGULAR_TEST' )
  print ( '  Test the Monte Carlo rule on the regular integrand.' )
  print ( '' )
  print ( '        N        Estimate           Error' )
  print ( '' )
  seed = 123456789
  n = 17
  v2 = integral_regular ( )
  for nlog in range ( 5, 21 ):
    n = ( n - 1 ) * 2 + 1
    h = 1.0 / float ( n )
    x, seed = r8vec_uniform_01 ( n, seed )
    f = integrand_regular ( n, x )
    v1 = h * sum ( f )
    print ( '  %7d  %14.6g  %14.6g' % ( n, v1, abs ( v1 - v2 ) ) )
  print ( '' )
  print ( '    Exact: %14.6g' % ( v2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'MONTE_CARLO_REGULAR_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( ) 
  monte_carlo_log_test ( )
  monte_carlo_power_test ( )
  monte_carlo_regular_test ( )
  timestamp ( )