1 double precision function randlc(x, a)
3 c---------------------------------------------------------------------
5 c This routine returns a uniform pseudorandom double precision number in the
6 c range (0, 1) by using the linear congruential generator
8 c x_{k+1} = a x_k (mod 2^46)
10 c where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers
11 c before repeating. The argument A is the same as 'a' in the above formula,
12 c and X is the same as x_0. A and X must be odd double precision integers
13 c in the range (1, 2^46). The returned value RANDLC is normalized to be
14 c between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain
15 c the new seed x_1, so that subsequent calls to RANDLC using the same
16 c arguments will generate a continuous sequence.
20 integer*8 Lx, La, a1, a2, x1, x2, xa
21 double precision d2m46
22 parameter(d2m46=0.5d0**46)
26 a1 = ibits(La, 23, 23)
28 x1 = ibits(Lx, 23, 23)
30 xa = ishft(ibits(a1*x2+a2*x1, 0, 23), 23) + a2*x2
38 c---------------------------------------------------------------------
39 c---------------------------------------------------------------------
42 SUBROUTINE VRANLC (N, X, A, Y)
45 double precision x, a, y(*)
46 integer*8 Lx, La, a1, a2, x1, x2, xa
47 double precision d2m46
48 parameter(d2m46=0.5d0**46)
52 a1 = ibits(La, 23, 23)
55 x1 = ibits(Lx, 23, 23)
57 xa = ishft(ibits(a1*x2+a2*x1, 0, 23), 23) + a2*x2