2 * FUNCTION RANDLC (X, A)
4 * This routine returns a uniform pseudorandom double precision number in the
5 * range (0, 1) by using the linear congruential generator
7 * x_{k+1} = a x_k (mod 2^46)
9 * where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
10 * before repeating. The argument A is the same as 'a' in the above formula,
11 * and X is the same as x_0. A and X must be odd double precision integers
12 * in the range (1, 2^46). The returned value RANDLC is normalized to be
13 * between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain
14 * the new seed x_1, so that subsequent calls to RANDLC using the same
15 * arguments will generate a continuous sequence.
17 * This routine should produce the same results on any computer with at least
18 * 48 mantissa bits in double precision floating point data. On Cray systems,
19 * double precision should be disabled.
21 * David H. Bailey October 26, 1990
23 * IMPLICIT DOUBLE PRECISION (A-H, O-Z)
24 * SAVE KS, R23, R46, T23, T46
27 * If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
28 * T23 = 2 ^ 23, and T46 = 2 ^ 46. These are computed in loops, rather than
29 * by merely using the ** operator, in order to insure that the results are
30 * exact on all systems. This code assumes that 0.5D0 is represented exactly.
32 double randlc(double *X, double*A)
35 static double R23, R46, T23, T46;
36 double T1, T2, T3, T4;
48 for (i=1; i<=23; i++) {
52 for (i=1; i<=46; i++) {
59 /* Break A into two parts such that A = 2^23 * A1 + A2 and set X = N. */
65 /* Break X into two parts such that X = 2^23 * X1 + X2, compute
66 Z = A1 * X2 + A2 * X1 (mod 2^23), and then X = 2^23 * Z + A2 * X2 (mod 2^46). */
71 T1 = A1 * X2 + A2 * X1;
76 T3 = T23 * Z + A2 * X2;