1 | // file: ranlux.xpp |
---|
2 | #include "ranlux.h" |
---|
3 | #include <stdlib.h> |
---|
4 | #include <stdio.h> |
---|
5 | |
---|
6 | /* This is a lagged fibonacci generator with skipping developed by Luescher. |
---|
7 | The sequence is a series of 24-bit integers, x_n, |
---|
8 | |
---|
9 | x_n = d_n + b_n |
---|
10 | |
---|
11 | where d_n = x_{n-10} - x_{n-24} - c_{n-1}, b_n = 0 if d_n >= 0 and |
---|
12 | b_n = 2^24 if d_n < 0, c_n = 0 if d_n >= 0 and c_n = 1 if d_n < 0, |
---|
13 | where after 24 samples a group of p integers are "skipped", to |
---|
14 | reduce correlations. By default p = 199, but can be increased to |
---|
15 | 365. |
---|
16 | |
---|
17 | The period of the generator is around 10^171. |
---|
18 | |
---|
19 | From: M. Luescher, "A portable high-quality random number generator |
---|
20 | for lattice field theory calculations", Computer Physics |
---|
21 | Communications, 79 (1994) 100-110. |
---|
22 | |
---|
23 | Available on the net as hep-lat/9309020 at http://xxx.lanl.gov/ |
---|
24 | |
---|
25 | See also, |
---|
26 | |
---|
27 | F. James, "RANLUX: A Fortran implementation of the high-quality |
---|
28 | pseudo-random number generator of Luscher", Computer Physics |
---|
29 | Communications, 79 (1994) 111-114 |
---|
30 | |
---|
31 | Kenneth G. Hamilton, F. James, "Acceleration of RANLUX", Computer |
---|
32 | Physics Communications, 101 (1997) 241-248 |
---|
33 | |
---|
34 | Kenneth G. Hamilton, "Assembler RANLUX for PCs", Computer Physics |
---|
35 | Communications, 101 (1997) 249-253 */ |
---|
36 | |
---|
37 | namespace siscone{ |
---|
38 | |
---|
39 | static const unsigned long int mask_lo = 0x00ffffffUL; // 2^24 - 1 |
---|
40 | static const unsigned long int mask_hi = ~0x00ffffffUL; |
---|
41 | static const unsigned long int two24 = 16777216; // 2^24 |
---|
42 | |
---|
43 | |
---|
44 | // internal generator structure |
---|
45 | //------------------------------ |
---|
46 | typedef struct { |
---|
47 | unsigned int i; |
---|
48 | unsigned int j; |
---|
49 | unsigned int n; |
---|
50 | unsigned int skip; |
---|
51 | unsigned int carry; |
---|
52 | unsigned long int u[24]; |
---|
53 | } ranlux_state_t; |
---|
54 | |
---|
55 | |
---|
56 | // internal generator state |
---|
57 | //-------------------------- |
---|
58 | ranlux_state_t local_ranlux_state; |
---|
59 | |
---|
60 | |
---|
61 | // incrementation of the generator state |
---|
62 | //--------------------------------------- |
---|
63 | static inline unsigned long int increment_state(){ |
---|
64 | unsigned int i = local_ranlux_state.i; |
---|
65 | unsigned int j = local_ranlux_state.j; |
---|
66 | long int delta = local_ranlux_state.u[j] - local_ranlux_state.u[i] |
---|
67 | - local_ranlux_state.carry; |
---|
68 | |
---|
69 | if (delta & mask_hi){ |
---|
70 | local_ranlux_state.carry = 1; |
---|
71 | delta &= mask_lo; |
---|
72 | } else { |
---|
73 | local_ranlux_state.carry = 0; |
---|
74 | } |
---|
75 | |
---|
76 | local_ranlux_state.u[i] = delta; |
---|
77 | |
---|
78 | if (i==0) |
---|
79 | i = 23; |
---|
80 | else |
---|
81 | i--; |
---|
82 | |
---|
83 | local_ranlux_state.i = i; |
---|
84 | |
---|
85 | if (j == 0) |
---|
86 | j = 23; |
---|
87 | else |
---|
88 | j--; |
---|
89 | |
---|
90 | local_ranlux_state.j = j; |
---|
91 | |
---|
92 | return delta; |
---|
93 | } |
---|
94 | |
---|
95 | |
---|
96 | // set generator state |
---|
97 | //--------------------- |
---|
98 | static void ranlux_set(unsigned long int s){ |
---|
99 | int i; |
---|
100 | long int seed; |
---|
101 | |
---|
102 | if (s==0) |
---|
103 | s = 314159265; /* default seed is 314159265 */ |
---|
104 | |
---|
105 | seed = s; |
---|
106 | |
---|
107 | /* This is the initialization algorithm of F. James, widely in use |
---|
108 | for RANLUX. */ |
---|
109 | |
---|
110 | for (i=0;i<24;i++){ |
---|
111 | unsigned long int k = seed/53668; |
---|
112 | seed = 40014*(seed-k*53668)-k*12211; |
---|
113 | if (seed<0){ |
---|
114 | seed += 2147483563; |
---|
115 | } |
---|
116 | local_ranlux_state.u[i] = seed%two24; |
---|
117 | } |
---|
118 | |
---|
119 | local_ranlux_state.i = 23; |
---|
120 | local_ranlux_state.j = 9; |
---|
121 | local_ranlux_state.n = 0; |
---|
122 | local_ranlux_state.skip = 389-24; // 389 => best decorrelation |
---|
123 | |
---|
124 | if (local_ranlux_state.u[23]&mask_hi){ |
---|
125 | local_ranlux_state.carry = 1; |
---|
126 | } else { |
---|
127 | local_ranlux_state.carry = 0; |
---|
128 | } |
---|
129 | } |
---|
130 | |
---|
131 | |
---|
132 | // generator initialization |
---|
133 | //-------------------------- |
---|
134 | void ranlux_init(){ |
---|
135 | // seed the generator |
---|
136 | ranlux_set(0); |
---|
137 | } |
---|
138 | |
---|
139 | |
---|
140 | // get random number |
---|
141 | //------------------- |
---|
142 | unsigned long int ranlux_get(){ |
---|
143 | const unsigned int skip = local_ranlux_state.skip; |
---|
144 | unsigned long int r = increment_state(); |
---|
145 | |
---|
146 | local_ranlux_state.n++; |
---|
147 | |
---|
148 | if (local_ranlux_state.n == 24){ |
---|
149 | unsigned int i; |
---|
150 | local_ranlux_state.n = 0; |
---|
151 | for (i = 0; i < skip; i++) |
---|
152 | increment_state(); |
---|
153 | } |
---|
154 | |
---|
155 | return r; |
---|
156 | } |
---|
157 | |
---|
158 | // print generator state |
---|
159 | //----------------------- |
---|
160 | void ranlux_print_state(){ |
---|
161 | size_t i; |
---|
162 | unsigned char *p = (unsigned char *) (&local_ranlux_state); |
---|
163 | const size_t n = sizeof (ranlux_state_t); |
---|
164 | |
---|
165 | for (i=0;i<n;i++){ |
---|
166 | /* FIXME: we're assuming that a char is 8 bits */ |
---|
167 | printf("%.2x", *(p+i)); |
---|
168 | } |
---|
169 | } |
---|
170 | |
---|
171 | } |
---|