1 | /* RAND functions using 64b INTEGERs |
---|
2 | |
---|
3 | F. CANTONNET - HPCL - GWU */ |
---|
4 | |
---|
5 | double randlc (double *x, double a) |
---|
6 | { |
---|
7 | /* This routine returns a uniform pseudorandom double precision number in the |
---|
8 | range (0, 1) by using the linear congruential generator |
---|
9 | |
---|
10 | x_{k+1} = a x_k (mod 2^46) |
---|
11 | |
---|
12 | where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers |
---|
13 | before repeating. The argument A is the same as 'a' in the above formula, |
---|
14 | and X is the same as x_0. A and X must be odd double precision integers |
---|
15 | in the range (1, 2^46). The returned value RANDLC is normalized to be |
---|
16 | between 0 and 1, i.e. RANDLC = 2^(-46) * x_1. X is updated to contain |
---|
17 | the new seed x_1, so that subsequent calls to RANDLC using the same |
---|
18 | arguments will generate a continuous sequence. |
---|
19 | |
---|
20 | This routine should produce the same results on any computer with at least |
---|
21 | 48 mantissa bits in double precision floating point data. On 64 bit |
---|
22 | systems, double precision should be disabled. |
---|
23 | |
---|
24 | David H. Bailey October 26, 1990 */ |
---|
25 | |
---|
26 | unsigned long long i246m1, Lx, La; |
---|
27 | double d2m46; |
---|
28 | |
---|
29 | d2m46 = 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
30 | 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
31 | 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
32 | 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
33 | 0.5*0.5*0.5*0.5*0.5*0.5; |
---|
34 | //d2m46 = pow( 0.5, 46 ); |
---|
35 | |
---|
36 | i246m1 = 0x00003FFFFFFFFFFFLL; |
---|
37 | |
---|
38 | Lx = *x; |
---|
39 | La = a; |
---|
40 | |
---|
41 | Lx = (Lx*La)&i246m1; |
---|
42 | *x = (double) Lx; |
---|
43 | return (d2m46 * (*x)); |
---|
44 | } |
---|
45 | |
---|
46 | void vranlc (int n, double *x_seed, double a, double y[]) { |
---|
47 | /* This routine generates N uniform pseudorandom double precision numbers in |
---|
48 | the range (0, 1) by using the linear congruential generator |
---|
49 | |
---|
50 | x_{k+1} = a x_k (mod 2^46) |
---|
51 | |
---|
52 | where 0 < x_k < 2^46 and 0 < a < 2^46. This scheme generates 2^44 numbers |
---|
53 | before repeating. The argument A is the same as 'a' in the above formula, |
---|
54 | and X is the same as x_0. A and X must be odd double precision integers |
---|
55 | in the range (1, 2^46). The N results are placed in Y and are normalized |
---|
56 | to be between 0 and 1. X is updated to contain the new seed, so that |
---|
57 | subsequent calls to VRANLC using the same arguments will generate a |
---|
58 | continuous sequence. If N is zero, only initialization is performed, and |
---|
59 | the variables X, A and Y are ignored. |
---|
60 | |
---|
61 | This routine is the standard version designed for scalar or RISC systems. |
---|
62 | However, it should produce the same results on any single processor |
---|
63 | computer with at least 48 mantissa bits in double precision floating point |
---|
64 | data. On 64 bit systems, double precision should be disabled. |
---|
65 | */ |
---|
66 | int i; |
---|
67 | double x; |
---|
68 | unsigned long long i246m1, Lx, La; |
---|
69 | double d2m46; |
---|
70 | |
---|
71 | d2m46 = 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
72 | 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
73 | 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
74 | 0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5*0.5* |
---|
75 | 0.5*0.5*0.5*0.5*0.5*0.5; |
---|
76 | // d2m46 = pow( 0.5, 46.0 ); |
---|
77 | i246m1 = 0x00003FFFFFFFFFFFLL; |
---|
78 | |
---|
79 | x = *x_seed; |
---|
80 | Lx = x; |
---|
81 | La = a; |
---|
82 | |
---|
83 | for (i = 1; i <= n; i++) |
---|
84 | { |
---|
85 | Lx = ((Lx*La)&i246m1); |
---|
86 | x = (double) Lx; |
---|
87 | y[i] = d2m46 * x; |
---|
88 | } |
---|
89 | *x_seed = x; |
---|
90 | } |
---|