1/*
2*/
3#if defined(USE_POW)
4#define r23 pow(0.5, 23.0)
5#define r46 (r23*r23)
6#define t23 pow(2.0, 23.0)
7#define t46 (t23*t23)
8#else
9#define r23 (0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5)
10#define r46 (r23*r23)
11#define t23 (2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0*2.0)
12#define t46 (t23*t23)
13#endif
14
15/*c---------------------------------------------------------------------
16c---------------------------------------------------------------------*/
17
18double randlc (double *x, double a) {
19
20/*c---------------------------------------------------------------------
21c---------------------------------------------------------------------*/
22
23/*c---------------------------------------------------------------------
24c
25c   This routine returns a uniform pseudorandom double precision number in the
26c   range (0, 1) by using the linear congruential generator
27c
28c   x_{k+1} = a x_k  (mod 2^46)
29c
30c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
31c   before repeating.  The argument A is the same as 'a' in the above formula,
32c   and X is the same as x_0.  A and X must be odd double precision integers
33c   in the range (1, 2^46).  The returned value RANDLC is normalized to be
34c   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
35c   the new seed x_1, so that subsequent calls to RANDLC using the same
36c   arguments will generate a continuous sequence.
37c
38c   This routine should produce the same results on any computer with at least
39c   48 mantissa bits in double precision floating point data.  On 64 bit
40c   systems, double precision should be disabled.
41c
42c   David H. Bailey     October 26, 1990
43c
44c---------------------------------------------------------------------*/
45
46    double t1,t2,t3,t4,a1,a2,x1,x2,z;
47
48/*c---------------------------------------------------------------------
49c   Break A into two parts such that A = 2^23 * A1 + A2.
50c---------------------------------------------------------------------*/
51    t1 = r23 * a;
52    a1 = (int)t1;
53    a2 = a - t23 * a1;
54
55/*c---------------------------------------------------------------------
56c   Break X into two parts such that X = 2^23 * X1 + X2, compute
57c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
58c   X = 2^23 * Z + A2 * X2  (mod 2^46).
59c---------------------------------------------------------------------*/
60    t1 = r23 * (*x);
61    x1 = (int)t1;
62    x2 = (*x) - t23 * x1;
63    t1 = a1 * x2 + a2 * x1;
64    t2 = (int)(r23 * t1);
65    z = t1 - t23 * t2;
66    t3 = t23 * z + a2 * x2;
67    t4 = (int)(r46 * t3);
68    (*x) = t3 - t46 * t4;
69
70    return (r46 * (*x));
71}
72
73/*c---------------------------------------------------------------------
74c---------------------------------------------------------------------*/
75
76void vranlc (int n, double *x_seed, double a, double y[]) {
77
78/*c---------------------------------------------------------------------
79c---------------------------------------------------------------------*/
80
81/*c---------------------------------------------------------------------
82c
83c   This routine generates N uniform pseudorandom double precision numbers in
84c   the range (0, 1) by using the linear congruential generator
85c
86c   x_{k+1} = a x_k  (mod 2^46)
87c
88c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
89c   before repeating.  The argument A is the same as 'a' in the above formula,
90c   and X is the same as x_0.  A and X must be odd double precision integers
91c   in the range (1, 2^46).  The N results are placed in Y and are normalized
92c   to be between 0 and 1.  X is updated to contain the new seed, so that
93c   subsequent calls to VRANLC using the same arguments will generate a
94c   continuous sequence.  If N is zero, only initialization is performed, and
95c   the variables X, A and Y are ignored.
96c
97c   This routine is the standard version designed for scalar or RISC systems.
98c   However, it should produce the same results on any single processor
99c   computer with at least 48 mantissa bits in double precision floating point
100c   data.  On 64 bit systems, double precision should be disabled.
101c
102c---------------------------------------------------------------------*/
103
104    int i;
105    double x,t1,t2,t3,t4,a1,a2,x1,x2,z;
106
107/*c---------------------------------------------------------------------
108c   Break A into two parts such that A = 2^23 * A1 + A2.
109c---------------------------------------------------------------------*/
110    t1 = r23 * a;
111    a1 = (int)t1;
112    a2 = a - t23 * a1;
113    x = *x_seed;
114
115/*c---------------------------------------------------------------------
116c   Generate N results.   This loop is not vectorizable.
117c---------------------------------------------------------------------*/
118    for (i = 1; i <= n; i++) {
119
120/*c---------------------------------------------------------------------
121c   Break X into two parts such that X = 2^23 * X1 + X2, compute
122c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
123c   X = 2^23 * Z + A2 * X2  (mod 2^46).
124c---------------------------------------------------------------------*/
125        t1 = r23 * x;
126        x1 = (int)t1;
127        x2 = x - t23 * x1;
128        t1 = a1 * x2 + a2 * x1;
129        t2 = (int)(r23 * t1);
130        z = t1 - t23 * t2;
131        t3 = t23 * z + a2 * x2;
132        t4 = (int)(r46 * t3);
133        x = t3 - t46 * t4;
134        y[i] = r46 * x;
135    }
136    *x_seed = x;
137}
138