# include # include # include # include # include # include # include using namespace std; # include "latin_center.hpp" //****************************************************************************80 int get_seed ( ) //****************************************************************************80 // // Purpose: // // GET_SEED returns a random seed for the random number generator. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 15 September 2003 // // Author: // // John Burkardt // // Parameters: // // Output, int GET_SEED, a random seed value. // { # define I_MAX 2147483647 time_t clock; 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. // // Licensing: // // This code is distributed under the MIT license. // // 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 integers. // // Licensing: // // This code is distributed under the MIT license. // // 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 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. // // Licensing: // // This code is distributed under the MIT license. // // 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_center ( int dim_num, int point_num, int *seed, double x[] ) //****************************************************************************80 // // Purpose: // // LATIN_CENTER returns center points in a Latin square. // // Discussion: // // In each spatial dimension, there will be exactly one // point with the coordinate value // // ( 1, 3, 5, ..., 2*point_num-1 ) / ( 2 * point_num ) // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 12 March 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 LATIN_CENTER[DIM_NUM*POINT_NUM], the points. // { int base = 0; 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. // k = 0; for ( i = 0; i < dim_num; i++ ) { perm = perm_uniform ( point_num, base, seed ); for ( j = 0; j < point_num; j++ ) { r = 0.5; x[k] = ( ( ( double ) perm[j] ) + r ) / ( ( double ) point_num ); k = k + 1; } delete [] perm; } return; } //****************************************************************************80 int *perm_uniform ( int n, int base, int *seed ) //****************************************************************************80 // // Purpose: // // PERM_UNIFORM selects a random permutation of N objects. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 31 October 2008 // // Author: // // John Burkardt // // Reference: // // Albert Nijenhuis, Herbert Wilf, // Combinatorial Algorithms, // Academic Press, 1978, second edition, // ISBN 0-12-519260-6. // // Parameters: // // Input, int N, the number of objects to be permuted. // // Input, int BASE, is 0 for a 0-based permutation and 1 for // a 1-based permutation. // // Input/output, int *SEED, a seed for the random number generator. // // Output, int PERM_UNIFORM[N], a permutation of (BASE, BASE+1, ..., BASE+N-1). // { int i; int j; int k; int *p; p = new int[n]; for ( i = 0; i < n; i++ ) { p[i] = i + base; } for ( i = 0; i < n; i++ ) { j = i4_uniform ( i, n - 1, seed ); k = p[i]; p[i] = p[j]; p[j] = k; } return p; } //****************************************************************************80 float r4_abs ( float x ) //****************************************************************************80 // // Purpose: // // R4_ABS returns the absolute value of an R4. // // Licensing: // // This code is distributed under the MIT license. // // 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. // // Example: // // 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 // // Licensing: // // This code is distributed under the MIT license. // // 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 void r8mat_write ( string output_filename, int m, int n, double table[] ) //****************************************************************************80 // // Purpose: // // R8MAT_WRITE writes an R8MAT file. // // Discussion: // // An R8MAT is an array of R8's. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 29 June 2009 // // Author: // // John Burkardt // // Parameters: // // Input, string OUTPUT_FILENAME, the output filename. // // Input, int M, the spatial dimension. // // Input, int N, the number of points. // // Input, double TABLE[M*N], the table data. // { int i; int j; ofstream output; // // Open the file. // output.open ( output_filename.c_str ( ) ); if ( !output ) { cerr << "\n"; cerr << "R8MAT_WRITE - Fatal error!\n"; cerr << " Could not open the output file.\n"; return; } // // Write the data. // for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { output << " " << setw(24) << setprecision(16) << table[i+j*m]; } output << "\n"; } // // Close the file. // output.close ( ); return; } //****************************************************************************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 MIT 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 }