1 c---------------------------------------------------------------------
2 double precision function randlc (x, a)
3 c---------------------------------------------------------------------
5 c---------------------------------------------------------------------
7 c This routine returns a uniform pseudorandom double precision number in the
8 c range (0, 1) by using the linear congruential generator
10 c x_{k+1} = a x_k (mod 2^46)
12 c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
13 c before repeating. The argument A is the same as 'a' in the above formula,
14 c and X is the same as x_0. A and X must be odd double precision integers
15 c in the range (1, 2^46). The returned value RANDLC is normalized to be
16 c between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain
17 c the new seed x_1, so that subsequent calls to RANDLC using the same
18 c arguments will generate a continuous sequence.
20 c This routine should produce the same results on any computer with at least
21 c 48 mantissa bits in double precision floating point data. On 64 bit
22 c systems, double precision should be disabled.
24 c David H. Bailey October 26, 1990
26 c---------------------------------------------------------------------
30 double precision r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
31 parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
34 c---------------------------------------------------------------------
35 c Break A into two parts such that A = 2^23 * A1 + A2.
36 c---------------------------------------------------------------------
41 c---------------------------------------------------------------------
42 c Break X into two parts such that X = 2^23 * X1 + X2, compute
43 c Z = A1 * X2 + A2 * X1 (mod 2^23), and then
44 c X = 2^23 * Z + A2 * X2 (mod 2^46).
45 c---------------------------------------------------------------------
51 t1 = a1 * x2 + a2 * x1
54 t3 = t23 * z + a2 * x2
63 c---------------------------------------------------------------------
64 c---------------------------------------------------------------------
66 subroutine vranlc (n, x, a, y)
68 c---------------------------------------------------------------------
69 c---------------------------------------------------------------------
71 c---------------------------------------------------------------------
72 c This routine generates N uniform pseudorandom double precision numbers in
73 c the range (0, 1) by using the linear congruential generator
75 c x_{k+1} = a x_k (mod 2^46)
77 c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
78 c before repeating. The argument A is the same as 'a' in the above formula,
79 c and X is the same as x_0. A and X must be odd double precision integers
80 c in the range (1, 2^46). The N results are placed in Y and are normalized
81 c to be between 0 and 1. X is updated to contain the new seed, so that
82 c subsequent calls to RANDLC using the same arguments will generate a
83 c continuous sequence.
85 c This routine generates the output sequence in batches of length NV, for
86 c convenience on vector computers. This routine should produce the same
87 c results on any computer with at least 48 mantissa bits in double precision
88 c floating point data. On Cray systems, double precision should be disabled.
90 c David H. Bailey August 30, 1990
91 c---------------------------------------------------------------------
94 double precision x, a, y(*)
96 double precision r23, r46, t23, t46
98 parameter (r23 = 2.d0 ** (-23), r46 = r23 * r23, t23 = 2.d0 ** 23,
99 > t46 = t23 * t23, nv = 64)
100 double precision xv(nv), t1, t2, t3, t4, an, a1, a2, x1, x2, yy
103 double precision randlc
105 c---------------------------------------------------------------------
106 c Compute the first NV elements of the sequence using RANDLC.
107 c---------------------------------------------------------------------
112 xv(i) = t46 * randlc (t1, a)
115 c---------------------------------------------------------------------
116 c It is not necessary to compute AN, A1 or A2 unless N is greater than NV.
117 c---------------------------------------------------------------------
120 c---------------------------------------------------------------------
121 c Compute AN = AA ^ NV (mod 2^46) using successive calls to RANDLC.
122 c---------------------------------------------------------------------
132 c---------------------------------------------------------------------
133 c Break AN into two parts such that AN = 2^23 * A1 + A2.
134 c---------------------------------------------------------------------
140 c---------------------------------------------------------------------
141 c Compute N pseudorandom results in batches of size NV.
142 c---------------------------------------------------------------------
146 c---------------------------------------------------------------------
147 c Compute up to NV results based on the current seed vector XV.
148 c---------------------------------------------------------------------
153 c---------------------------------------------------------------------
154 c If this is the last pass through the 140 loop, it is not necessary to
155 c update the XV vector.
156 c---------------------------------------------------------------------
157 if (j + n1 .eq. n) goto 150
159 c---------------------------------------------------------------------
160 c Update the XV vector by multiplying each element by AN (mod 2^46).
161 c---------------------------------------------------------------------
165 x2 = xv(i) - t23 * x1
166 t1 = a1 * x2 + a2 * x1
169 t3 = t23 * yy + a2 * x2
171 xv(i) = t3 - t46 * t4
176 c---------------------------------------------------------------------
177 c Save the last seed in X so that subsequent calls to VRANLC will generate
178 c a continuous sequence.
179 c---------------------------------------------------------------------
185 c----- end of program ------------------------------------------------