Line data Source code
1 : // Random.cpp is based on Random.h
2 : // Mersenne Twister random number generator -- a C++ class Random
3 : // Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
4 : // Richard J. Wagner v1.0 15 May 2003 rjwagner@writeme.com
5 :
6 : // The Mersenne Twister is an algorithm for generating random numbers. It
7 : // was designed with consideration of the flaws in various other generators.
8 : // The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
9 : // are far greater. The generator is also fast; it avoids multiplication and
10 : // division, and it benefits from caches and pipelines. For more information
11 : // see the inventors' web page at http://www.math.keio.ac.jp/~matumoto/emt.html
12 :
13 : // Reference
14 : // M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
15 : // Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
16 : // Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
17 :
18 : // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
19 : // Copyright (C) 2000 - 2003, Richard J. Wagner
20 : // All rights reserved.
21 : //
22 : // Redistribution and use in source and binary forms, with or without
23 : // modification, are permitted provided that the following conditions
24 : // are met:
25 : //
26 : // 1. Redistributions of source code must retain the above copyright
27 : // notice, this list of conditions and the following disclaimer.
28 : //
29 : // 2. Redistributions in binary form must reproduce the above copyright
30 : // notice, this list of conditions and the following disclaimer in the
31 : // documentation and/or other materials provided with the distribution.
32 : //
33 : // 3. The names of its contributors may not be used to endorse or promote
34 : // products derived from this software without specific prior written
35 : // permission.
36 : //
37 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
38 : // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
39 : // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
40 : // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
42 : // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
43 : // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
44 : // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
45 : // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
46 : // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
47 : // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
48 :
49 : // The original code included the following notice:
50 : //
51 : // When you use this, send an email to: matumoto@math.keio.ac.jp
52 : // with an appropriate reference to your work.
53 : //
54 : // It would be nice to CC: rjwagner@writeme.com and Cokus@math.washington.edu
55 : // when you write.
56 :
57 : // Parts of this file are modified beginning in 29.10.09 for adaption in PXL.
58 : // Parts of this file are modified beginning in 10.02.12 for adaption in CRPropa.
59 :
60 : #include "crpropa/Random.h"
61 :
62 : #include "crpropa/base64.h"
63 :
64 : #include <cstdio>
65 :
66 : namespace crpropa {
67 :
68 0 : Random::Random(const uint32_t& oneSeed) {
69 0 : seed(oneSeed);
70 0 : }
71 :
72 0 : Random::Random(uint32_t * const bigSeed, const uint32_t seedLength) {
73 0 : seed(bigSeed, seedLength);
74 0 : }
75 :
76 4625 : Random::Random() {
77 4625 : seed();
78 4625 : }
79 :
80 22417505 : double Random::rand() {
81 22417505 : return double(randInt()) * (1.0 / 4294967295.0);
82 : }
83 :
84 0 : double Random::rand(const double& n) {
85 0 : return rand() * n;
86 : }
87 :
88 1715954 : double Random::randExc() {
89 1715954 : return double(randInt()) * (1.0 / 4294967296.0);
90 : }
91 :
92 0 : double Random::randExc(const double& n) {
93 0 : return randExc() * n;
94 : }
95 :
96 1715954 : double Random::randDblExc() {
97 1715954 : return (double(randInt()) + 0.5) * (1.0 / 4294967296.0);
98 : }
99 :
100 0 : double Random::randDblExc(const double& n) {
101 0 : return randDblExc() * n;
102 : }
103 :
104 0 : double Random::rand53() {
105 0 : uint32_t a = randInt() >> 5, b = randInt() >> 6;
106 0 : return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); // by Isaku Wada
107 : }
108 :
109 : // Return a real number from a normal (Gaussian) distribution with given
110 : // mean and variance by Box-Muller method
111 1715954 : double Random::randNorm(const double& mean, const double& variance) {
112 1715954 : double r = sqrt(-2.0 * log(1.0 - randDblExc())) * variance;
113 1715954 : double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
114 1715954 : return mean + r * cos(phi);
115 : }
116 :
117 2404687 : double Random::randUniform(double min, double max) {
118 2404687 : return min + (max - min) * rand();
119 : }
120 :
121 0 : double Random::randRayleigh(double sigma) {
122 0 : return sigma * sqrt(-2.0 * log(1 - rand()));
123 : }
124 :
125 1000 : double Random::randFisher(double kappa) {
126 1000 : return acos(1. + 1. / kappa * log(1 - rand() * (1 - exp(-2 * kappa))));
127 : }
128 :
129 10102 : size_t Random::randBin(const std::vector<float> &cdf) {
130 10102 : std::vector<float>::const_iterator it = std::lower_bound(cdf.begin(),
131 10102 : cdf.end(), rand() * cdf.back());
132 10102 : return it - cdf.begin();
133 : }
134 :
135 3744701 : size_t Random::randBin(const std::vector<double> &cdf) {
136 3744701 : std::vector<double>::const_iterator it = std::lower_bound(cdf.begin(),
137 3744701 : cdf.end(), rand() * cdf.back());
138 3744701 : return it - cdf.begin();
139 : }
140 :
141 962203 : Vector3d Random::randVector() {
142 962203 : double z = randUniform(-1.0, 1.0);
143 962203 : double t = randUniform(-1.0 * M_PI, M_PI);
144 962203 : double r = sqrt(1 - z * z);
145 962203 : return Vector3d(r * cos(t), r * sin(t), z);
146 : }
147 :
148 481001 : Vector3d Random::randVectorAroundMean(const Vector3d &meanDirection,
149 : double angle) {
150 481001 : Vector3d axis = meanDirection.cross(randVector());
151 : Vector3d v = meanDirection;
152 481001 : return v.getRotated(axis, angle);
153 : }
154 :
155 1000 : Vector3d Random::randFisherVector(const Vector3d &meanDirection, double kappa) {
156 1000 : return randVectorAroundMean(meanDirection, randFisher(kappa));
157 : }
158 :
159 480001 : Vector3d Random::randConeVector(const Vector3d &meanDirection, double angularRadius) {
160 480001 : const double theta = acos(randUniform(1, cos(angularRadius)));
161 480001 : return randVectorAroundMean(meanDirection, theta);
162 : }
163 :
164 0 : Vector3d Random::randVectorLamberts() {
165 : // random vector following Lamberts cosine law (https://en.wikipedia.org/wiki/Lambert%27s_cosine_law)
166 : // for a surface element with normal vector pointing in positive z-axis (0, 0, 1)
167 0 : double phi = randUniform(-1.0 * M_PI, M_PI);
168 0 : double theta = M_PI / 2.0 - acos(sqrt(randUniform(0, 1)));
169 0 : return Vector3d(cos(phi) * cos(theta), sin(phi) * cos(theta), sin(theta));
170 : }
171 :
172 0 : Vector3d Random::randVectorLamberts(const Vector3d &normalVector) {
173 : // random vector following Lamberts cosine law for a surface element described by normalVector
174 0 : Vector3d vLambertz = randVectorLamberts();
175 : // find rotation axis that rotates the z-axis to the normalVector of the surface element
176 : Vector3d axis = normalVector.cross(Vector3d(0, 0, 1));
177 0 : if (axis.getR() < std::numeric_limits<double>::epsilon()) {
178 : axis = Vector3d(0, 0, 1);
179 : }
180 0 : double angle = normalVector.getAngleTo(Vector3d(0, 0, 1));
181 : // rotate the random Lamberts vector from z-axis to respective surface element
182 0 : return vLambertz.getRotated(axis / axis.getR(), -angle);
183 : }
184 :
185 3634120 : Vector3d Random::randomInterpolatedPosition(const Vector3d &a, const Vector3d &b) {
186 3634120 : return a + rand() * (b - a);
187 : }
188 :
189 1104 : double Random::randPowerLaw(double index, double min, double max) {
190 1104 : if ((min < 0) || (max < min)) {
191 : throw std::runtime_error(
192 0 : "Power law distribution only possible for 0 <= min <= max");
193 : }
194 : //check for index -1!
195 1104 : if ((std::abs(index + 1.0)) < std::numeric_limits<double>::epsilon()) {
196 2 : double part1 = log(max);
197 2 : double part2 = log(min);
198 2 : return exp((part1 - part2) * rand() + part2);
199 : } else {
200 1102 : double part1 = pow(max, index + 1);
201 1102 : double part2 = pow(min, index + 1);
202 1102 : double ex = 1 / (index + 1);
203 1102 : return pow((part1 - part2) * rand() + part2, ex);
204 : }
205 : }
206 :
207 0 : double Random::randBrokenPowerLaw(double index1, double index2,
208 : double breakpoint, double min, double max) {
209 0 : if ((min <= 0) || (max < min)) {
210 : throw std::runtime_error(
211 0 : "Power law distribution only possible for 0 < min <= max");
212 : }
213 0 : if (min >= breakpoint) {
214 0 : return this->randPowerLaw(index2, min, max);
215 0 : } else if (max <= breakpoint) {
216 0 : return this->randPowerLaw(index2, min, max);
217 : } else {
218 : double intPL1;
219 : // check if index1 = -1
220 0 : if ((std::abs(index1 + 1.0)) < std::numeric_limits<double>::epsilon()) {
221 0 : intPL1 = log(breakpoint / min);
222 : } else {
223 0 : intPL1 = (pow(breakpoint, index1 + 1) - pow(min, index1 + 1))
224 : / (index1 + 1);
225 : }
226 : double intPL2;
227 : // check if index2 = -1
228 0 : if ((std::abs(index2 + 1.0)) < std::numeric_limits<double>::epsilon()) {
229 0 : intPL2 = log(max / breakpoint) * pow(breakpoint, index1 - index2);
230 : } else {
231 0 : intPL2 = (pow(max, index2 + 1) - pow(breakpoint, index2 + 1))
232 0 : * pow(breakpoint, index1 - index2) / (index2 + 1);
233 : }
234 0 : if (rand() > intPL1 / (intPL1 + intPL2))
235 0 : return randPowerLaw(index2, breakpoint, max);
236 : else
237 0 : return randPowerLaw(index1, min, breakpoint);
238 : }
239 : }
240 :
241 0 : double Random::randExponential() {
242 : double dum;
243 : do {
244 0 : dum = rand();
245 0 : } while (dum < std::numeric_limits<double>::epsilon());
246 0 : return -1.0 * log(dum);
247 : }
248 :
249 25854363 : uint32_t Random::randInt() {
250 25854363 : if (left == 0)
251 41427 : reload();
252 25854363 : --left;
253 :
254 : uint32_t s1;
255 25854363 : s1 = *pNext++;
256 25854363 : s1 ^= (s1 >> 11);
257 25854363 : s1 ^= (s1 << 7) & 0x9d2c5680UL;
258 25854363 : s1 ^= (s1 << 15) & 0xefc60000UL;
259 25854363 : return (s1 ^ (s1 >> 18));
260 : }
261 :
262 0 : uint32_t Random::randInt(const uint32_t& n) {
263 : // Find which bits are used in n
264 : // Optimized by Magnus Jonsson (magnus@smartelectronix.com)
265 0 : uint32_t used = n;
266 0 : used |= used >> 1;
267 0 : used |= used >> 2;
268 0 : used |= used >> 4;
269 0 : used |= used >> 8;
270 0 : used |= used >> 16;
271 :
272 : // Draw numbers until one is found in [0,n]
273 : uint32_t i;
274 : do
275 0 : i = randInt() & used; // toss unused bits to shorten search
276 0 : while (i > n);
277 0 : return i;
278 : }
279 :
280 :
281 0 : uint64_t Random::randInt64()
282 : {
283 0 : int64_t a = randInt();
284 0 : int64_t b = randInt();
285 : // a is shifted to the left to create a proper 64 bit integer
286 0 : return (b + (a << 32));
287 : }
288 :
289 :
290 0 : uint64_t Random::randInt64(const uint64_t &n)
291 : {
292 0 : uint64_t used = n;
293 0 : used |= used >> 1;
294 0 : used |= used >> 2;
295 0 : used |= used >> 4;
296 0 : used |= used >> 8;
297 0 : used |= used >> 16;
298 0 : used |= used >> 32;
299 :
300 : // Draw numbers until one is found in [0,n]
301 : uint64_t i;
302 : do
303 0 : i = randInt64() & used; // toss unused bits to shorten search
304 0 : while (i > n);
305 0 : return i;
306 : }
307 :
308 :
309 :
310 13 : void Random::seed(const uint32_t oneSeed) {
311 13 : initial_seed.resize(1);
312 13 : initial_seed[0] = oneSeed;
313 13 : initialize(oneSeed);
314 13 : reload();
315 13 : }
316 :
317 4628 : void Random::seed(uint32_t * const bigSeed, const uint32_t seedLength) {
318 :
319 4628 : initial_seed.resize(seedLength);
320 2891260 : for (size_t i =0; i< seedLength; i++)
321 : {
322 2886632 : initial_seed[i] = bigSeed[i];
323 : }
324 :
325 4628 : initialize(19650218UL);
326 : int i = 1;
327 : uint32_t j = 0;
328 4628 : int k = (N > seedLength ? N : seedLength);
329 2892500 : for (; k; --k) {
330 2887872 : state[i] = state[i]
331 2887872 : ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL);
332 2887872 : state[i] += (bigSeed[j] & 0xffffffffUL) + j;
333 : state[i] &= 0xffffffffUL;
334 2887872 : ++i;
335 2887872 : ++j;
336 2887872 : if (i >= N) {
337 4628 : state[0] = state[N - 1];
338 : i = 1;
339 : }
340 2887872 : if (j >= seedLength)
341 : j = 0;
342 : }
343 2887872 : for (k = N - 1; k; --k) {
344 2883244 : state[i] = state[i]
345 2883244 : ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL);
346 2883244 : state[i] -= i;
347 : state[i] &= 0xffffffffUL;
348 2883244 : ++i;
349 2883244 : if (i >= N) {
350 4628 : state[0] = state[N - 1];
351 : i = 1;
352 : }
353 : }
354 4628 : state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
355 4628 : reload();
356 4628 : }
357 :
358 4625 : void Random::seed() {
359 : // First try getting an array from /dev/urandom
360 4625 : FILE* urandom = std::fopen("/dev/urandom", "rb");
361 4625 : if (urandom) {
362 : uint32_t bigSeed[N];
363 : uint32_t *s = bigSeed;
364 : int i = N;
365 : bool success = true;
366 2890625 : while (success && i--)
367 5772000 : success = std::fread(s++, sizeof(uint32_t), 1, urandom) != 0;
368 4625 : std::fclose(urandom);
369 4625 : if (success) {
370 4625 : seed(bigSeed, N);
371 4625 : return;
372 : }
373 : }
374 :
375 : // Was not successful, so use time() and clock() instead
376 0 : seed(hash(time(NULL), clock()));
377 : }
378 :
379 :
380 4641 : void Random::initialize(const uint32_t seed) {
381 : uint32_t *s = state;
382 4641 : uint32_t *r = state;
383 : int i = 1;
384 4641 : *s++ = seed & 0xffffffffUL;
385 2895984 : for (; i < N; ++i) {
386 2891343 : *s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
387 2891343 : r++;
388 : }
389 4641 : }
390 :
391 46068 : void Random::reload() {
392 46068 : uint32_t *p = state;
393 : int i;
394 10503504 : for (i = +N - M; i--; ++p)
395 10457436 : *p = twist(p[M], p[0], p[1]);
396 18288996 : for (i = M; --i; ++p)
397 18242928 : *p = twist(p[+M - N], p[0], p[1]);
398 46068 : *p = twist(p[+M - N], p[0], state[0]);
399 :
400 46068 : left = N, pNext = state;
401 46068 : }
402 :
403 0 : uint32_t Random::hash(time_t t, clock_t c) {
404 : static uint32_t differ = 0; // guarantee time-based seeds will change
405 :
406 : uint32_t h1 = 0;
407 : unsigned char *p = (unsigned char *) &t;
408 0 : for (size_t i = 0; i < sizeof(t); ++i) {
409 0 : h1 *= std::numeric_limits<unsigned char>::max() + 2U;
410 0 : h1 += p[i];
411 : }
412 : uint32_t h2 = 0;
413 : p = (unsigned char *) &c;
414 0 : for (size_t j = 0; j < sizeof(c); ++j) {
415 0 : h2 *= std::numeric_limits<unsigned char>::max() + 2U;
416 0 : h2 += p[j];
417 : }
418 0 : return (h1 + differ++) ^ h2;
419 : }
420 :
421 0 : void Random::save(uint32_t* saveArray) const {
422 : uint32_t *sa = saveArray;
423 0 : const uint32_t *s = state;
424 : int i = N;
425 0 : for (; i--; *sa++ = *s++) {
426 : }
427 0 : *sa = left;
428 0 : }
429 :
430 2 : const std::vector<uint32_t> &Random::getSeed() const
431 : {
432 2 : return initial_seed;
433 : }
434 :
435 0 : void Random::load(uint32_t * const loadArray) {
436 0 : uint32_t *s = state;
437 : uint32_t *la = loadArray;
438 : int i = N;
439 0 : for (; i--; *s++ = *la++) {
440 : }
441 0 : left = *la;
442 0 : pNext = &state[N - left];
443 0 : }
444 :
445 0 : std::ostream& operator<<(std::ostream& os, const Random& mtrand) {
446 0 : const uint32_t *s = mtrand.state;
447 : int i = mtrand.N;
448 0 : for (; i--; os << *s++ << "\t") {
449 : }
450 0 : return os << mtrand.left;
451 : }
452 :
453 0 : std::istream& operator>>(std::istream& is, Random& mtrand) {
454 0 : uint32_t *s = mtrand.state;
455 : int i = mtrand.N;
456 0 : for (; i--; is >> *s++) {
457 : }
458 0 : is >> mtrand.left;
459 0 : mtrand.pNext = &mtrand.state[mtrand.N - mtrand.left];
460 0 : return is;
461 : }
462 :
463 : #ifdef _OPENMP
464 : #include <omp.h>
465 : #include <stdexcept>
466 :
467 : // see http://stackoverflow.com/questions/8051108/using-the-openmp-threadprivate-directive-on-static-instances-of-c-stl-types
468 : const static int MAX_THREAD = 256;
469 :
470 : struct RANDOM_TLS_ITEM {
471 : Random r;
472 : char padding[(sizeof(Random) / 64 + 1) * 64 - sizeof(Random)];
473 : };
474 :
475 : #ifdef _MSC_VER
476 : __declspec(align(64)) static RANDOM_TLS_ITEM _tls[MAX_THREAD];
477 : #else
478 : __attribute__ ((aligned(64))) static RANDOM_TLS_ITEM _tls[MAX_THREAD];
479 : #endif
480 :
481 2629024 : Random &Random::instance() {
482 2629024 : int i = omp_get_thread_num();
483 2629024 : if (i >= MAX_THREAD)
484 0 : throw std::runtime_error("crpropa::Random: more than MAX_THREAD threads!");
485 2629024 : return _tls[i].r;
486 : }
487 :
488 0 : void Random::seedThreads(const uint32_t oneSeed) {
489 0 : for(size_t i = 0; i < MAX_THREAD; ++i)
490 0 : _tls[i].r.seed(oneSeed + i);
491 0 : }
492 :
493 0 : std::vector< std::vector<uint32_t> > Random::getSeedThreads()
494 : {
495 : std::vector< std::vector<uint32_t> > seeds;
496 0 : for(size_t i = 0; i < omp_get_num_threads(); ++i)
497 0 : seeds.push_back(_tls[i].r.getSeed() );
498 0 : return seeds;
499 0 : }
500 :
501 : #else
502 : static Random _random;
503 : Random &Random::instance() {
504 : return _random;
505 : }
506 : void Random::seedThreads(const uint32_t oneSeed) {
507 : _random.seed(oneSeed);
508 : }
509 : std::vector< std::vector<uint32_t> > Random::getSeedThreads()
510 : {
511 : std::vector< std::vector<uint32_t> > seeds;
512 : seeds.push_back(_random.getSeed() );
513 : return seeds;
514 : }
515 : #endif
516 :
517 0 : const std::string Random::getSeed_base64() const
518 : {
519 0 : return Base64::encode((unsigned char*) &initial_seed[0], sizeof(initial_seed[0]) * initial_seed.size() / sizeof(unsigned char));
520 : }
521 :
522 1 : void Random::seed(const std::string &b64Seed)
523 : {
524 1 : std::string decoded_data = Base64::decode(b64Seed);
525 1 : size_t seedSize = decoded_data.size() * sizeof(decoded_data[0]) / sizeof(uint32_t);
526 1 : seed((uint32_t*)decoded_data.c_str(), seedSize );
527 1 : }
528 :
529 : } // namespace crpropa
530 :
|
