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