#! /usr/bin/env python # def gamma_values ( n_data ): #*****************************************************************************80 # ## GAMMA_VALUES returns some values of the Gamma function. # # Discussion: # # The Gamma function is defined as: # # Gamma(Z) = Integral ( 0 <= T < Infinity) T^(Z-1) exp(-T) dT # # It satisfies the recursion: # # Gamma(X+1) = X * Gamma(X) # # Gamma is undefined for nonpositive integral X. # Gamma(0.5) = sqrt(PI) # For N a positive integer, Gamma(N+1) = N!, the standard factorial. # # In Mathematica, the function can be evaluated by: # # Gamma[x] # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 25 July 2014 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Parameters: # # Input/output, integer N_DATA. The user sets N_DATA to 0 before the # first call. On each call, the routine increments N_DATA by 1, and # returns the corresponding data; when there is no more data, the # output value of N_DATA will be 0 again. # # Output, real X, the argument of the function. # # Output, real FX, the value of the function. # import numpy as np n_max = 25 fx_vec = np.array ( ( \ -0.3544907701811032E+01, \ -0.1005871979644108E+03, \ 0.9943258511915060E+02, \ 0.9513507698668732E+01, \ 0.4590843711998803E+01, \ 0.2218159543757688E+01, \ 0.1772453850905516E+01, \ 0.1489192248812817E+01, \ 0.1164229713725303E+01, \ 0.1000000000000000E+01, \ 0.9513507698668732E+00, \ 0.9181687423997606E+00, \ 0.8974706963062772E+00, \ 0.8872638175030753E+00, \ 0.8862269254527580E+00, \ 0.8935153492876903E+00, \ 0.9086387328532904E+00, \ 0.9313837709802427E+00, \ 0.9617658319073874E+00, \ 0.1000000000000000E+01, \ 0.2000000000000000E+01, \ 0.6000000000000000E+01, \ 0.3628800000000000E+06, \ 0.1216451004088320E+18, \ 0.8841761993739702E+31 ) ) x_vec = np.array ( ( \ -0.50E+00, \ -0.01E+00, \ 0.01E+00, \ 0.10E+00, \ 0.20E+00, \ 0.40E+00, \ 0.50E+00, \ 0.60E+00, \ 0.80E+00, \ 1.00E+00, \ 1.10E+00, \ 1.20E+00, \ 1.30E+00, \ 1.40E+00, \ 1.50E+00, \ 1.60E+00, \ 1.70E+00, \ 1.80E+00, \ 1.90E+00, \ 2.00E+00, \ 3.00E+00, \ 4.00E+00, \ 10.00E+00, \ 20.00E+00, \ 30.00E+00 ) ) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 x = 0.0 fx = 0.0 else: x = x_vec[n_data] fx = fx_vec[n_data] n_data = n_data + 1 return n_data, x, fx def gamma_values_test ( ): #*****************************************************************************80 # ## GAMMA_VALUES_TEST demonstrates the use of GAMMA_VALUES. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 02 February 2009 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'GAMMA_VALUES_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' GAMMA_VALUES stores values of the Gamma function.' ) print ( '' ) print ( ' X GAMMA(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx = gamma_values ( n_data ) if ( n_data == 0 ): break print ( ' %12f %24.16f' % ( x, fx ) ) # # Terminate. # print ( '' ) print ( 'GAMMA_VALUES_TEST:' ) print ( ' Normal end of execution.' ) return def hyperball_monte_carlo_test01 ( ): #*****************************************************************************80 # ## HYPERBALL_MONTE_CARLO_TEST01 uses HYPERBALL01_SAMPLE in 3D. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 10 November 2016 # # Author: # # John Burkardt # import numpy as np import platform m = 3 e_test = np.array ( [ \ [ 0, 0, 0 ], \ [ 2, 0, 0 ], \ [ 0, 2, 0 ], \ [ 0, 0, 2 ], \ [ 4, 0, 0 ], \ [ 2, 2, 0 ], \ [ 0, 0, 4 ] ] ) print ( '' ) print ( 'HYPERBALL_MONTE_CARLO_TEST01' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Use the Monte Carlo method to estimate integrals' ) print ( ' over the interior of the unit hyperball in M dimensions.' ) print ( '' ) print ( ' Spatial dimension M = %d' % ( m ) ) seed = 123456789 print ( '' ) print ( ' N 1 X^2 Y^2' ), print ( ' Z^2 X^4 X^2Y^2 Z^4' ) print ( '' ) n = 1 e = np.zeros ( 7 ) while ( n <= 65536 ): x, seed = hyperball01_sample ( m, n, seed ) print ( ' %8d' % ( n ) ), for j in range ( 0, 7 ): e[0:m] = e_test[j,0:m] value = monomial_value ( m, n, e, x ) result = hyperball01_volume ( m ) * np.sum ( value[0:n] ) / float ( n ) print ( ' %14.6g' % ( result ) ), print ( '' ) n = 2 * n print ( '' ) print ( ' Exact' ), for j in range ( 0, 7 ): e[0:m] = e_test[j,0:m] result = hyperball01_monomial_integral ( m, e ) print ( ' %14.6g' % ( result ) ), print ( '' ) # # Terminate. # print ( '' ) print ( 'HYPERBALL_MONTE_CARLO_TEST01' ) print ( ' Normal end of execution.' ) return def hyperball_monte_carlo_test02 ( ): #*****************************************************************************80 # ## HYPERBALL_MONTE_CARLO_TEST02 uses HYPERBALL01_SAMPLE in 6D. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 10 November 2016 # # Author: # # John Burkardt # import numpy as np import platform m = 6 e_test = np.array ( [ \ [ 0, 0, 0, 0, 0, 0 ], \ [ 1, 0, 0, 0, 0, 0 ], \ [ 0, 2, 0, 0, 0, 0 ], \ [ 0, 2, 2, 0, 0, 0 ], \ [ 0, 0, 0, 4, 0, 0 ], \ [ 2, 0, 0, 0, 2, 2 ], \ [ 0, 0, 0, 0, 0, 6 ] ] ) print ( '' ) print ( 'HYPERBALL_MONTE_CARLO_TEST02' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Use the Monte Carlo method to estimate integrals' ) print ( ' over the interior of the unit hyperball in M dimensions.' ) print ( '' ) print ( ' Spatial dimension M = %d' % ( m ) ) seed = 123456789 print ( '' ) print ( ' N' ), print ( ' 1 ' ), print ( ' U ' ), print ( ' V^2 ' ), print ( ' V^2W^2' ), print ( ' X^4 ' ), print ( ' Y^2Z^2' ), print ( ' Z^6\n' ), print ( '' ) n = 1 e = np.zeros ( m ) while ( n <= 65536 ): x, seed = hyperball01_sample ( m, n, seed ) print ( ' %8d' % ( n ) ), for j in range ( 0, 7 ): e[0:m] = e_test[j,0:m] value = monomial_value ( m, n, e, x ) result = hyperball01_volume ( m ) * np.sum ( value[0:n] ) / float ( n ) print ( ' %14.6g' % ( result ) ), print ( '' ) n = 2 * n print ( '' ) print ( ' Exact' ), for j in range ( 0, 7 ): e[0:m] = e_test[j,0:m] result = hyperball01_monomial_integral ( m, e ) print ( ' %14.6g' % ( result ) ), print ( '' ) # # Terminate. # print ( '' ) print ( 'HYPERBALL_MONTE_CARLO_TEST02' ) print ( ' Use the Monte Carlo method to estimate integrals' ) return def hyperball01_monomial_integral ( m, e ): #*****************************************************************************80 # ## HYPERBALL01_MONOMIAL_INTEGRAL: integrals in unit hyperball in M dimensions. # # 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 # # Reference: # # Gerald Folland, # How to Integrate a Polynomial Over a Sphere, # American Mathematical Monthly, # Volume 108, Number 5, May 2001, pages 446-448. # # Parameters: # # Input, integer M, the spatial dimension. # # Input, integer E(M), the exponents. Each exponent must be nonnegative. # # Output, real INTEGRAL, the integral. # from sys import exit for i in range ( 0, m ): if ( e[i] < 0 ): print ( '' ) print ( 'HYPERBALL01_MONOMIAL_INTEGRAL - Fatal error!' ) print ( ' All exponents must be nonnegative.' ) exit ( 'HYPERBALL01_MONOMIAL_INTEGRAL - Fatal error!' ) # # Integrate over the surface. # 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 ): integral = integral * r8_gamma ( 0.5 * float ( e[i] + 1 ) ) s = 0 for i in range ( 0, m ): s = s + e[i] + 1 integral = integral / r8_gamma ( 0.5 * float ( s ) ) # # The surface integral is now adjusted to give the volume integral. # r = 1.0 integral = integral * r ** s / float ( s ) return integral def hyperball01_monomial_integral_test ( ): #*****************************************************************************80 # ## HYPERBALL01_MONOMIAL_INTEGRAL_TEST tests HYPERBALL01_MONOMIAL_INTEGRAL. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 22 June 2015 # # Author: # # John Burkardt # import numpy as np import platform m = 3 n = 4192 test_num = 20 print ( '' ) print ( 'HYPERBALL01_MONOMIAL_INTEGRAL_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' HYPERBALL01_MONOMIAL_INTEGRAL computes the integral of a monomial' ) print ( ' over the interior of the unit hyperball in M dimensions.' ) print ( ' Compare with a Monte Carlo estimate.' ) print ( '' ) print ( ' Spatial dimension M = %d' % ( m ) ) # # Get sample points. # seed = 123456789 x, seed = hyperball01_sample ( m, n, seed ) print ( '' ) print ( ' Number of sample points used is %d' % ( n ) ) # # Randomly choose exponents between 0 and 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 = hyperball01_volume ( m ) * np.sum ( value ) / float ( n ) exact = hyperball01_monomial_integral ( m, e ) error = abs ( result - exact ) for i in range ( 0, m ): print ( ' %2d' % ( e[i] ) ), print ( ' %14.6g %14.6g %10.2g' % ( result, exact, error ) ) # # Terminate. # print ( '' ) print ( 'HYPERBALL01_MONOMIAL_INTEGRAL_TEST' ) print ( ' Normal end of execution.' ) return def hyperball01_sample ( m, n, seed ): #*****************************************************************************80 # ## HYPERBALL01_SAMPLE uniformly samples the unit hyperball in M dimensions. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 22 June 2015 # # Author: # # John Burkardt # # Reference: # # Russell Cheng, # Random Variate Generation, # in Handbook of Simulation, # edited by Jerry Banks, # Wiley, 1998, pages 168. # # Reuven Rubinstein, # Monte Carlo Optimization, Simulation, and Sensitivity # of Queueing Networks, # Krieger, 1992, # ISBN: 0894647644, # LC: QA298.R79. # # Parameters: # # Input, integer M, the spatial dimension. # # Input, integer N, the number of points. # # Input/output, integer SEED, a seed for the random # number generator. # # Output, real X(M,N), the points. # import numpy as np x = np.zeros ( [ m, n ] ) exponent = 1.0 / float ( m ) for j in range ( 0, n ): # # Fill a vector with normally distributed values. # v, seed = r8vec_normal_01 ( m, seed ) # # Compute the length of the vector. # norm = r8vec_norm ( m, v ) # # Normalize the vector. # for i in range ( 0, m ): v[i] = v[i] / norm # # Now compute a value to map the point ON the sphere INTO the sphere. # r, seed = r8_uniform_01 ( seed ) for i in range ( 0, m ): x[i,j] = r ** exponent * v[i] return x, seed def hyperball01_sample_test ( ): #*****************************************************************************80 # ## HYPERBALL01_SAMPLE_TEST tests HYPERBALL01_SAMPLE. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 22 June 2015 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'HYPERBALL01_SAMPLE_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' HYPERBALL01_SAMPLE samples the unit hyperball.' ) m = 3 n = 10 seed = 123456789 x, seed = hyperball01_sample ( m, n, seed ) r8mat_transpose_print ( m, n, x, ' Sample points in the unit hyperball.' ) # # Terminate. # print ( '' ) print ( 'HYPERBALL01_SAMPLE_TEST' ) print ( ' Normal end of execution.' ) return def hyperball01_volume ( m ): #*****************************************************************************80 # ## HYPERBALL01_VOLUME returns the volume of the unit hyperball in M dimensions. # # Discussion: # # M Volume # # 1 2 # 2 1 * PI # 3 ( 4 / 3) * PI # 4 ( 1 / 2) * PI^2 # 5 ( 8 / 15) * PI^2 # 6 ( 1 / 6) * PI^3 # 7 (16 / 105) * PI^3 # 8 ( 1 / 24) * PI^4 # 9 (32 / 945) * PI^4 # 10 ( 1 / 120) * PI^5 # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 22 June 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the spatial dimension. # # Output, real VOLUME, the volume of the unit ball. # import numpy as np if ( ( m % 2 ) == 0 ): m_half = ( m // 2 ); volume = np.pi ** m_half for i in range ( 1, m_half + 1 ): volume = volume / float ( i ) else: m_half = ( ( m - 1 ) // 2 ) volume = np.pi ** m_half * 2.0 ** m for i in range ( m_half + 1, 2 * m_half + 2 ): volume = volume / float ( i ) return volume def hyperball01_volume_test ( ): #*****************************************************************************80 # ## HYPERBALL01_VOLUME tests HYPERBALL01_VOLUME. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 22 June 2015 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'HYPERBALL01_VOLUME_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' HYPERBALL01_VOLUME returns the volume of the unit hyperball' ) print ( ' in M dimensions.' ) print ( '' ) print ( ' M Volume' ) print ( '' ) for m in range ( 1, 11 ): value = hyperball01_volume ( m ) print ( ' %2d %g' % ( m, value ) ) # # Terminate. # print ( '' ) print ( 'HYPERBALL01_VOLUME_TEST' ) print ( ' Normal end of execution.' ) return def i4vec_print ( n, a, title ): #*****************************************************************************80 # ## I4VEC_PRINT prints an I4VEC. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the dimension of the vector. # # Input, integer A(N), the vector to be printed. # # Input, string TITLE, a title. # print ( '' ) print ( title ) print ( '' ) for i in range ( 0, n ): print ( '%6d %6d' % ( i, a[i] ) ) return def i4vec_print_test ( ): #*****************************************************************************80 # ## I4VEC_PRINT_TEST tests I4VEC_PRINT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 25 September 2016 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'I4VEC_PRINT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' I4VEC_PRINT prints an I4VEC.' ) n = 4 v = np.array ( [ 91, 92, 93, 94 ], dtype = np.int32 ) i4vec_print ( n, v, ' Here is an I4VEC:' ) # # Terminate. # print ( '' ) print ( 'I4VEC_PRINT_TEST:' ) print ( ' Normal end of execution.' ) return def i4vec_transpose_print ( n, a, title ): #*****************************************************************************80 # ## I4VEC_TRANSPOSE_PRINT prints an I4VEC "transposed". # # Example: # # A = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 /) # TITLE = 'My vector: ' # # My vector: # # 1 2 3 4 5 # 6 7 8 9 10 # 11 # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 02 June 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the number of components of the vector. # # Input, integer A(N), the vector to be printed. # # Input, string TITLE, a title. # if ( 0 < len ( title ) ): print ( '' ) print ( title ) if ( 0 < n ): for i in range ( 0, n ): print ( '%8d' % ( a[i] ) ), if ( ( i + 1 ) % 10 == 0 or i == n - 1 ): print ( '' ) else: print ( ' (empty vector)' ) return def i4vec_transpose_print_test ( ): #*****************************************************************************80 # ## I4VEC_TRANSPOSE_PRINT_TEST tests I4VEC_TRANSPOSE_PRINT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 07 April 2015 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'I4VEC_TRANSPOSE_PRINT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' I4VEC_TRANSPOSE_PRINT prints an I4VEC' ) print ( ' with 5 entries to a row, and an optional title.' ) n = 12 a = np.zeros ( n, dtype = np.int32 ) for i in range ( 0, n ): a[i] = i + 1 i4vec_transpose_print ( n, a, ' My array: ' ) # # Terminate. # print ( '' ) print ( 'I4VEC_TRANSPOSE_PRINT_TEST:' ) print ( ' Normal end of execution.' ) return def i4vec_uniform_ab ( n, a, b, seed ): #*****************************************************************************80 # ## I4VEC_UNIFORM_AB returns a scaled pseudorandom I4VEC. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 April 2013 # # Author: # # John Burkardt # # Reference: # # Paul Bratley, Bennett Fox, Linus Schrage, # A Guide to Simulation, # Second Edition, # Springer, 1987, # ISBN: 0387964673, # LC: QA76.9.C65.B73. # # Bennett Fox, # Algorithm 647: # Implementation and Relative Efficiency of Quasirandom # Sequence Generators, # ACM Transactions on Mathematical Software, # Volume 12, Number 4, December 1986, pages 362-376. # # Pierre L'Ecuyer, # Random Number Generation, # in Handbook of Simulation, # edited by Jerry Banks, # Wiley, 1998, # ISBN: 0471134031, # LC: T57.62.H37. # # Peter Lewis, Allen Goodman, James Miller, # A Pseudo-Random Number Generator for the System/360, # IBM Systems Journal, # Volume 8, Number 2, 1969, pages 136-143. # # Parameters: # # Input, integer N, the number of entries in the vector. # # Input, integer A, B, the minimum and maximum acceptable values. # # Input, integer SEED, a seed for the random number generator. # # Output, integer C(N), the randomly chosen integer vector. # # Output, integer SEED, the updated seed. # import numpy as np from sys import exit i4_huge = 2147483647 seed = int ( seed ) if ( seed < 0 ): seed = seed + i4_huge if ( seed == 0 ): print ( '' ) print ( 'I4VEC_UNIFORM_AB - Fatal error!' ) print ( ' Input SEED = 0!' ) exit ( 'I4VEC_UNIFORM_AB - Fatal error!' ) a = round ( a ) b = round ( b ) c = np.zeros ( n, dtype = np.int32 ) for i in range ( 0, n ): k = ( seed // 127773 ) seed = 16807 * ( seed - k * 127773 ) - k * 2836 seed = ( seed % i4_huge ) if ( seed < 0 ): seed = seed + i4_huge r = seed * 4.656612875E-10 # # Scale R to lie between A-0.5 and B+0.5. # r = ( 1.0 - r ) * ( min ( a, b ) - 0.5 ) \ + r * ( max ( a, b ) + 0.5 ) # # Use rounding to convert R to an integer between A and B. # value = round ( r ) value = max ( value, min ( a, b ) ) value = min ( value, max ( a, b ) ) c[i] = value return c, seed def i4vec_uniform_ab_test ( ): #*****************************************************************************80 # ## I4VEC_UNIFORM_AB_TEST tests I4VEC_UNIFORM_AB. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 27 October 2014 # # Author: # # John Burkardt # import platform n = 20 a = -100 b = 200 seed = 123456789 print ( '' ) print ( 'I4VEC_UNIFORM_AB_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' I4VEC_UNIFORM_AB computes pseudorandom values' ) print ( ' in an interval [A,B].' ) print ( '' ) print ( ' The lower endpoint A = %d' % ( a ) ) print ( ' The upper endpoint B = %d' % ( b ) ) print ( ' The initial seed is %d' % ( seed ) ) print ( '' ) v, seed = i4vec_uniform_ab ( n, a, b, seed ) i4vec_print ( n, v, ' The random vector:' ) # # Terminate. # print ( '' ) print ( 'I4VEC_UNIFORM_AB_TEST:' ) print ( ' Normal end of execution.' ) return def monomial_value ( m, n, e, x ): #*****************************************************************************80 # ## MONOMIAL_VALUE evaluates a monomial. # # Discussion: # # This routine evaluates a monomial of the form # # product ( 1 <= i <= m ) x(i)^e(i) # # The combination 0.0^0, if encountered, is treated as 1.0. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 07 April 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the spatial dimension. # # Input, integer N, the number of evaluation points. # # Input, integer E(M), the exponents. # # Input, real X(M,N), the point coordinates. # # Output, real V(N), the monomial values. # import numpy as np v = np.ones ( n ) for i in range ( 0, m ): if ( 0 != e[i] ): for j in range ( 0, n ): v[j] = v[j] * x[i,j] ** e[i] return v def monomial_value_test ( ): #*****************************************************************************80 # ## MONOMIAL_VALUE_TEST tests MONOMIAL_VALUE on sets of data in various dimensions. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 07 April 2015 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'MONOMIAL_VALUE_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Use monomial_value() to evaluate some monomials' ) print ( ' in dimensions 1 through 3.' ) e_min = -3 e_max = 6 n = 5 seed = 123456789 x_min = -2.0 x_max = +10.0 for m in range ( 1, 4 ): print ( '' ) print ( ' Spatial dimension M = %d' % ( m ) ) e, seed = i4vec_uniform_ab ( m, e_min, e_max, seed ) i4vec_transpose_print ( m, e, ' Exponents:' ) x, seed = r8mat_uniform_ab ( m, n, x_min, x_max, seed ) # # To make checking easier, make the X values integers. # for i in range ( 0, m ): for j in range ( 0, n ): x[i,j] = round ( x[i,j] ) v = monomial_value ( m, n, e, x ) print ( '' ) print ( ' V(X) ' ), for i in range ( 0, m ): print ( ' X(%d)' % ( i ) ), print ( '' ) print ( '' ) for j in range ( 0, n ): print ( '%14.6g ' % ( v[j] ) ), for i in range ( 0, m ): print ( '%10.4f' % ( x[i,j] ) ), print ( '' ) # # Terminate. # print ( '' ) print ( 'MONOMIAL_VALUE_TEST' ) print ( ' Normal end of execution.' ) return def r8_gamma ( x ): #*****************************************************************************80 # ## R8_GAMMA evaluates Gamma(X) for a real argument. # # Discussion: # # This routine calculates the gamma function for a real argument X. # # Computation is based on an algorithm outlined in reference 1. # The program uses rational functions that approximate the gamma # function to at least 20 significant decimal digits. Coefficients # for the approximation over the interval (1,2) are unpublished. # Those for the approximation for 12 <= X are from reference 2. # # PYTHON provides a GAMMA function, which is likely to be faster, and more # accurate. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 24 July 2014 # # Author: # # Original FORTRAN77 version by William Cody, Laura Stoltz. # PYTHON version by John Burkardt. # # Reference: # # William Cody, # An Overview of Software Development for Special Functions, # in Numerical Analysis Dundee, 1975, # edited by GA Watson, # Lecture Notes in Mathematics 506, # Springer, 1976. # # John Hart, Ward Cheney, Charles Lawson, Hans Maehly, # Charles Mesztenyi, John Rice, Henry Thatcher, # Christoph Witzgall, # Computer Approximations, # Wiley, 1968, # LC: QA297.C64. # # Parameters: # # Input, real X, the argument of the function. # # Output, real VALUE, the value of the function. # import numpy as np # # Coefficients for minimax approximation over (12, INF). # c = np.array ( [ -1.910444077728E-03, \ 8.4171387781295E-04, \ -5.952379913043012E-04, \ 7.93650793500350248E-04, \ -2.777777777777681622553E-03, \ 8.333333333333333331554247E-02, \ 5.7083835261E-03 ] ) # # Mathematical constants # sqrtpi = 0.9189385332046727417803297 # # Machine dependent parameters # xbig = 171.624 xminin = 2.23E-308 eps = 2.22E-16 xinf = 1.79E+308 # # Numerator and denominator coefficients for rational minimax # approximation over (1,2). # p = np.array ( [ \ -1.71618513886549492533811E+00, \ 2.47656508055759199108314E+01, \ -3.79804256470945635097577E+02, \ 6.29331155312818442661052E+02, \ 8.66966202790413211295064E+02, \ -3.14512729688483675254357E+04, \ -3.61444134186911729807069E+04, \ 6.64561438202405440627855E+04 ] ) q = np.array ( [ \ -3.08402300119738975254353E+01, \ 3.15350626979604161529144E+02, \ -1.01515636749021914166146E+03, \ -3.10777167157231109440444E+03, \ 2.25381184209801510330112E+04, \ 4.75584627752788110767815E+03, \ -1.34659959864969306392456E+05, \ -1.15132259675553483497211E+05 ] ) parity = 0 fact = 1.0 n = 0 y = x # # Argument is negative. # if ( y <= 0.0 ): y = - x y1 = int ( y ) res = y - y1 if ( res != 0.0 ): if ( y1 != int ( y1 * 0.5 ) * 2.0 ): parity = 1 fact = - np.pi / np.sin ( np.pi * res ) y = y + 1.0 else: res = xinf value = res return value # # Argument is positive. # if ( y < eps ): # # Argument < EPS. # if ( xminin <= y ): res = 1.0 / y else: res = xinf value = res return value elif ( y < 12.0 ): y1 = y # # 0.0 < argument < 1.0. # if ( y < 1.0 ): z = y y = y + 1.0 # # 1.0 < argument < 12.0. # Reduce argument if necessary. # else: n = int ( y ) - 1 y = y - n z = y - 1.0 # # Evaluate approximation for 1.0 < argument < 2.0. # xnum = 0.0 xden = 1.0 for i in range ( 0, 8 ): xnum = ( xnum + p[i] ) * z xden = xden * z + q[i] res = xnum / xden + 1.0 # # Adjust result for case 0.0 < argument < 1.0. # if ( y1 < y ): res = res / y1 # # Adjust result for case 2.0 < argument < 12.0. # elif ( y < y1 ): for i in range ( 0, n ): res = res * y y = y + 1.0 else: # # Evaluate for 12.0 <= argument. # if ( y <= xbig ): ysq = y * y sum = c[6] for i in range ( 0, 6 ): sum = sum / ysq + c[i] sum = sum / y - y + sqrtpi sum = sum + ( y - 0.5 ) * np.log ( y ) res = np.exp ( sum ) else: res = xinf value = res return value # # Final adjustments and return. # if ( parity ): res = - res if ( fact != 1.0 ): res = fact / res value = res return value def r8_gamma_test ( ): #*****************************************************************************80 # ## R8_GAMMA_TEST demonstrates the use of R8_GAMMA. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 25 July 2014 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'R8_GAMMA_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8_GAMMA evaluates the Gamma function.' ) print ( '' ) print ( ' X GAMMA(X) R8_GAMMA(X)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, fx1 = gamma_values ( n_data ) if ( n_data == 0 ): break fx2 = r8_gamma ( x ) print ( ' %12g %24.16g %24.16g' % ( x, fx1, fx2 ) ) # # Terminate. # print ( '' ) print ( 'R8_GAMMA_TEST' ) print ( ' Normal end of execution.' ) return def r8_normal_01 ( seed ): #*****************************************************************************80 # ## R8_NORMAL_01 samples the standard normal probability distribution. # # Discussion: # # The standard normal probability distribution function (PDF) has # mean 0 and standard deviation 1. # # The Box-Muller method is used. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 21 January 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer SEED, a seed for the random number generator. # # Output, real X, a sample of the standard normal PDF. # # Output, integer SEED, an updated seed for the random number generator. # import numpy as np r1, seed = r8_uniform_01 ( seed ) r2, seed = r8_uniform_01 ( seed ) x = np.sqrt ( - 2.0 * np.log ( r1 ) ) * np.cos ( 2.0 * np.pi * r2 ) return x, seed def r8_normal_01_test ( ): #*****************************************************************************80 # ## R8_NORMAL_01_TEST tests R8_NORMAL_01. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 26 July 2014 # # Author: # # John Burkardt # import platform seed = 123456789 test_num = 20 print ( '' ) print ( 'R8_NORMAL_01_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8_NORMAL_01 generates normally distributed' ) print ( ' random values.' ) print ( ' Using initial random number seed = %d' % ( seed ) ) print ( '' ) for test in range ( 0, test_num ): x, seed = r8_normal_01 ( seed ) print ( ' %f' % ( x ) ) # # Terminate. # print ( '' ) print ( 'R8_NORMAL_01_TEST' ) print ( ' Normal end of execution.' ) return def r8_uniform_01 ( seed ): #*****************************************************************************80 # ## R8_UNIFORM_01 returns a unit pseudorandom R8. # # Discussion: # # This routine implements the recursion # # seed = 16807 * seed mod ( 2^31 - 1 ) # r8_uniform_01 = seed / ( 2^31 - 1 ) # # The integer arithmetic never requires more than 32 bits, # including a sign bit. # # If the initial seed is 12345, then the first three computations are # # Input Output R8_UNIFORM_01 # SEED SEED # # 12345 207482415 0.096616 # 207482415 1790989824 0.833995 # 1790989824 2035175616 0.947702 # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 17 March 2013 # # Author: # # John Burkardt # # Reference: # # Paul Bratley, Bennett Fox, Linus Schrage, # A Guide to Simulation, # Second Edition, # Springer, 1987, # ISBN: 0387964673, # LC: QA76.9.C65.B73. # # Bennett Fox, # Algorithm 647: # Implementation and Relative Efficiency of Quasirandom # Sequence Generators, # ACM Transactions on Mathematical Software, # Volume 12, Number 4, December 1986, pages 362-376. # # Pierre L'Ecuyer, # Random Number Generation, # in Handbook of Simulation, # edited by Jerry Banks, # Wiley, 1998, # ISBN: 0471134031, # LC: T57.62.H37. # # Peter Lewis, Allen Goodman, James Miller, # A Pseudo-Random Number Generator for the System/360, # IBM Systems Journal, # Volume 8, Number 2, 1969, pages 136-143. # # Parameters: # # Input, integer SEED, the integer "seed" used to generate # the output random number. SEED should not be 0. # # Output, real R, a random value between 0 and 1. # # Output, integer SEED, the updated seed. This would # normally be used as the input seed on the next call. # from sys import exit i4_huge = 2147483647 seed = int ( seed ) seed = ( seed % i4_huge ) if ( seed < 0 ): seed = seed + i4_huge if ( seed == 0 ): print ( '' ) print ( 'R8_UNIFORM_01 - Fatal error!' ) print ( ' Input SEED = 0!' ) exit ( 'R8_UNIFORM_01 - Fatal error!' ) k = ( seed // 127773 ) seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ): seed = seed + i4_huge r = seed * 4.656612875E-10 return r, seed def r8_uniform_01_test ( ): #*****************************************************************************80 # ## R8_UNIFORM_01_TEST tests R8_UNIFORM_01. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 26 July 2014 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'R8_UNIFORM_01_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8_UNIFORM_01 produces a sequence of random values.' ) seed = 123456789 print ( '' ) print ( ' Using random seed %d' % ( seed ) ) print ( '' ) print ( ' SEED R8_UNIFORM_01(SEED)' ) print ( '' ) for i in range ( 0, 10 ): seed_old = seed x, seed = r8_uniform_01 ( seed ) print ( ' %12d %14f' % ( seed, x ) ) print ( '' ) print ( ' Verify that the sequence can be restarted.' ) print ( ' Set the seed back to its original value, and see that' ) print ( ' we generate the same sequence.' ) seed = 123456789 print ( '' ) print ( ' SEED R8_UNIFORM_01(SEED)' ) print ( '' ) for i in range ( 0, 10 ): seed_old = seed x, seed = r8_uniform_01 ( seed ) print ( ' %12d %14f' % ( seed, x ) ) # # Terminate. # print ( '' ) print ( 'R8_UNIFORM_01_TEST' ) print ( ' Normal end of execution.' ) return def r8mat_print ( m, n, a, title ): #*****************************************************************************80 # ## R8MAT_PRINT prints an R8MAT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows in A. # # Input, integer N, the number of columns in A. # # Input, real A(M,N), the matrix. # # Input, string TITLE, a title. # r8mat_print_some ( m, n, a, 0, 0, m - 1, n - 1, title ) return def r8mat_print_test ( ): #*****************************************************************************80 # ## R8MAT_PRINT_TEST tests R8MAT_PRINT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 10 February 2015 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'R8MAT_PRINT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8MAT_PRINT prints an R8MAT.' ) m = 4 n = 6 v = np.array ( [ \ [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ], [ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ], [ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ], [ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 ) r8mat_print ( m, n, v, ' Here is an R8MAT:' ) # # Terminate. # print ( '' ) print ( 'R8MAT_PRINT_TEST:' ) print ( ' Normal end of execution.' ) return def r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## R8MAT_PRINT_SOME prints out a portion of an R8MAT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 10 February 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, N, the number of rows and columns of the matrix. # # Input, real A(M,N), an M by N matrix to be printed. # # Input, integer ILO, JLO, the first row and column to print. # # Input, integer IHI, JHI, the last row and column to print. # # Input, string TITLE, a title. # incx = 5 print ( '' ) print ( title ) if ( m <= 0 or n <= 0 ): print ( '' ) print ( ' (None)' ) return for j2lo in range ( max ( jlo, 0 ), min ( jhi + 1, n ), incx ): j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) print ( '' ) print ( ' Col: ' ), for j in range ( j2lo, j2hi + 1 ): print ( '%7d ' % ( j ) ), print ( '' ) print ( ' Row' ) i2lo = max ( ilo, 0 ) i2hi = min ( ihi, m ) for i in range ( i2lo, i2hi + 1 ): print ( '%7d :' % ( i ) ), for j in range ( j2lo, j2hi + 1 ): print ( '%12g ' % ( a[i,j] ) ), print ( '' ) return def r8mat_print_some_test ( ): #*****************************************************************************80 # ## R8MAT_PRINT_SOME_TEST tests R8MAT_PRINT_SOME. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 October 2014 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'R8MAT_PRINT_SOME_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8MAT_PRINT_SOME prints some of an R8MAT.' ) m = 4 n = 6 v = np.array ( [ \ [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ], [ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ], [ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ], [ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 ) r8mat_print_some ( m, n, v, 0, 3, 2, 5, ' Here is an R8MAT:' ) # # Terminate. # print ( '' ) print ( 'R8MAT_PRINT_SOME_TEST:' ) print ( ' Normal end of execution.' ) return def r8mat_transpose_print ( m, n, a, title ): #*****************************************************************************80 # ## R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, the number of rows in A. # # Input, integer N, the number of columns in A. # # Input, real A(M,N), the matrix. # # Input, string TITLE, a title. # r8mat_transpose_print_some ( m, n, a, 0, 0, m - 1, n - 1, title ) return def r8mat_transpose_print_test ( ): #*****************************************************************************80 # ## R8MAT_TRANSPOSE_PRINT_TEST tests R8MAT_TRANSPOSE_PRINT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 October 2014 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'R8MAT_TRANSPOSE_PRINT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8MAT_TRANSPOSE_PRINT prints an R8MAT.' ) m = 4 n = 3 v = np.array ( [ \ [ 11.0, 12.0, 13.0 ], [ 21.0, 22.0, 23.0 ], [ 31.0, 32.0, 33.0 ], [ 41.0, 42.0, 43.0 ] ], dtype = np.float64 ) r8mat_transpose_print ( m, n, v, ' Here is an R8MAT, transposed:' ) # # Terminate. # print ( '' ) print ( 'R8MAT_TRANSPOSE_PRINT_TEST:' ) print ( ' Normal end of execution.' ) return def r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ): #*****************************************************************************80 # ## R8MAT_TRANSPOSE_PRINT_SOME prints a portion of an R8MAT, transposed. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 13 November 2014 # # Author: # # John Burkardt # # Parameters: # # Input, integer M, N, the number of rows and columns of the matrix. # # Input, real A(M,N), an M by N matrix to be printed. # # Input, integer ILO, JLO, the first row and column to print. # # Input, integer IHI, JHI, the last row and column to print. # # Input, string TITLE, a title. # incx = 5 print ( '' ) print ( title ) if ( m <= 0 or n <= 0 ): print ( '' ) print ( ' (None)' ) return for i2lo in range ( max ( ilo, 0 ), min ( ihi, m - 1 ), incx ): i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m - 1 ) i2hi = min ( i2hi, ihi ) print ( '' ) print ( ' Row: ' ), for i in range ( i2lo, i2hi + 1 ): print ( '%7d ' % ( i ) ), print ( '' ) print ( ' Col' ) j2lo = max ( jlo, 0 ) j2hi = min ( jhi, n - 1 ) for j in range ( j2lo, j2hi + 1 ): print ( '%7d :' % ( j ) ), for i in range ( i2lo, i2hi + 1 ): print ( '%12g ' % ( a[i,j] ) ), print ( '' ) return def r8mat_transpose_print_some_test ( ): #*****************************************************************************80 # ## R8MAT_TRANSPOSE_PRINT_SOME_TEST tests R8MAT_TRANSPOSE_PRINT_SOME. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 October 2014 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'R8MAT_TRANSPOSE_PRINT_SOME_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.' ) m = 4 n = 6 v = np.array ( [ \ [ 11.0, 12.0, 13.0, 14.0, 15.0, 16.0 ], [ 21.0, 22.0, 23.0, 24.0, 25.0, 26.0 ], [ 31.0, 32.0, 33.0, 34.0, 35.0, 36.0 ], [ 41.0, 42.0, 43.0, 44.0, 45.0, 46.0 ] ], dtype = np.float64 ) r8mat_transpose_print_some ( m, n, v, 0, 3, 2, 5, ' R8MAT, rows 0:2, cols 3:5:' ) # # Terminate. # print ( '' ) print ( 'R8MAT_TRANSPOSE_PRINT_SOME_TEST:' ) print ( ' Normal end of execution.' ) return def r8mat_uniform_01 ( m, n, seed ): #*****************************************************************************80 # ## R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 08 April 2013 # # Author: # # John Burkardt # # Reference: # # Paul Bratley, Bennett Fox, Linus Schrage, # A Guide to Simulation, # Second Edition, # Springer, 1987, # ISBN: 0387964673, # LC: QA76.9.C65.B73. # # Bennett Fox, # Algorithm 647: # Implementation and Relative Efficiency of Quasirandom # Sequence Generators, # ACM Transactions on Mathematical Software, # Volume 12, Number 4, December 1986, pages 362-376. # # Pierre L'Ecuyer, # Random Number Generation, # in Handbook of Simulation, # edited by Jerry Banks, # Wiley, 1998, # ISBN: 0471134031, # LC: T57.62.H37. # # Peter Lewis, Allen Goodman, James Miller, # A Pseudo-Random Number Generator for the System/360, # IBM Systems Journal, # Volume 8, Number 2, 1969, pages 136-143. # # Parameters: # # Input, integer M, N, the number of rows and columns in the array. # # Input, integer SEED, the integer "seed" used to generate # the output random number. # # Output, real R(M,N), an array of random values between 0 and 1. # # Output, integer SEED, the updated seed. This would # normally be used as the input seed on the next call. # import numpy from sys import exit i4_huge = 2147483647 seed = int ( seed ) if ( seed < 0 ): seed = seed + i4_huge if ( seed == 0 ): print ( '' ) print ( 'R8MAT_UNIFORM_01 - Fatal error!' ) print ( ' Input SEED = 0!' ) exit ( 'R8MAT_UNIFORM_01 - Fatal error!' ) r = numpy.zeros ( ( m, n ) ) for j in range ( 0, n ): for i in range ( 0, m ): k = ( seed // 127773 ) seed = 16807 * ( seed - k * 127773 ) - k * 2836 seed = ( seed % i4_huge ) if ( seed < 0 ): seed = seed + i4_huge r[i][j] = seed * 4.656612875E-10 return r, seed def r8mat_uniform_01_test ( ): #*****************************************************************************80 # ## R8MAT_UNIFORM_01_TEST tests R8MAT_UNIFORM_01. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 October 2014 # # Author: # # John Burkardt # import numpy as np import platform m = 5 n = 4 seed = 123456789 print ( '' ) print ( 'R8MAT_UNIFORM_01_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8MAT_UNIFORM_01 computes a random R8MAT.' ) print ( '' ) print ( ' 0 <= X <= 1' ) print ( ' Initial seed is %d' % ( seed ) ) v, seed = r8mat_uniform_01 ( m, n, seed ) r8mat_print ( m, n, v, ' Random R8MAT:' ) # # Terminate. # print ( '' ) print ( 'R8MAT_UNIFORM_01_TEST:' ) print ( ' Normal end of execution.' ) return def r8mat_uniform_ab ( m, n, a, b, seed ): #*****************************************************************************80 # ## R8MAT_UNIFORM_AB returns a scaled pseudorandom R8MAT. # # Discussion: # # An R8MAT is an array of R8's. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 08 April 2013 # # Author: # # John Burkardt # # Reference: # # Paul Bratley, Bennett Fox, Linus Schrage, # A Guide to Simulation, # Second Edition, # Springer, 1987, # ISBN: 0387964673, # LC: QA76.9.C65.B73. # # Bennett Fox, # Algorithm 647: # Implementation and Relative Efficiency of Quasirandom # Sequence Generators, # ACM Transactions on Mathematical Software, # Volume 12, Number 4, December 1986, pages 362-376. # # Pierre L'Ecuyer, # Random Number Generation, # in Handbook of Simulation, # edited by Jerry Banks, # Wiley, 1998, # ISBN: 0471134031, # LC: T57.62.H37. # # Peter Lewis, Allen Goodman, James Miller, # A Pseudo-Random Number Generator for the System/360, # IBM Systems Journal, # Volume 8, Number 2, 1969, pages 136-143. # # Parameters: # # Input, integer M, N, the number of rows and columns in the array. # # Input, real A, B, the range of the pseudorandom values. # # Input, integer SEED, the integer "seed" used to generate # the output random number. # # Output, real R(M,N), an array of random values between 0 and 1. # # Output, integer SEED, the updated seed. This would # normally be used as the input seed on the next call. # import numpy from sys import exit i4_huge = 2147483647 seed = int ( seed ) if ( seed < 0 ): seed = seed + i4_huge if ( seed == 0 ): print ( '' ) print ( 'R8MAT_UNIFORM_AB - Fatal error!' ) print ( ' Input SEED = 0!' ) exit ( 'R8MAT_UNIFORM_AB - Fatal error!' ) r = numpy.zeros ( ( m, n ) ) for j in range ( 0, n ): for i in range ( 0, m ): k = ( seed // 127773 ) seed = 16807 * ( seed - k * 127773 ) - k * 2836 seed = ( seed % i4_huge ) if ( seed < 0 ): seed = seed + i4_huge r[i][j] = a + ( b - a ) * seed * 4.656612875E-10 return r, seed def r8mat_uniform_ab_test ( ): #*****************************************************************************80 # ## R8MAT_UNIFORM_AB_TEST tests R8MAT_UNIFORM_AB. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 October 2014 # # Author: # # John Burkardt # import numpy as np m = 5 n = 4 a = -1.0 b = +5.0 seed = 123456789 print ( '' ) print ( 'R8MAT_UNIFORM_AB_TEST' ) print ( ' R8MAT_UNIFORM_AB computes a random R8MAT.' ) print ( '' ) print ( ' %g <= X <= %g' % ( a, b ) ) print ( ' Initial seed is %d' % ( seed ) ) v, seed = r8mat_uniform_ab ( m, n, a, b, seed ) r8mat_print ( m, n, v, ' Random R8MAT:' ) print ( '' ) print ( 'R8MAT_UNIFORM_AB_TEST:' ) print ( ' Normal end of execution.' ) return def r8vec_norm ( n, a ): #*****************************************************************************80 # ## R8VEC_NORM returns the L2 norm of an R8VEC. # # Discussion: # # The vector L2 norm is defined as: # # value = sqrt ( sum ( 1 <= I <= N ) A(I)^2 ). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 02 January 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the number of entries in A. # # Input, real A(N), the vector whose L2 norm is desired. # # Output, real VALUE, the L2 norm of A. # import numpy as np value = 0.0 for i in range ( 0, n ): value = value + a[i] * a[i] value = np.sqrt ( value ) return value def r8vec_norm_test ( ): #*****************************************************************************80 # ## R8VEC_NORM_TEST tests R8VEC_NORM. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 22 June 2015 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'R8VEC_NORM_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8VEC_NORM computes the L2 norm of an R8VEC.' ) n = 10 seed = 123456789 a, seed = r8vec_uniform_01 ( n, seed ) r8vec_print ( n, a, ' Input vector:' ) a_norm = r8vec_norm ( n, a ) print ( '' ) print ( ' L2 norm = %g' % ( a_norm ) ) # # Terminate. # print ( '' ) print ( 'R8VEC_NORM_TEST:' ) print ( ' Normal end of execution.' ) return def r8vec_normal_01 ( n, seed ): #*****************************************************************************80 # ## R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 04 March 2015 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the number of entries in the vector. # # Input, integer SEED, a seed for the random number generator. # # Output, real X(N), the vector of pseudorandom values. # # Output, integer SEED, an updated seed for the random number generator. # import numpy as np x = np.zeros ( n ) for i in range ( 0, n ): x[i], seed = r8_normal_01 ( seed ) return x, seed def r8vec_normal_01_test ( ): #*****************************************************************************80 # ## R8VEC_NORMAL_01_TEST tests R8VEC_NORMAL_01. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 04 March 2015 # # Author: # # John Burkardt # import numpy as np import platform n = 10 seed = 123456789 print ( '' ) print ( 'R8VEC_NORMAL_01_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8VEC_NORMAL_01 returns a vector of Normal 01 values' ) print ( '' ) print ( ' SEED = %d' % ( seed ) ) r, seed = r8vec_normal_01 ( n, seed ) r8vec_print ( n, r, ' Vector:' ) # # Terminate. # print ( '' ) print ( 'R8VEC_NORMAL_01_TEST:' ) print ( ' Normal end of execution.' ) return def r8vec_print ( n, a, title ): #*****************************************************************************80 # ## R8VEC_PRINT prints an R8VEC. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 31 August 2014 # # Author: # # John Burkardt # # Parameters: # # Input, integer N, the dimension of the vector. # # Input, real A(N), the vector to be printed. # # Input, string TITLE, a title. # print ( '' ) print ( title ) print ( '' ) for i in range ( 0, n ): print ( '%6d: %12g' % ( i, a[i] ) ) def r8vec_print_test ( ): #*****************************************************************************80 # ## R8VEC_PRINT_TEST tests R8VEC_PRINT. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 29 October 2014 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'R8VEC_PRINT_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8VEC_PRINT prints an R8VEC.' ) n = 4 v = np.array ( [ 123.456, 0.000005, -1.0E+06, 3.14159265 ], dtype = np.float64 ) r8vec_print ( n, v, ' Here is an R8VEC:' ) # # Terminate. # print ( '' ) print ( 'R8VEC_PRINT_TEST:' ) print ( ' Normal end of execution.' ) return def r8vec_uniform_01 ( n, seed ): #*****************************************************************************80 # ## R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # # Reference: # # Paul Bratley, Bennett Fox, Linus Schrage, # A Guide to Simulation, # Second Edition, # Springer, 1987, # ISBN: 0387964673, # LC: QA76.9.C65.B73. # # Bennett Fox, # Algorithm 647: # Implementation and Relative Efficiency of Quasirandom # Sequence Generators, # ACM Transactions on Mathematical Software, # Volume 12, Number 4, December 1986, pages 362-376. # # Pierre L'Ecuyer, # Random Number Generation, # in Handbook of Simulation, # edited by Jerry Banks, # Wiley, 1998, # ISBN: 0471134031, # LC: T57.62.H37. # # Peter Lewis, Allen Goodman, James Miller, # A Pseudo-Random Number Generator for the System/360, # IBM Systems Journal, # Volume 8, Number 2, 1969, pages 136-143. # # Parameters: # # Input, integer N, the number of entries in the vector. # # Input, integer SEED, a seed for the random number generator. # # Output, real X(N), the vector of pseudorandom values. # # Output, integer SEED, an updated seed for the random number generator. # import numpy from sys import exit i4_huge = 2147483647 seed = int ( seed ) if ( seed < 0 ): seed = seed + i4_huge if ( seed == 0 ): print ( '' ) print ( 'R8VEC_UNIFORM_01 - Fatal error!' ) print ( ' Input SEED = 0!' ) exit ( 'R8VEC_UNIFORM_01 - Fatal error!' ) x = numpy.zeros ( n ); for i in range ( 0, n ): k = ( seed // 127773 ) seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ): seed = seed + i4_huge x[i] = seed * 4.656612875E-10 return x, seed def r8vec_uniform_01_test ( ): #*****************************************************************************80 # ## R8VEC_UNIFORM_01_TEST tests R8VEC_UNIFORM_01. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 29 October 2014 # # Author: # # John Burkardt # import numpy as np import platform n = 10 seed = 123456789 print ( '' ) print ( 'R8VEC_UNIFORM_01_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' R8VEC_UNIFORM_01 computes a random R8VEC.' ) print ( '' ) print ( ' Initial seed is %d' % ( seed ) ) v, seed = r8vec_uniform_01 ( n, seed ) r8vec_print ( n, v, ' Random R8VEC:' ) # # Terminate. # print ( '' ) print ( 'R8VEC_UNIFORM_01_TEST:' ) print ( ' Normal end of execution.' ) return def timestamp ( ): #*****************************************************************************80 # ## TIMESTAMP prints the date as a timestamp. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # # Parameters: # # None # import time t = time.time ( ) print ( time.ctime ( t ) ) return None def timestamp_test ( ): #*****************************************************************************80 # ## TIMESTAMP_TEST tests TIMESTAMP. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 03 December 2014 # # Author: # # John Burkardt # # Parameters: # # None # import platform print ( '' ) print ( 'TIMESTAMP_TEST:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' TIMESTAMP prints a timestamp of the current date and time.' ) print ( '' ) timestamp ( ) # # Terminate. # print ( '' ) print ( 'TIMESTAMP_TEST:' ) print ( ' Normal end of execution.' ) return def hyperball_monte_carlo_test ( ): #*****************************************************************************80 # ## HYPERBALL_MONTE_CARLO_TEST tests the HYPERBALL_MONTE_CARLO library. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 January 2014 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'HYPERBALL_MONTE_CARLO_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Test the HYPERBALL_MONTE_CARLO library.' ) gamma_values_test ( ) hyperball_monte_carlo_test01 ( ) hyperball_monte_carlo_test02 ( ) hyperball01_monomial_integral_test ( ) hyperball01_sample_test ( ) hyperball01_volume_test ( ) i4vec_print_test ( ) i4vec_transpose_print_test ( ) i4vec_uniform_ab_test ( ) monomial_value_test ( ) r8_gamma_test ( ) r8_normal_01_test ( ) r8_uniform_01_test ( ) r8mat_print_test ( ) r8mat_print_some_test ( ) r8mat_transpose_print_test ( ) r8mat_transpose_print_some_test ( ) r8mat_uniform_01_test ( ) r8mat_uniform_ab_test ( ) r8vec_norm_test ( ) r8vec_normal_01_test ( ) r8vec_print_test ( ) r8vec_uniform_01_test ( ) # # Terminate. # print ( '' ) print ( 'HYPERBALL_MONTE_CARLO_TEST' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): timestamp ( ) hyperball_monte_carlo_test ( ) timestamp ( )