NAMD
Random.h
Go to the documentation of this file.
1 
7 /*
8  * Copyright (c) 1993 Martin Birgmeier
9  * All rights reserved.
10  *
11  * You may redistribute unmodified or modified versions of this source
12  * code provided that the above copyright notice and this and the
13  * following conditions are retained.
14  *
15  * This software is provided ``as is'', and comes with no warranties
16  * of any kind. I shall in no event be liable for anything that happens
17  * to anyone/anything when using this software.
18  */
19 
20 #ifndef RANDOM_H
21 #define RANDOM_H
22 
23 #include "common.h"
24 #include "Vector.h"
25 
26 #ifdef _MSC_VER
27 #define UINT64_LITERAL(X) X ## ui64
28 #else
29 #define UINT64_LITERAL(X) X ## ULL
30 #endif
31 
32 #define RAND48_SEED UINT64_LITERAL(0x00001234abcd330e)
33 #define RAND48_MULT UINT64_LITERAL(0x00000005deece66d)
34 #define RAND48_ADD UINT64_LITERAL(0x000000000000000b)
35 #define RAND48_MASK UINT64_LITERAL(0x0000ffffffffffff)
36 
37 class Random {
38 
39 private:
40 
41  double second_gaussian;
42  uint64_t second_gaussian_waiting;
43  uint64_t rand48_seed;
44  uint64_t rand48_mult;
45  uint64_t rand48_add;
46 
47 public:
48 
49  // default constructor
50  Random(void) {
51  init(0);
52  rand48_seed = RAND48_SEED;
53  }
54 
55  // constructor with seed
56  Random(unsigned long seed) {
57  init(seed);
58  }
59 
60  // reinitialize with seed
61  void init(unsigned long seed) {
62  second_gaussian = 0;
63  second_gaussian_waiting = 0;
64  rand48_seed = seed & UINT64_LITERAL(0x00000000ffffffff);
65  rand48_seed = rand48_seed << 16;
66  rand48_seed |= RAND48_SEED & UINT64_LITERAL(0x0000ffff);
67  rand48_mult = RAND48_MULT;
68  rand48_add = RAND48_ADD;
69  }
70 
71  // advance generator by one (seed = seed * mult + add, to 48 bits)
72  void skip(void) {
73  rand48_seed = ( rand48_seed * rand48_mult + rand48_add ) & RAND48_MASK;
74  }
75 
76  // split into numStreams different steams and take stream iStream
77  void split(int iStream, int numStreams) {
78 
79  int i;
80 
81  // make sure that numStreams is odd to ensure maximum period
82  numStreams |= 1;
83 
84  // iterate to get to the correct stream
85  for ( i = 0; i < iStream; ++i ) skip();
86 
87  // save seed and add so we can use skip() for our calculations
88  uint64_t save_seed = rand48_seed;
89 
90  // calculate c *= ( 1 + a + ... + a^(numStreams-1) )
91  rand48_seed = rand48_add;
92  for ( i = 1; i < numStreams; ++i ) skip();
93  uint64_t new_add = rand48_seed;
94 
95  // calculate a = a^numStreams
96  rand48_seed = rand48_mult;
97  rand48_add = 0;
98  for ( i = 1; i < numStreams; ++i ) skip();
99  rand48_mult = rand48_seed;
100 
101  rand48_add = new_add;
102  rand48_seed = save_seed;
103 
104  second_gaussian = 0;
105  second_gaussian_waiting = 0;
106  }
107 
108  // return a number uniformly distributed between 0 and 1
109  BigReal uniform(void) {
110  skip();
111  const double exp48 = ( 1.0 / (double)(UINT64_LITERAL(1) << 48) );
112  return ( (double) rand48_seed * exp48 );
113  }
114 
115  // return a number from a standard gaussian distribution
117 #ifdef DEBUG_RANDOM
118  static int count = 0;
119 #endif
120  BigReal fac, r, v1, v2;
121 
122  if (second_gaussian_waiting) {
123  second_gaussian_waiting = 0;
124 #ifdef DEBUG_RANDOM
125  count++;
126  fprintf(stderr, "RANDOM %6d %12.9f\n", count, second_gaussian);
127 #endif
128  return second_gaussian;
129  } else {
130  r = 2.; // r >= 1.523e-8 ensures abs result < 6
131  while (r >=1. || r < 1.523e-8) { // make sure we are within unit circle
132  v1 = 2.0 * uniform() - 1.0;
133  v2 = 2.0 * uniform() - 1.0;
134  r = v1*v1 + v2*v2;
135  }
136  fac = sqrt(-2.0 * log(r)/r);
137  // now make the Box-Muller transformation to get two normally
138  // distributed random numbers. Save one and return the other.
139  second_gaussian_waiting = 1;
140  second_gaussian = v1 * fac;
141 #ifdef DEBUG_RANDOM
142  count++;
143  fprintf(stderr, "RANDOM %6d %12.9f\n", count, v2 * fac);
144 #endif
145  return v2 * fac;
146  }
147  }
148 
189  BigReal sum_of_squared_gaussians(int64_t num_gaussians) {
190  BigReal z, u, v;
191 
192  if (num_gaussians <= 2) {
193  v = 0.0;
194  for(int i=0; i<num_gaussians; ++i) {
195  z = gaussian();
196  v += z*z;
197  }
198  return v;
199  } else if (2 < num_gaussians && num_gaussians <= 200) {
200  const BigReal d = 0.5*num_gaussians - 1./3;
201  const BigReal c = 1. / (3*sqrt(d));
202  const BigReal zmin = -1. / c;
203  do {
204  do {
205  z = gaussian();
206  }
207  while(z <= zmin);
208  u = uniform();
209  v = (1 + c*z); v = v*v*v; // v = (1 + c*z)^3
210  }
211  while(log(u) >= (0.5*z*z + d*(1 - v + log(v))));
212  return 2*d*v;
213  } else {
214  return num_gaussians + sqrt(2*num_gaussians)*gaussian();
215  }
216  }
217 
218  // return a vector of gaussian random numbers
220  return Vector( gaussian(), gaussian(), gaussian() );
221  }
222 
223  // return a random long
224  // signed int32 ranges from 0 to (2^31)-1
225  long integer(void) {
226  skip();
227  return ( ( rand48_seed >> 17 ) & UINT64_LITERAL(0x000000007fffffff) );
228  }
229 
230  // randomly order an array of whatever
231  template <class Elem> void reorder(Elem *a, int n) {
232  for ( int i = 0; i < (n-1); ++i ) {
233  int ie = i + ( integer() % (n-i) );
234  if ( ie == i ) continue;
235  const Elem e = a[ie];
236  a[ie] = a[i];
237  a[i] = e;
238  }
239  }
240 
244  void uniform_array_f(float *a, int n) {
245  // for now just call uniform()
246  // ultimately we will provide faster implementation
247  for (int i=0; i < n; i++) {
248  a[i] = (float)uniform();
249  }
250  }
251 
255  void gaussian_array_f(float *a, int n) {
256  // for now just call gaussian()
257  // ultimately we will provide faster implementation
258  for (int i=0; i < n; i++) {
259  a[i] = (float)gaussian();
260  }
261  }
262 
263 };
264 
265 #endif // RANDOM_H
266 
Vector gaussian_vector(void)
Definition: Random.h:219
BigReal uniform(void)
Definition: Random.h:109
BigReal gaussian(void)
Definition: Random.h:116
Random(void)
Definition: Random.h:50
void skip(void)
Definition: Random.h:72
Definition: Vector.h:72
BigReal sum_of_squared_gaussians(int64_t num_gaussians)
Return the sum of i.i.d. squared Gaussians.
Definition: Random.h:189
void uniform_array_f(float *a, int n)
Definition: Random.h:244
void split(int iStream, int numStreams)
Definition: Random.h:77
void reorder(Elem *a, int n)
Definition: Random.h:231
#define RAND48_ADD
Definition: Random.h:34
void init(unsigned long seed)
Definition: Random.h:61
Definition: Random.h:37
long integer(void)
Definition: Random.h:225
Random(unsigned long seed)
Definition: Random.h:56
void gaussian_array_f(float *a, int n)
Definition: Random.h:255
#define RAND48_SEED
Definition: Random.h:32
#define RAND48_MULT
Definition: Random.h:33
#define UINT64_LITERAL(X)
Definition: Random.h:29
#define RAND48_MASK
Definition: Random.h:35
double BigReal
Definition: common.h:123