# include # include # include # include # include using namespace std; # include "lattice_rule.hpp" //****************************************************************************80 double e_01_2d ( int dim_num, double a[], double b[] ) //****************************************************************************80 // // Purpose: // // E_01_2D is the exact integral of 2d test function #1. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 18 April 2003 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, double A[DIM_NUM], B[DIM_NUM], the integration limits. // // Output, double E_01_2D, the integral of the function // over the limits. // { double value; value = 1.0; return value; } //****************************************************************************80 double f_01_2d ( int dim_num, double x[] ) //****************************************************************************80 // // Purpose: // // F_01_2D is the 2D test function #1. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, double X[DIM_NUM], the point where the function // is to be evaluated. // // Output, double F_01_2D, the value of the function at X. // { double e = 2.718281828459045; double value; value = x[1] * exp ( x[0] * x[1] ) / ( e - 2.0 ); return value; } //****************************************************************************80 double f2 ( double x ) //****************************************************************************80 // // Purpose: // // F2 evaluates a function of a scalar used in defining P2(Q). // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, double X, the value of the argument. // // Output, double F2, the value of F2(X). // { double pi = 3.141592653589793; double value; value = 1.0 + 2.0 * pi * pi * ( x * x - x + 1.0 / 6.0 ); return value; } //****************************************************************************80 double f20_s ( int dim_num, double x[] ) //****************************************************************************80 // // Purpose: // // F20_S evaluates a function of a vector used in defining P2(Q). // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, double X[DIM_NUM], the value of the argument. // // Output, double F20_S, the value of F20_S(X). // { int i; double value; value = 1.0; for ( i = 0; i < dim_num; i++ ) { value = value * ( 1.0 + ( f2 ( x[i] ) - 1.0 ) ); } value = value - 1.0; return value; } //****************************************************************************80 int fibonacci ( int k ) //****************************************************************************80 // // Purpose: // // FIBONACCI returns the Fibonacci number of given index. // // Example: // // K Fibonacci // // 0 0 // 1 1 // 2 1 // 3 2 // 4 3 // 5 5 // 6 8 // 7 13 // 8 21 // 9 34 // 10 55 // 11 89 // 12 144 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 April 2003 // // Author: // // John Burkardt // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 1. // // Output, int FIBONACCI, the value of the K-th Fibonacci number. // { int a; int b; int c; int kk; if ( k < 0 ) { a = - i4_huge ( ); return a; } else if ( k == 0 ) { a = 0; return a; } else if ( k == 1 ) { a = 1; return a; } c = 0; b = 0; a = 1; for ( kk = 2; kk <= k; kk++ ) { c = b; b = a; a = c + b; } return a; } //****************************************************************************80 double fibonacci_lattice_b ( int k, double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_B applies an optimal Fibonacci lattice integration rule in 2D. // // Discussion: // // This routine may be applied to integrands which are not periodic. // // When K is odd, this is the same as the symmetric Fibonacci lattice // integration rule. But when K is even, a correction is made to the // corner weights which is expected to improve the results. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double FIBONACCI_LATTICE_B, the estimated integral. // { double delta; int dim; int dim_num = 2; int j; int m; int n; double quad; double quad1; double quad2; int rank; double w[2*2]; double *x; int *z; x = new double[dim_num]; z = new int[dim_num]; quad = 0.0; m = fibonacci ( k ); n = fibonacci ( k - 1 ); // // Get the corner weights. // if ( ( k % 2 ) == 1 ) { w[0+0*2] = 1.0 / ( double ) ( 4 * m ); w[1+0*2] = 1.0 / ( double ) ( 4 * m ); w[0+1*2] = 1.0 / ( double ) ( 4 * m ); w[1+1*2] = 1.0 / ( double ) ( 4 * m ); } else { delta = 0.0; for ( j = 1; j <= m - 1; j++ ) { delta = delta + ( double ) ( j * ( ( j * n ) % m ) ) / ( double ) ( m * m ); } w[0+0*2] = 0.25 - delta / ( double ) ( m ); delta = 0.0; for ( j = 1; j <= m - 1; j++ ) { delta = delta + ( double ) ( j * ( m - ( ( j * n ) % m ) ) ) / ( double ) ( m * m ); } w[0+1*2] = 0.25 - delta / ( double ) ( m ); w[1+0*2] = w[0+1*2]; w[1+1*2] = w[0+0*2]; } // // Get all the corner values. // rank = 0; quad1 = 0.0; for ( ; ; ) { tuple_next ( 0, 1, dim_num, &rank, z ); if ( rank == 0 ) { break; } for ( dim = 0; dim < dim_num; dim++ ) { x[dim] = ( double ) z[dim]; } quad1 = quad1 + w[z[0]+z[1]*2] * f ( dim_num, x ); } // // Get the interior values. // z[0] = 1; z[1] = fibonacci ( k - 1 ); quad2 = 0.0; for ( j = 1; j <= m - 1; j++ ) { for ( dim = 0; dim < dim_num; dim++ ) { x[dim] = fmod ( ( double ) ( j * z[dim] ) / ( double ) ( m ), 1.0 ); } quad2 = quad2 + f ( dim_num, x ); } quad = quad1 + quad2 / ( double ) ( m ); delete [] x; delete [] z; return quad; } //****************************************************************************80 double fibonacci_lattice_q ( int k, double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_Q applies a Fibonacci lattice integration rule in 2D. // // Discussion: // // Because this is a standard lattice rule, it is really only suited // for functions which are periodic, of period 1, in both X and Y. // // The related routines FIBONACCI_LATTICE_S and FIBONACCI_LATTICE_B // may be helpful in cases where the integrand does not satisfy this // requirement. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double FIBONACCI_LATTICE_Q, the estimated integral. // { int dim; int dim_num = 2; int j; int m; double quad; double *x; int *z; x = new double[dim_num]; z = new int[dim_num]; quad = 0.0; m = fibonacci ( k ); z[0] = 1; z[1] = fibonacci ( k - 1 ); for ( j = 0; j <= m - 1; j++ ) { for ( dim = 0; dim < dim_num; dim++ ) { x[dim] = fmod ( ( double ) ( j * z[dim] ) / ( double ) ( m ), 1.0 ); } quad = quad + f ( dim_num, x ); } quad = quad / ( double ) ( m ); delete [] z; delete [] x; return quad; } //****************************************************************************80 double *fibonacci_lattice_q_nodes ( int k ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_Q_NODES returns Fibonacci lattice nodes in 2D. // // Discussion: // // Because this is a standard lattice rule, it is really only suited // for functions which are periodic, of period 1, in both X and Y. // // The number of nodes returned is // // M = fibonacci ( k ). // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Output, double X[2*M], the nodes. // { int dim; int dim_num = 2; int j; int m; double *x; int *z; m = fibonacci ( k ); x = new double[2*m]; z = new int[dim_num]; z[0] = 1; z[1] = fibonacci ( k - 1 ); for ( j = 0; j <= m - 1; j++ ) { for ( dim = 0; dim < dim_num; dim++ ) { x[dim+j*dim_num] = fmod ( ( double ) ( j * z[dim] ) / ( double ) ( m ), 1.0 ); } } delete [] z; return x; } //****************************************************************************80 double fibonacci_lattice_q1 ( int k, double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_Q1 applies a Fibonacci lattice integration rule in 2D. // // Discussion: // // This is a modification of the algorithm in FIBONACCI_LATTICE_Q. // It uses a nonlinear transformation on the integrand, which makes // the lattice rule more suitable for nonperiodic integrands. // // The transformation replaces the integration variable X by // // PHI(X) = 3*X^2 - 2*X**3 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double FIBONACCI_LATTICE_Q1, the estimated integral. // { int dim_num = 2; double dphi; int i; int j; int m; double quad; double *x; double *y; int *z; x = new double[dim_num]; y = new double[dim_num]; z = new int[dim_num]; quad = 0.0; m = fibonacci ( k ); z[0] = 1; z[1] = fibonacci ( k - 1 ); for ( j = 0; j <= m - 1; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i] = fmod ( ( double ) ( j * z[i] ) / ( double ) ( m ), 1.0 ); } dphi = 1.0; for ( i = 0; i < dim_num; i++ ) { y[i] = ( 3.0 - 2.0 * x[i] ) * x[i] * x[i]; dphi = dphi * 6.0 * ( 1.0 - x[i] ) * x[i]; } quad = quad + f ( dim_num, y ) * dphi; } quad = quad / ( double ) ( m ); delete [] x; delete [] y; delete [] z; return quad; } //****************************************************************************80 double fibonacci_lattice_q2 ( int k, double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_Q2 applies a Fibonacci lattice integration rule in 2D. // // Discussion: // // This is a modification of the algorithm in FIBONACCI_LATTICE_Q. // It uses a nonlinear transformation on the integrand, which makes // the lattice rule more suitable for nonperiodic integrands. // // The transformation replaces the integration variable X by // // PHI(X) = 3*X^3 *( 10 - 15 * X + 6 * X^2 ) // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double FIBONACCI_LATTICE_Q2, the estimated integral. // { int dim_num = 2; double dphi; int i; int j; int m; double quad; double *x; double *y; int *z; x = new double[dim_num]; y = new double[dim_num]; z = new int[dim_num]; quad = 0.0; m = fibonacci ( k ); z[0] = 1; z[1] = fibonacci ( k - 1 ); for ( j = 0; j <= m - 1; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i] = fmod ( ( double ) ( j * z[i] ) / ( double ) ( m ), 1.0 ); } dphi = 1.0; for ( i = 0; i < dim_num; i++ ) { y[i] = ( 10.0 - 15.0 * x[i] + 6.0 * pow ( x[i], 2 ) ) * pow ( x[i], 3 ); dphi = dphi * 30.0 * pow ( 1.0 - x[i], 2 ) * pow ( x[i], 2 ); } quad = quad + f ( dim_num, y ) * dphi; } quad = quad / ( double ) ( m ); delete [] x; delete [] y; delete [] z; return quad; } //****************************************************************************80 double fibonacci_lattice_q3 ( int k, double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_Q3 applies a Fibonacci lattice integration rule in 2D. // // Discussion: // // This is a modification of the algorithm in FIBONACCI_LATTICE_Q. // It uses a nonlinear transformation on the integrand, which makes // the lattice rule more suitable for nonperiodic integrands. // // The transformation replaces the integration variable X by // // PHI(X) = X - sin ( 2 * PI * X ) / ( 2 * PI ) // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double FIBONACCI_LATTICE_Q3, the estimated integral. // { int dim_num = 2; double dphi; int i; int j; int m; double quad; double pi = 3.141592653589793; double two_pi; double *x; double *y; int *z; x = new double[dim_num]; y = new double[dim_num]; z = new int[dim_num]; quad = 0.0; m = fibonacci ( k ); z[0] = 1; z[1] = fibonacci ( k - 1 ); two_pi = 2.0 * pi; for ( j = 0; j <= m - 1; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i] = fmod ( ( double ) ( j * z[i] ) / ( double ) ( m ), 1.0 ); } dphi = 1.0; for ( i = 0; i < dim_num; i++ ) { y[i] = x[i] - sin ( two_pi * x[i] ) / two_pi; dphi = dphi * ( 1.0 - cos ( two_pi * x[i] ) ); } quad = quad + f ( dim_num, y ) * dphi; } quad = quad / ( double ) ( m ); delete [] x; delete [] y; delete [] z; return quad; } //****************************************************************************80 double fibonacci_lattice_t ( int k, double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // FIBONACCI_LATTICE_T applies a symmetric Fibonacci lattice integration rule in 2D. // // Discussion: // // This routine may be applied to integrands which are not periodic. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int K, the index of the Fibonacci number to be used. // K must be at least 3. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double FIBONACCI_LATTICE_T, the estimated integral. // { int dim_num = 2; int i; int j; int m; double quad; double quad1; double quad2; int rank; double w; double *x; int *z; x = new double[dim_num]; z = new int[dim_num]; quad = 0.0; m = fibonacci ( k ); // // Get all the corner values. // rank = 0; quad1 = 0.0; w = 1.0 / ( double ) i4_power ( 2, dim_num ); for ( ; ; ) { tuple_next ( 0, 1, dim_num, &rank, z ); if ( rank == 0 ) { break; } for ( i = 0; i < dim_num; i++ ) { x[i] = ( double ) z[i]; } quad1 = quad1 + w * f ( dim_num, x ); } // // Get the interior values. // z[0] = 1; z[1] = fibonacci ( k - 1 ); quad2 = 0.0; for ( j = 1; j <= m - 1; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i] = fmod ( ( double ) ( j * z[i] ) / ( double ) ( m ), 1.0 ); } quad2 = quad2 + f ( dim_num, x ); } quad = ( quad1 + quad2 ) / ( double ) ( m ); delete [] x; delete [] z; return quad; } //****************************************************************************80 int *find_z20 ( int dim_num, int m ) //****************************************************************************80 // // Purpose: // // FIND_Z20 finds the appropriate Z vector to minimize P2(QS). // // Discussion: // // For the method of good lattice points, a number of points M, and // a single generator vector Z is chosen. The integrand is assumed // to be periodic of period 1 in each argument, and is evaluated at // each of the points X(I)(1:DIM_NUM) = I * Z(1:DIM_NUM) / M, // for I = 0 to M-1. The integral is then approximated by the average // of these values. // // Assuming that DIM_NUM and M are known, and that the integrand is not // known beforehand, the accuracy of the method depends entirely // on the choice of Z. One method of choosing Z is to search for // the Z among all candidates which minimizes a particular quantity // P_ALPHA(QS). // // We only need to look at vectors Z of the form // (1, L, L^2, ..., L^(DIM_NUM-1)), // for L = 1 to M/2. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int M, the number of points to be used. // // Output, int FIND_Z20[DIM_NUM], the optimal vector. // { int dim; int i; double q0; double q0_min; int value; int *z; int *z_min; z = new int[dim_num]; z_min = new int[dim_num]; q0_min = HUGE_VAL; for ( i = 1; i <= m / 2; i++ ) { value = 1; for ( dim = 0; dim < dim_num; dim++ ) { z[dim] = value; value = ( value * i ) % m; } // // Use this Z and the lattice integral method Q0 of order M, // to approximate the integral of P2. // q0 = lattice ( dim_num, m, z, f20_s ); // // If this result is the smallest so far, save the corresponding Z. // if ( q0 < q0_min ) { q0_min = q0; for ( dim = 0; dim < dim_num; dim++ ) { z_min[dim] = z[dim]; } } } delete [] z; // // Return the best Z. // return z_min; } //****************************************************************************80 void gray_next ( int n, int *change ) //****************************************************************************80 // // Purpose: // // GRAY_NEXT generates the next Gray code by switching one item at a time. // // Discussion: // // On the first call only, the user must set CHANGE = -N. // This initializes the routine to the Gray code for N zeroes. // // Each time it is called thereafter, it returns in CHANGE the index // of the item to be switched in the Gray code. The sign of CHANGE // indicates whether the item is to be added or subtracted (or // whether the corresponding bit should become 1 or 0). When // CHANGE is equal to N+1 on output, all the Gray codes have been // generated. // // The routine has internal memory that is set up on call with // CHANGE = -N, and released on final return. // // Example: // // N CHANGE Subset in/out Binary Number // Interpretation Interpretation // 1 2 4 8 // -- --------- -------------- -------------- // // 4 -4 / 0 0 0 0 0 0 // // +1 1 0 0 0 1 // +2 1 1 0 0 3 // -1 0 1 0 0 2 // +3 0 1 1 0 6 // +1 1 1 1 0 7 // -2 1 0 1 0 5 // -1 0 0 1 0 4 // +4 0 0 1 1 12 // +1 1 0 1 1 13 // +2 1 1 1 1 15 // -1 0 1 1 1 14 // -3 0 1 0 1 10 // +1 1 1 0 1 11 // -2 1 0 0 1 9 // -1 0 0 0 1 8 // -4 0 0 0 0 0 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 14 May 2003 // // Author: // // John Burkardt // // Reference: // // Albert Nijenhuis, Herbert Wilf, // Combinatorial Algorithms for Computers and Calculators, // Second Edition, // Academic Press, 1978, // ISBN: 0-12-519260-6, // LC: QA164.N54. // // Parameters: // // Input, int N, the order of the total set from which // subsets will be drawn. // // Input/output, int *CHANGE. This item is used for input only // on the first call for a particular sequence of Gray codes, // at which time it must be set to -N. This corresponds to // all items being excluded, or all bits being 0, in the Gray code. // On output, CHANGE indicates which of the N items must be "changed", // and the sign indicates whether the item is to be added or removed // (or the bit is to become 1 or 0). Note that on return from the // first call, CHANGE is set to 0, indicating that we begin with // the empty set. // { static int *a = NULL; int i; static int k = 0; static int n_save = -1; if ( n <= 0 ) { cout << "\n"; cout << "GRAY_NEXT - Fatal error!\n"; cout << " Input value of N <= 0.\n"; exit ( 1 ); } if ( *change == -n ) { if ( a ) { delete [] a; } a = new int[n]; for ( i = 0; i < n; i++ ) { a[i] = 0; } n_save = n; k = 1; *change = 0; return; } if ( n != n_save ) { cout << "\n"; cout << "GRAY_NEXT - Fatal error!\n"; cout << " Input value of N has changed from definition value.\n"; exit ( 1 ); } // // First determine WHICH item is to be changed. // if ( ( k % 2 ) == 1 ) { *change = 1; } else { for ( i = 1; i <= n_save; i++ ) { if ( a[i-1] == 1 ) { *change = i + 1; break; } } } // // Take care of the terminal case CHANGE = N + 1. // if ( *change == n + 1 ) { *change = n; } // // Now determine HOW the item is to be changed. // if ( a[*change-1] == 0 ) { a[*change-1] = 1; } else if ( a[*change-1] == 1 ) { a[*change-1] = 0; *change = -( *change ); } // // Update the counter. // k = k + 1; // // If the output CHANGE = -N_SAVE, then we're done. // if ( *change == -n_save ) { delete [] a; a = NULL; n_save = 0; k = 0; } return; } //****************************************************************************80 int i4_huge ( ) //****************************************************************************80 // // Purpose: // // I4_HUGE returns a "huge" I4. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 16 May 2003 // // Author: // // John Burkardt // // Parameters: // // Output, int I4_HUGE, a "huge" I4. // { return 2147483647; } //****************************************************************************80 int i4_max ( int i1, int i2 ) //****************************************************************************80 // // Purpose: // // I4_MAX returns the maximum of two I4's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 13 October 1998 // // Author: // // John Burkardt // // Parameters: // // Input, int I1, I2, are two integers to be compared. // // Output, int I4_MAX, the larger of I1 and I2. // { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } //****************************************************************************80 int i4_min ( int i1, int i2 ) //****************************************************************************80 // // Purpose: // // I4_MIN returns the minimum of two I4's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 13 October 1998 // // Author: // // John Burkardt // // Parameters: // // Input, int I1, I2, two integers to be compared. // // Output, int I4_MIN, the smaller of I1 and I2. // { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } //****************************************************************************80 int i4_power ( int i, int j ) //****************************************************************************80 // // Purpose: // // I4_POWER returns the value of I^J. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 01 April 2004 // // Author: // // John Burkardt // // Parameters: // // Input, int I, J, the base and the power. J should be nonnegative. // // Output, int I4_POWER, the value of I^J. // { int k; int value; if ( j < 0 ) { if ( i == 1 ) { value = 1; } else if ( i == 0 ) { cout << "\n"; cout << "I4_POWER - Fatal error!\n"; cout << " I^J requested, with I = 0 and J negative.\n"; exit ( 1 ); } else { value = 0; } } else if ( j == 0 ) { if ( i == 0 ) { cout << "\n"; cout << "I4_POWER - Fatal error!\n"; cout << " I^J requested, with I = 0 and J = 0.\n"; exit ( 1 ); } else { value = 1; } } else if ( j == 1 ) { value = i; } else { value = 1; for ( k = 1; k <= j; k++ ) { value = value * i; } } return value; } //****************************************************************************80 void i4vec_print ( int n, int a[], string title ) //****************************************************************************80 // // Purpose: // // I4VEC_PRINT prints an I4VEC. // // Discussion: // // An I4VEC is a vector of I4's. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 14 November 2003 // // Author: // // John Burkardt // // Parameters: // // Input, int N, the number of components of the vector. // // Input, int A[N], the vector to be printed. // // Input, string TITLE, a title. // { int i; cout << "\n"; cout << title << "\n"; cout << "\n"; for ( i = 0; i < n; i++ ) { cout << " " << setw(8) << i << ": " << setw(8) << a[i] << "\n"; } return; } //****************************************************************************80 double lattice ( int dim_num, int m, int z[], double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // LATTICE applies a lattice integration rule. // // Discussion: // // Because this is a standard lattice rule, it is really only suited // for functions which are periodic, of period 1, in both X and Y. // // For a suitable F, and a given value of M (the number of lattice points), // the performance of the routine is affected by the choice of the // generator vector Z. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 20 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int M, the order (number of points) of the rule. // // Input, int Z[DIM_NUM], the generator vector. Typically, the elements // of Z satisfy 1 <= Z[*] < M, and are relatively prime to M. // This is easy to guarantee if M is itself a prime number. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double LATTICE, the estimated integral. // { int i; int j; double quad; double *x; x = new double[dim_num]; quad = 0.0; for ( j = 0; j <= m - 1; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i] = fmod ( ( double ) ( j * z[i] ) / ( double ) ( m ), 1.0 ); } quad = quad + f ( dim_num, x ); } quad = quad / ( double ) ( m ); delete [] x; return quad; } //****************************************************************************80 double lattice_np0 ( int dim_num, int m, int z[], double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // LATTICE_NP0 applies a lattice integration rule to a nonperiodic function. // // Discussion: // // This routine attempts to modify a lattice rule, suitable for use // with a periodic function, for use with a nonperiodic function F(X), // essentially by applying the lattice rule to the function // // G(X) = ( F(X) + F(1-X) ) / 2 // // This is the rule in 1 dimension. In two dimensions, we have // // G(X,Y) = ( F(X,Y) + F(X,1-Y) + F(1-X,Y) + F(1-X,1-Y) ) / 4 // // with the obvious generalizations to higher dimensions. // // Drawbacks of this approach include: // // * in dimension DIM_NUM, we must evaluate the function F at // 2**DIM_NUM points for every single evaluation of G; // // * the function G, regarded as a periodic function, is continuous, // but not generally differentiable, at 0 and 1. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 November 2008 // // Author: // // John Burkardt // // Reference: // // Seymour Haber, // Parameters for Integrating Periodic Functions of Several Variables, // Mathematics of Computation, // Volume 41, Number 163, July 1983, pages 115-129. // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int M, the order (number of points) of the rule. // // Input, int Z[DIM_NUM], the generator vector. Typically, the elements // of Z satisfy 1 <= Z[*] < M, and are relatively prime to M. // This is easy to guarantee if M is itself a prime number. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double LATTICE_NP0, the estimated integral. // { int dim; int gray_bit; int j; double quad; double *x; double *y; x = new double[dim_num]; y = new double[dim_num]; quad = 0.0; for ( j = 0; j <= m - 1; j++ ) { for ( dim = 0; dim < dim_num; dim++ ) { x[dim] = fmod ( ( double ) ( j * z[dim] ) / ( double ) ( m ), 1.0 ); } // // Generate all DIM_NUM-tuples for which the I-th element is X(I) or 1-X(I). // gray_bit = - dim_num; for ( ; ; ) { gray_next ( dim_num, &gray_bit ); if ( gray_bit == - dim_num ) { break; } if ( gray_bit == 0 ) { for ( dim = 0; dim < dim_num; dim++ ) { y[dim] = x[dim]; } } else { dim = abs ( gray_bit ) - 1; y[dim] = 1.0 - y[dim]; } quad = quad + f ( dim_num, y ); } } quad = quad / ( double ) ( i4_power ( 2, dim_num ) * m ); delete [] x; delete [] y; return quad; } //****************************************************************************80 double lattice_np1 ( int dim_num, int m, int z[], double f ( int dim_num, double x[] ) ) //****************************************************************************80 // // Purpose: // // LATTICE_NP1 applies a lattice integration rule to a nonperiodic function. // // Discussion: // // This routine applies the transformation function // // PHI(T) = 3*T^2 - 2*T^3 // // to a nonperiodic integrand to make it suitable for treatment // by a lattice rule. // // For a suitable F, and a given value of M (the number of lattice points), // the performance of the routine is affected by the choice of the // generator vector Z. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int M, the order (number of points) of the rule. // // Input, int Z[DIM_NUM], the generator vector. Typically, the elements // of Z satisfy 1 <= Z(1:DIM_NUM) < M, and are relatively prime to M. // This is easy to guarantee if M is itself a prime number. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Output, double LATTICE_NP1, the estimated integral. // { int dim; double dphi; int j; double quad; double *x; double *y; x = new double[dim_num]; y = new double[dim_num]; quad = 0.0; for ( j = 0; j <= m - 1; j++ ) { for ( dim = 0; dim < dim_num; dim++ ) { x[dim] = fmod ( ( double ) ( j * z[dim] ) / ( double ) ( m ), 1.0 ); } dphi = 1.0; for ( dim = 0; dim < dim_num; dim++ ) { y[dim] = ( 3.0 - 2.0 * x[dim] ) * x[dim] * x[dim]; dphi = dphi * 6.0 * ( 1.0 - x[dim] ) * x[dim]; } quad = quad + f ( dim_num, y ) * dphi; } quad = quad / ( double ) ( m ); delete [] x; delete [] y; return quad; } //****************************************************************************80 void lattice_print ( int dim_num, int m, int z[], string title ) //****************************************************************************80 // // Purpose: // // LATTICE_PRINT prints the points in a lattice rule. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int M, the number of points to use. // // Input, int Z[DIM_NUM], the generator vector. // // Input, string TITLE, a title. // { int dim; int i; int *y; y = new int[dim_num]; cout << "\n"; cout << title << "\n"; cout << "\n"; for ( i = 0; i <= m - 1; i++ ) { cout << setw(4) << i + 1 << " "; for ( dim = 0; dim < dim_num; dim++ ) { y[dim] = ( i * z[dim] ) % m; cout << setw(4) << y[dim]; } cout << "\n"; } delete [] y; return; } //****************************************************************************80 double monte_carlo ( int dim_num, int m, double f ( int dim_num, double x[] ), int *seed ) //****************************************************************************80 // // Purpose: // // MONTE_CARLO applies a Monte Carlo integration rule. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 November 2008 // // Author: // // John Burkardt // // Reference: // // Ian Sloan, Stephen Joe, // Lattice Methods for Multiple Integration, // Oxford, 1994, // ISBN: 0198534728, // LC: QA311.S56 // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int M, the number of points to use. // // Input, double F ( int DIM_NUM, double X[] ), the name of the // user-supplied routine which evaluates the function. // // Input/output, int *SEED, a seed for the random number generator. // // Output, double MONTE_CARLO, the estimated integral. // { int j; double quad; double *x; quad = 0.0; for ( j = 1; j <= m; j++ ) { x = r8vec_uniform_01 ( dim_num, seed ); quad = quad + f ( dim_num, x ); delete [] x; } quad = quad / ( double ) m; return quad; } //****************************************************************************80 int prime ( int n ) //****************************************************************************80 // // Purpose: // // PRIME returns any of the first PRIME_MAX prime numbers. // // Discussion: // // PRIME_MAX is 1600, and the largest prime stored is 13499. // // Thanks to Bart Vandewoestyne for pointing out a typo, 18 February 2005. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 18 February 2005 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Daniel Zwillinger, // CRC Standard Mathematical Tables and Formulae, // 30th Edition, // CRC Press, 1996, pages 95-98. // // Parameters: // // Input, int N, the index of the desired prime number. // In general, is should be true that 0 <= N <= PRIME_MAX. // N = -1 returns PRIME_MAX, the index of the largest prime available. // N = 0 is legal, returning PRIME = 1. // // Output, int PRIME, the N-th prime. If N is out of range, PRIME // is returned as -1. // { # define PRIME_MAX 1600 int npvec[PRIME_MAX] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007, 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099, 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177, 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271, 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343, 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459, 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567, 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657, 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739, 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859, 10861,10867,10883,10889,10891,10903,10909,10937,10939,10949, 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059, 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149, 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251, 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329, 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443, 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527, 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657, 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777, 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833, 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933, 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011, 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109, 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211, 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289, 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401, 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487, 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553, 12569,12577,12583,12589,12601,12611,12613,12619,12637,12641, 12647,12653,12659,12671,12689,12697,12703,12713,12721,12739, 12743,12757,12763,12781,12791,12799,12809,12821,12823,12829, 12841,12853,12889,12893,12899,12907,12911,12917,12919,12923, 12941,12953,12959,12967,12973,12979,12983,13001,13003,13007, 13009,13033,13037,13043,13049,13063,13093,13099,13103,13109, 13121,13127,13147,13151,13159,13163,13171,13177,13183,13187, 13217,13219,13229,13241,13249,13259,13267,13291,13297,13309, 13313,13327,13331,13337,13339,13367,13381,13397,13399,13411, 13417,13421,13441,13451,13457,13463,13469,13477,13487,13499 }; if ( n == -1 ) { return PRIME_MAX; } else if ( n == 0 ) { return 1; } else if ( n <= PRIME_MAX ) { return npvec[n-1]; } else { cout << "\n"; cout << "PRIME - Fatal error//\n"; cout << " Unexpected input value of n = " << n << "\n"; exit ( 1 ); } return 0; # undef PRIME_MAX } //****************************************************************************80 void r8mat_transpose_print ( int m, int n, double a[], string title ) //****************************************************************************80 // // Purpose: // // R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. // // Discussion: // // An R8MAT is a doubly dimensioned array of R8 values, stored as a vector // in column-major order. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 10 September 2009 // // Author: // // John Burkardt // // Parameters: // // Input, int M, N, the number of rows and columns. // // Input, double A[M*N], an M by N matrix to be printed. // // Input, string TITLE, a title. // { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } //****************************************************************************80 void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, string title ) //****************************************************************************80 // // Purpose: // // R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. // // Discussion: // // An R8MAT is a doubly dimensioned array of R8 values, stored as a vector // in column-major order. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 20 August 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int M, N, the number of rows and columns. // // Input, double A[M*N], an M by N matrix to be printed. // // Input, int ILO, JLO, the first row and column to print. // // Input, int IHI, JHI, the last row and column to print. // // Input, string TITLE, a title. // { # define INCX 5 int i; int i2; int i2hi; int i2lo; int inc; int j; int j2hi; int j2lo; cout << "\n"; cout << title << "\n"; for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); inc = i2hi + 1 - i2lo; cout << "\n"; cout << " Row: "; for ( i = i2lo; i <= i2hi; i++ ) { cout << setw(7) << i - 1 << " "; } cout << "\n"; cout << " Col\n"; cout << "\n"; j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { cout << setw(5) << j - 1 << ":"; for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; cout << setw(14) << a[(i-1)+(j-1)*m]; } cout << "\n"; } } return; # undef INCX } //****************************************************************************80 double *r8vec_uniform_01 ( int n, int *seed ) //****************************************************************************80 // // Purpose: // // R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. // // Discussion: // // An R8VEC is a vector of R8's. // // This routine implements the recursion // // seed = 16807 * seed mod ( 2**31 - 1 ) // unif = seed / ( 2**31 - 1 ) // // The integer arithmetic never requires more than 32 bits, // including a sign bit. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 August 2004 // // Author: // // John Burkardt // // Reference: // // Paul Bratley, Bennett Fox, Linus Schrage, // A Guide to Simulation, // Springer Verlag, pages 201-202, 1983. // // Bennett Fox, // Algorithm 647: // Implementation and Relative Efficiency of Quasirandom // Sequence Generators, // ACM Transactions on Mathematical Software, // Volume 12, Number 4, pages 362-376, 1986. // // Parameters: // // Input, int N, the number of entries in the vector. // // Input/output, int *SEED, a seed for the random number generator. // // Output, double R8VEC_UNIFORM_01[N], the vector of pseudorandom values. // { int i; int k; double *r; r = new double[n]; for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } //****************************************************************************80 void timestamp ( ) //****************************************************************************80 // // Purpose: // // TIMESTAMP prints the current YMDHMS date as a time stamp. // // Example: // // May 31 2001 09:45:54 AM // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 03 October 2003 // // Author: // // John Burkardt // // Parameters: // // None // { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); cout << time_buffer << "\n"; return; # undef TIME_SIZE } //****************************************************************************80 void tuple_next ( int m1, int m2, int n, int *rank, int x[] ) //****************************************************************************80 // // Purpose: // // TUPLE_NEXT computes the next element of a tuple space. // // Discussion: // // The elements are N vectors. Each entry is constrained to lie // between M1 and M2. The elements are produced one at a time. // The first element is // (M1,M1,...,M1), // the second element is // (M1,M1,...,M1+1), // and the last element is // (M2,M2,...,M2) // Intermediate elements are produced in lexicographic order. // // Example: // // N = 2, M1 = 1, M2 = 3 // // INPUT OUTPUT // ------- ------- // Rank X Rank X // ---- --- ----- --- // 0 * * 1 1 1 // 1 1 1 2 1 2 // 2 1 2 3 1 3 // 3 1 3 4 2 1 // 4 2 1 5 2 2 // 5 2 2 6 2 3 // 6 2 3 7 3 1 // 7 3 1 8 3 2 // 8 3 2 9 3 3 // 9 3 3 0 0 0 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 April 2003 // // Author: // // John Burkardt // // Parameters: // // Input, int M1, M2, the minimum and maximum entries. // // Input, int N, the number of components. // // Input/output, int *RANK, counts the elements. // On first call, set RANK to 0. Thereafter, the output value of RANK // will indicate the order of the element returned. When there are no // more elements, RANK will be returned as 0. // // Input/output, int X[N], on input the previous tuple. // On output, the next tuple. // { int i; int j; if ( m2 < m1 ) { *rank = 0; return; } if ( *rank <= 0 ) { for ( i = 0; i < n; i++ ) { x[i] = m1; } *rank = 1; } else { *rank = *rank + 1; i = n - 1; for ( ; ; ) { if ( x[i] < m2 ) { x[i] = x[i] + 1; break; } x[i] = m1; if ( i == 0 ) { *rank = 0; for ( j = 0; j < n; j++ ) { x[j] = m1; } break; } i = i - 1; } } return; }