1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
4 double precision function randlc (x, a)
6 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
9 c---------------------------------------------------------------------
11 c This routine returns a uniform pseudorandom double precision number in the
12 c range (0, 1) by using the linear congruential generator
14 c x_{k+1} = a x_k (mod 2^46)
16 c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
17 c before repeating. The argument A is the same as 'a' in the above formula,
18 c and X is the same as x_0. A and X must be odd double precision integers
19 c in the range (1, 2^46). The returned value RANDLC is normalized to be
20 c between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain
21 c the new seed x_1, so that subsequent calls to RANDLC using the same
22 c arguments will generate a continuous sequence.
24 c This routine should produce the same results on any computer with at least
25 c 48 mantissa bits in double precision floating point data. On 64 bit
26 c systems, double precision should be disabled.
28 c David H. Bailey October 26, 1990
30 c---------------------------------------------------------------------
34 double precision r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
35 parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
38 c---------------------------------------------------------------------
39 c Break A into two parts such that A = 2^23 * A1 + A2.
40 c---------------------------------------------------------------------
45 c---------------------------------------------------------------------
46 c Break X into two parts such that X = 2^23 * X1 + X2, compute
47 c Z = A1 * X2 + A2 * X1 (mod 2^23), and then
48 c X = 2^23 * Z + A2 * X2 (mod 2^46).
49 c---------------------------------------------------------------------
53 t1 = a1 * x2 + a2 * x1
56 t3 = t23 * z + a2 * x2
67 c---------------------------------------------------------------------
68 c---------------------------------------------------------------------
70 subroutine vranlc (n, x, a, y)
72 c---------------------------------------------------------------------
73 c---------------------------------------------------------------------
75 c---------------------------------------------------------------------
77 c This routine generates N uniform pseudorandom double precision numbers in
78 c the range (0, 1) by using the linear congruential generator
80 c x_{k+1} = a x_k (mod 2^46)
82 c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
83 c before repeating. The argument A is the same as 'a' in the above formula,
84 c and X is the same as x_0. A and X must be odd double precision integers
85 c in the range (1, 2^46). The N results are placed in Y and are normalized
86 c to be between 0 and 1. X is updated to contain the new seed, so that
87 c subsequent calls to VRANLC using the same arguments will generate a
88 c continuous sequence. If N is zero, only initialization is performed, and
89 c the variables X, A and Y are ignored.
91 c This routine is the standard version designed for scalar or RISC systems.
92 c However, it should produce the same results on any single processor
93 c computer with at least 48 mantissa bits in double precision floating point
94 c data. On 64 bit systems, double precision should be disabled.
96 c---------------------------------------------------------------------
101 double precision y,r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
103 parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
107 c---------------------------------------------------------------------
108 c Break A into two parts such that A = 2^23 * A1 + A2.
109 c---------------------------------------------------------------------
114 c---------------------------------------------------------------------
115 c Generate N results. This loop is not vectorizable.
116 c---------------------------------------------------------------------
119 c---------------------------------------------------------------------
120 c Break X into two parts such that X = 2^23 * X1 + X2, compute
121 c Z = A1 * X2 + A2 * X1 (mod 2^23), and then
122 c X = 2^23 * Z + A2 * X2 (mod 2^46).
123 c---------------------------------------------------------------------
127 t1 = a1 * x2 + a2 * x1
130 t3 = t23 * z + a2 * x2