# include # include # include # include # include using namespace std; # include "latin_random.H" //****************************************************************************80 int get_seed ( void ) //****************************************************************************80 // // Purpose: // // GET_SEED returns a random seed for the random number generator. // // Modified: // // 15 September 2003 // // Author: // // John Burkardt // // Parameters: // // Output, int GET_SEED, a random seed value. // { # define I_MAX 2147483647 time_t clock; int i; int ihour; int imin; int isec; int seed; struct tm *lt; time_t tloc; // // If the internal seed is 0, generate a value based on the time. // clock = time ( &tloc ); lt = localtime ( &clock ); // // Hours is 1, 2, ..., 12. // ihour = lt->tm_hour; if ( 12 < ihour ) { ihour = ihour - 12; } // // Move Hours to 0, 1, ..., 11 // ihour = ihour - 1; imin = lt->tm_min; isec = lt->tm_sec; seed = isec + 60 * ( imin + 60 * ihour ); // // We want values in [1,43200], not [0,43199]. // seed = seed + 1; // // Remap SEED from [1,43200] to [1,IMAX]. // seed = ( int ) ( ( ( double ) seed ) * ( ( double ) I_MAX ) / ( 60.0 * 60.0 * 12.0 ) ); // // Never use a seed of 0. // if ( seed == 0 ) { seed = 1; } return seed; # undef I_MAX } //****************************************************************************80 int i4_max ( int i1, int i2 ) //****************************************************************************80 // // Purpose: // // I4_MAX returns the maximum of two I4's. // // Modified: // // 05 May 2003 // // Author: // // John Burkardt // // Parameters: // // Input, int I1, I2, two integers to be compared. // // Output, int I4_MAX, the larger of I1 and I2. // { if ( i2 < i1 ) { return i1; } else { return i2; } } //****************************************************************************80 int i4_min ( int i1, int i2 ) //****************************************************************************80 // // Purpose: // // I4_MIN returns the smaller of two I4's. // // Modified: // // 05 May 2003 // // Author: // // John Burkardt // // Parameters: // // Input, int I1, I2, two integers to be compared. // // Output, int I4_MIN, the smaller of I1 and I2. // { if ( i1 < i2 ) { return i1; } else { return i2; } } //****************************************************************************80 void i4_swap ( int *i, int *j ) //****************************************************************************80 // // Purpose: // // I4_SWAP switches two I4's. // // Modified: // // 07 January 2002 // // Author: // // John Burkardt // // Parameters: // // Input/output, int *I, *J. On output, the values of I and // J have been interchanged. // { int k; k = *i; *i = *j; *j = k; return; } //****************************************************************************80 int i4_uniform ( int a, int b, int *seed ) //****************************************************************************80 // // Purpose: // // I4_UNIFORM returns a scaled pseudorandom I4. // // Discussion: // // The pseudorandom number should be uniformly distributed // between A and B. // // Modified: // // 12 November 2006 // // Author: // // John Burkardt // // Reference: // // Paul Bratley, Bennett Fox, Linus Schrage, // A Guide to Simulation, // Springer Verlag, pages 201-202, 1983. // // Pierre L'Ecuyer, // Random Number Generation, // in Handbook of Simulation, // edited by Jerry Banks, // Wiley Interscience, page 95, 1998. // // 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. // // Peter Lewis, Allen Goodman, James Miller // A Pseudo-Random Number Generator for the System/360, // IBM Systems Journal, // Volume 8, pages 136-143, 1969. // // Parameters: // // Input, int A, B, the limits of the interval. // // Input/output, int *SEED, the "seed" value, which should NOT be 0. // On output, SEED has been updated. // // Output, int I4_UNIFORM, a number between A and B. // { int k; float r; int value; if ( *seed == 0 ) { cerr << "\n"; cerr << "I4_UNIFORM - Fatal error!\n"; cerr << " Input value of SEED = 0.\n"; exit ( 1 ); } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r = ( float ) ( *seed ) * 4.656612875E-10; // // Scale R to lie between A-0.5 and B+0.5. // r = ( 1.0 - r ) * ( ( float ) ( i4_min ( a, b ) ) - 0.5 ) + r * ( ( float ) ( i4_max ( a, b ) ) + 0.5 ); // // Use rounding to convert R to an integer between A and B. // value = r4_nint ( r ); value = i4_max ( value, i4_min ( a, b ) ); value = i4_min ( value, i4_max ( a, b ) ); return value; } //****************************************************************************80 void latin_random ( int dim_num, int point_num, int *seed, double x[] ) //****************************************************************************80 // // Purpose: // // LATIN_RANDOM returns points in a Latin Random square. // // Discussion: // // In each spatial dimension, there will be exactly one // point whose coordinate value lies between consecutive // values in the list: // // ( 0, 1, 2, ..., point_num ) / point_num // // Modified: // // 08 April 2003 // // Author: // // John Burkardt // // Parameters: // // Input, int DIM_NUM, the spatial dimension. // // Input, int POINT_NUM, the number of points. // // Input/output, int *SEED, a seed for UNIFORM. // // Output, double X[DIM_NUM,POINT_NUM], the points. // { int i; int j; int k; int* perm; double r; // // For spatial dimension I, // pick a random permutation of 1 to POINT_NUM, // force the corresponding I-th components of X to lie in the // interval ( PERM[J]-1, PERM[J] ) / POINT_NUM. // perm = new int[point_num]; k = 0; for ( i = 0; i < dim_num; i++ ) { perm_random ( point_num, seed, perm ); for ( j = 0; j < point_num; j++ ) { r = r8_uniform_01 ( seed ); x[k] = ( ( ( double ) ( perm[j] - 1 ) ) + r ) / ( ( double ) point_num ); k = k + 1; } } delete [] perm; return; } //****************************************************************************80 void perm_random ( int n, int *seed, int p[] ) //****************************************************************************80 // // Purpose: // // PERM_RANDOM selects a random permutation of N objects. // // Reference: // // A Nijenhuis and H Wilf, // Combinatorial Algorithms, // Academic Press, 1978, second edition, // ISBN 0-12-519260-6. // // Modified: // // 25 May 1999 // // Parameters: // // Input, int N, the number of objects to be permuted. // // Input/output, int *SEED, a seed for UNIFORM. // // Output, int P[N], a permutation of (1, 2, ..., N). // { int i; int j; for ( i = 0; i < n; i++ ) { p[i] = i + 1; } for ( i = 1; i <= n; i++ ) { j = i4_uniform ( i, n, seed ); i4_swap ( &p[i-1], &p[j-1] ); } return; } //****************************************************************************80 float r4_abs ( float x ) //****************************************************************************80 // // Purpose: // // R4_ABS returns the absolute value of an R4. // // Modified: // // 01 December 2006 // // Author: // // John Burkardt // // Parameters: // // Input, float X, the quantity whose absolute value is desired. // // Output, float R4_ABS, the absolute value of X. // { float value; if ( 0.0 <= x ) { value = x; } else { value = -x; } return value; } //****************************************************************************80 int r4_nint ( float x ) //****************************************************************************80 // // Purpose: // // R4_NINT returns the nearest integer to an R4. // // Examples: // // X R4_NINT // // 1.3 1 // 1.4 1 // 1.5 1 or 2 // 1.6 2 // 0.0 0 // -0.7 -1 // -1.1 -1 // -1.6 -2 // // Modified: // // 14 November 2006 // // Author: // // John Burkardt // // Parameters: // // Input, float X, the value. // // Output, int R4_NINT, the nearest integer to X. // { int value; if ( x < 0.0 ) { value = - ( int ) ( r4_abs ( x ) + 0.5 ); } else { value = ( int ) ( r4_abs ( x ) + 0.5 ); } return value; } //****************************************************************************80 double r8_uniform_01 ( int *seed ) //****************************************************************************80 // // Purpose: // // 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 // // Modified: // // 11 August 2004 // // Author: // // John Burkardt // // Reference: // // Paul Bratley, Bennett Fox, Linus Schrage, // A Guide to Simulation, // Springer Verlag, pages 201-202, 1983. // // Pierre L'Ecuyer, // Random Number Generation, // in Handbook of Simulation // edited by Jerry Banks, // Wiley Interscience, page 95, 1998. // // 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. // // Peter Lewis, Allen Goodman, James Miller, // A Pseudo-Random Number Generator for the System/360, // IBM Systems Journal, // Volume 8, pages 136-143, 1969. // // Parameters: // // Input/output, int *SEED, the "seed" value. Normally, this // value should not be 0. On output, SEED has been updated. // // Output, double R8_UNIFORM_01, a new pseudorandom variate, // strictly between 0 and 1. // { int k; double r; if ( *seed == 0 ) { cerr << "\n"; cerr << "R8_UNIFORM_01 - Fatal error!\n"; cerr << " Input value of SEED = 0.\n"; exit ( 1 ); } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } // // Although SEED can be represented exactly as a 32 bit integer, // it generally cannot be represented exactly as a 32 bit real number! // r = ( double ) ( *seed ) * 4.656612875E-10; return r; } //****************************************************************************80 void timestamp ( void ) //****************************************************************************80 // // Purpose: // // TIMESTAMP prints the current YMDHMS date as a time stamp. // // Example: // // May 31 2001 09:45:54 AM // // Modified: // // 03 October 2003 // // Author: // // John Burkardt // // Parameters: // // None // { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; size_t len; time_t now; now = time ( NULL ); tm = localtime ( &now ); len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); cout << time_buffer << "\n"; return; # undef TIME_SIZE }