#! /usr/bin/env python
#
def hypersphere01_monomial_integral ( m, e ):

#*****************************************************************************80
#
## HYPERSPHERE01_MONOMIAL_INTEGRAL: monomial integrals on the unit hypersphere.
#
#  Discussion:
#
#    The integration region is 
#
#      sum ( 1 <= I <= M ) X(I)^2 = 1.
#
#    The monomial is F(X) = product ( 1 <= I <= M ) X(I)^E(I)
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    22 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer M, the spatial dimension.
#
#    Input, integer E(M), the exponents of X(1) through X(M). 
#    Each exponent must be nonnegative.
#
#    Output, real INTEGRAL, the integral.
#
  from r8_gamma import r8_gamma
  from sys import exit

  for i in range ( 0, m ):
    if ( e[i] < 0 ):
      print ( '' )
      print ( 'HYPERSPHERE01_MONOMIAL_INTEGRAL - Fatal error!' )
      print ( '  All exponents must be nonnegative.' )
      error ( 'HYPERSPHERE01_MONOMIAL_INTEGRAL - Fatal error!' )

  for i in range ( 0, m ):
    if ( ( e[i] % 2 ) == 1 ):
      integral = 0.0
      return integral

  integral = 2.0

  for i in range ( 0, m ):
    arg = 0.5 * float ( e[i] + 1 )
    integral = integral * r8_gamma ( arg )

  s = 0
  for i in range ( 0, m ):
    s = s + float ( e[i] + 1 )

  arg = 0.5 * float ( s )
  integral = integral / r8_gamma ( arg )

  return integral

def hypersphere01_monomial_integral_test ( ):

#*****************************************************************************80
#
## HYPERSPHERE01_MONOMIAL_INTEGRAL_TEST tests HYPERSPHERE01_MONOMIAL_INTEGRAL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    07 January 2014
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform
  from hypersphere01_area import hypersphere01_area
  from hypersphere01_sample import hypersphere01_sample
  from i4vec_uniform_ab import i4vec_uniform_ab
  from monomial_value import monomial_value

  m = 3
  n = 4192
  test_num = 20

  print ( '' )
  print ( 'HYPERSPHERE01_MONOMIAL_INTEGRAL' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  HYPERSPHERE01_MONOMIAL_INTEGRAL returns the integral of' )
  print ( '  a monomial over the surface of the unit hypersphere in 3D.' )
  print ( '  Compare with a Monte Carlo estimate.' )
#
#  Get sample points.
#
  seed = 123456789
  x, seed = hypersphere01_sample ( m, n, seed )

  print ( '' )
  print ( '  Number of sample points used is %d' % ( n ) )
#
#  Randomly choose X,Y,Z exponents between (0,0,0) and (8,8,8).
#
  print ( '' )
  print ( '  If any exponent is odd, the integral is zero.' )
  print ( '  We will restrict this test to randomly chosen even exponents.' )
  print ( '' )
  print ( '  Ex  Ey  Ez     MC-Estimate           Exact      Error' )
  print ( '' )

  for test in range ( 0, test_num ):

    e, seed = i4vec_uniform_ab ( m, 0, 4, seed )

    for i in range ( 0, m ):
      e[i] = e[i] * 2

    value = monomial_value ( m, n, e, x )

    result = hypersphere01_area ( m ) * np.sum ( value ) / float ( n )
    exact = hypersphere01_monomial_integral ( m, e )
    error = abs ( result - exact )

    for i in range ( 0, m ):
      print ( '  %2d' % ( e[i] ), end = '' )
    print ( '  %14.6g  %14.6g  %10.2g' % ( result, exact, error ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'HYPERSPHERE_MONOMIAL_INTEGRAL_TEST' )
  print ( '  Normal end of execution.' )
  return

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