# include # include # include # include "randlc.h" /******************************************************************************/ double randlc ( double *x ) /******************************************************************************/ /* Purpose: RANDLC returns a uniform pseudorandom double precision number. Discussion: This function, when called repeatedly, computes a sequence of pseudorandom double precision numbers in the range (0,1). The calculation involves a series of integer seeds X(0), X(1), ..., X(K+1), and a fixed multiplier A, so that X(K+1) = A * X(K) mod 2^46 and the pseudorandom value is then computed by RANDLC = X(K+1) / 2^46. This scheme generates 2^44 numbers before repeating. The multiplier A and the seed X must be odd double precision integers in the range (1, 2^46). Licensing: This code is distributed under the MIT license. Modified: 08 March 2010 Author: Original C version by David Bailey. This C version by John Burkardt. Reference: David Bailey, Eric Barszcz, John Barton, D Browning, Robert Carter, Leonardo Dagum, Rod Fatoohi, Samuel Fineberg, Paul Frederickson, Thomas Lasinski, Robert Schreiber, Horst Simon, V Venkatakrishnan, Sisira Weeratunga, The NAS Parallel Benchmarks, RNR Technical Report RNR-94-007, March 1994. Donald Knuth, The Art of Computer Programming, Volume 2, Seminumerical Algorithms, Third Edition, Addison Wesley, 1997, ISBN: 0201896842, LC: QA76.6.K64. Parameters: Input/output, double *X, the seed value. On input, X should be an odd integer value in the range 1, 2^46. On output, X has been updated by the linear congruential function. If the input value is zero, it is replaced by 314159265. If the input value is negative, it is replaced by its absolute value. Output, double RANDLC, the pseudorandom value corresponding to the input value of X. */ { static double a = 1220703125.00; static double a1; static double a2; int i; static int ks = 0; static double r23; static double r46; double t1; double t2; static double t23; double t3; double t4; static double t46; double value; double x1; double x2; double z; /* If this is the first call, compute R23 = 2 ^ -23, R46 = 2 ^ -46, T23 = 2 ^ 23, T46 = 2 ^ 46. These are computed in loops, rather than by merely using the power operator, in order to insure that the results are exact on all systems. */ if ( ks == 0 ) { r23 = 1.0; r46 = 1.0; t23 = 1.0; t46 = 1.0; for ( i = 1; i <= 23; i++ ) { r23 = 0.50 * r23; t23 = 2.0 * t23; } for ( i = 1; i <= 46; i++ ) { r46 = 0.50 * r46; t46 = 2.0 * t46; } /* Break A into two parts such that A = 2^23 * A1 + A2. */ t1 = r23 * a; a1 = ( double ) ( int ) t1; a2 = a - t23 * a1; ks = 1; } /* Deal with a 0 input value of X. */ if ( *x == 0 ) { *x = 314159265.0; } /* Deal somewhat arbitrarily with negative input X. */ if ( *x < 0 ) { *x = - ( *x ); } /* Break X into two parts X1 and X2 such that: X = 2^23 * X1 + X2, then compute Z = A1 * X2 + A2 * X1 (mod 2^23) X = 2^23 * Z + A2 * X2 (mod 2^46). */ t1 = r23 * *x; x1 = ( double ) ( int ) t1; x2 = *x - t23 * x1; t1 = a1 * x2 + a2 * x1; t2 = ( double ) ( int ) ( r23 * t1 ); z = t1 - t23 * t2; t3 = t23 * z + a2 * x2; t4 = ( double ) ( int ) ( r46 * t3 ); *x = t3 - t46 * t4; value = r46 * ( *x ); return value; } /******************************************************************************/ double randlc_jump ( double x, int k ) /******************************************************************************/ /* Purpose: RANDLC_JUMP returns the K-th element of a uniform pseudorandom sequence. Discussion: The sequence uses the linear congruential generator: X(K+1) = A * X(K) mod 2^46 The K-th element, which can be represented as X(K) = A^K * X(0) mod 2^46 is computed directly using the binary algorithm for exponentiation. Licensing: This code is distributed under the MIT license. Modified: 12 March 2010 Author: John Burkardt Reference: David Bailey, Eric Barszcz, John Barton, D Browning, Robert Carter, Leonardo Dagum, Rod Fatoohi, Samuel Fineberg, Paul Frederickson, Thomas Lasinski, Robert Schreiber, Horst Simon, V Venkatakrishnan, Sisira Weeratunga, The NAS Parallel Benchmarks, RNR Technical Report RNR-94-007, March 1994. Donald Knuth, The Art of Computer Programming, Volume 2, Seminumerical Algorithms, Third Edition, Addison Wesley, 1997, ISBN: 0201896842, LC: QA76.6.K64. Parameters: Input, double X, the initial seed (with index 0). Input, int K, the index of the desired value. Output, double RANDLC_JUMP, the K-th value in the sequence. */ { static double a = 1220703125.0; static double a1; static double a2; double b; double b1; double b2; int i; int j; int k2; static int ks = 0; int m; static double r23; static double r46; double t1; double t2; static double t23; double t3; double t4; static double t46; int twom; double value; double x1; double x2; double xk; double z; /* If this is the first call, compute R23 = 2 ^ -23, R46 = 2 ^ -46, T23 = 2 ^ 23, T46 = 2 ^ 46. These are computed in loops, rather than by merely using the power operator, in order to insure that the results are exact on all systems. */ if ( ks == 0 ) { r23 = 1.0; r46 = 1.0; t23 = 1.0; t46 = 1.0; for ( i = 1; i <= 23; i++ ) { r23 = 0.50 * r23; t23 = 2.0 * t23; } for ( i = 1; i <= 46; i++ ) { r46 = 0.50 * r46; t46 = 2.0 * t46; } /* Break A into two parts such that A = 2^23 * A1 + A2. */ t1 = r23 * a; a1 = ( double ) ( int ) t1; a2 = a - t23 * a1; ks = 1; } if ( k < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "RANDLC_JUMP - Fatal error!\n" ); fprintf ( stderr, " K < 0.\n" ); exit ( 1 ); } else if ( k == 0 ) { xk = x; } /* Find M so that K < 2^M. */ else { k2 = k; xk = x; m = 1; twom = 2; while ( twom <= k ) { twom = twom * 2; m = m + 1; } b = a; b1 = a1; b2 = a2; for ( i = 1; i <= m; i++ ) { j = k2 / 2; /* Replace X by A * X, if appropriate. */ if ( 2 * j != k2 ) { t1 = r23 * xk; x1 = ( double ) ( ( int ) ( t1 ) ); x2 = xk - t23 * x1; t1 = b1 * x2 + b2 * x1; t2 = ( double ) ( ( int ) ( r23 * t1 ) ); z = t1 - t23 * t2; t3 = t23 * z + b2 * x2; t4 = ( double ) ( ( int ) ( r46 * t3 ) ); xk = t3 - t46 * t4; } /* Replace A by A * A mod 2^46. */ t1 = r23 * b; x1 = ( double ) ( ( int ) ( t1 ) ); x2 = b - t23 * x1; t1 = b1 * x2 + b2 * x1; t2 = ( double ) ( ( int ) ( r23 * t1 ) ); z = t1 - t23 * t2; t3 = t23 * z + b2 * x2; t4 = ( double ) ( ( int ) ( r46 * t3 ) ); b = t3 - t46 * t4; /* Update A1, A2. */ t1 = r23 * b; b1 = ( double ) ( ( int ) ( t1 ) ); b2 = b - t23 * b1; k2 = j; } } value = r46 * xk; return value; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 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 ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }