SISCone  2.0.6
ranlux.cpp
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];
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 }
Definition: area.cpp:33
The SISCone project has been developed by Gavin Salam and Gregory Soyez
Documentation generated on Sun Sep 20 2015 01:33:13 for SISCone by  Doxygen 1.8.9.1