# 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 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. Input: double *X, the seed value: X should be an odd integer value in the range 1, 2^46. 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 *X: X has been updated by the linear congruential function. 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. Input: double X, the initial seed (with index 0). 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; }