3 * FUNCTION RANDLC (X, A)
5 * This routine returns a uniform pseudorandom double precision number in the
6 * range (0, 1) by using the linear congruential generator
8 * x_{k+1} = a x_k (mod 2^46)
10 * where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
11 * before repeating. The argument A is the same as 'a' in the above formula,
12 * and X is the same as x_0. A and X must be odd double precision integers
13 * in the range (1, 2^46). The returned value RANDLC is normalized to be
14 * between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain
15 * the new seed x_1, so that subsequent calls to RANDLC using the same
16 * arguments will generate a continuous sequence.
18 * This routine should produce the same results on any computer with at least
19 * 48 mantissa bits in double precision floating point data. On Cray systems,
20 * double precision should be disabled.
22 * David H. Bailey October 26, 1990
24 * IMPLICIT DOUBLE PRECISION (A-H, O-Z)
25 * SAVE KS, R23, R46, T23, T46
28 * If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
29 * T23 = 2 ^ 23, and T46 = 2 ^ 46. These are computed in loops, rather than
30 * by merely using the ** operator, in order to insure that the results are
31 * exact on all systems. This code assumes that 0.5D0 is represented exactly.
35 /*****************************************************************/
36 /************* R A N D L C ************/
37 /************* ************/
38 /************* portable random number generator ************/
39 /*****************************************************************/
41 double randlc( double *X, double *A )
44 static double R23, R46, T23, T46;
45 double T1, T2, T3, T4;
73 /* Break A into two parts such that A = 2^23 * A1 + A2 and set X = N. */
80 /* Break X into two parts such that X = 2^23 * X1 + X2, compute
81 Z = A1 * X2 + A2 * X1 (mod 2^23), and then
82 X = 2^23 * Z + A2 * X2 (mod 2^46). */
88 T1 = A1 * X2 + A2 * X1;
93 T3 = T23 * Z + A2 * X2;
102 /*****************************************************************/
103 /************ F I N D _ M Y _ S E E D ************/
104 /************ ************/
105 /************ returns parallel random number seq seed ************/
106 /*****************************************************************/