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 4881 : Random::Random() {
77 4881 : seed();
78 4881 : }
79 :
80 22563868 : double Random::rand() {
81 22563868 : 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 : 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 3742320 : size_t Random::randBin(const std::vector<double> &cdf) {
136 : std::vector<double>::const_iterator it = std::lower_bound(cdf.begin(),
137 3742320 : cdf.end(), rand() * cdf.back());
138 3742320 : 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 3631825 : Vector3d Random::randomInterpolatedPosition(const Vector3d &a, const Vector3d &b) {
186 3631825 : 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 0 : 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 0 : 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 26002486 : uint32_t Random::randInt() {
250 26002486 : if (left == 0)
251 41663 : reload();
252 26002486 : --left;
253 :
254 : uint32_t s1;
255 26002486 : s1 = *pNext++;
256 26002486 : s1 ^= (s1 >> 11);
257 26002486 : s1 ^= (s1 << 7) & 0x9d2c5680UL;
258 26002486 : s1 ^= (s1 << 15) & 0xefc60000UL;
259 26002486 : 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 880 : uint64_t Random::randInt64()
282 : {
283 880 : int64_t a = randInt();
284 880 : int64_t b = randInt();
285 880 : return (b + a << 32);
286 : }
287 :
288 :
289 421 : uint64_t Random::randInt64(const uint64_t &n)
290 : {
291 421 : uint64_t used = n;
292 421 : used |= used >> 1;
293 421 : used |= used >> 2;
294 421 : used |= used >> 4;
295 421 : used |= used >> 8;
296 421 : used |= used >> 16;
297 421 : used |= used >> 32;
298 :
299 : // Draw numbers until one is found in [0,n]
300 : uint64_t i;
301 : do
302 880 : i = randInt64() & used; // toss unused bits to shorten search
303 880 : while (i > n);
304 421 : return i;
305 : }
306 :
307 :
308 :
309 13 : void Random::seed(const uint32_t oneSeed) {
310 13 : initial_seed.resize(1);
311 13 : initial_seed[0] = oneSeed;
312 13 : initialize(oneSeed);
313 13 : reload();
314 13 : }
315 :
316 4884 : void Random::seed(uint32_t * const bigSeed, const uint32_t seedLength) {
317 :
318 4884 : initial_seed.resize(seedLength);
319 3051260 : for (size_t i =0; i< seedLength; i++)
320 : {
321 3046376 : initial_seed[i] = bigSeed[i];
322 : }
323 :
324 4884 : initialize(19650218UL);
325 : int i = 1;
326 : uint32_t j = 0;
327 4884 : int k = (N > seedLength ? N : seedLength);
328 3052500 : for (; k; --k) {
329 3047616 : state[i] = state[i]
330 3047616 : ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL);
331 3047616 : state[i] += (bigSeed[j] & 0xffffffffUL) + j;
332 : state[i] &= 0xffffffffUL;
333 3047616 : ++i;
334 3047616 : ++j;
335 3047616 : if (i >= N) {
336 4884 : state[0] = state[N - 1];
337 : i = 1;
338 : }
339 3047616 : if (j >= seedLength)
340 : j = 0;
341 : }
342 3047616 : for (k = N - 1; k; --k) {
343 3042732 : state[i] = state[i]
344 3042732 : ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL);
345 3042732 : state[i] -= i;
346 : state[i] &= 0xffffffffUL;
347 3042732 : ++i;
348 3042732 : if (i >= N) {
349 4884 : state[0] = state[N - 1];
350 : i = 1;
351 : }
352 : }
353 4884 : state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
354 4884 : reload();
355 4884 : }
356 :
357 4881 : void Random::seed() {
358 : // First try getting an array from /dev/urandom
359 4881 : FILE* urandom = std::fopen("/dev/urandom", "rb");
360 4881 : if (urandom) {
361 : uint32_t bigSeed[N];
362 : uint32_t *s = bigSeed;
363 : int i = N;
364 : bool success = true;
365 3050625 : while (success && i--)
366 6091488 : success = std::fread(s++, sizeof(uint32_t), 1, urandom) != 0;
367 4881 : std::fclose(urandom);
368 4881 : if (success) {
369 4881 : seed(bigSeed, N);
370 4881 : return;
371 : }
372 : }
373 :
374 : // Was not successful, so use time() and clock() instead
375 0 : seed(hash(time(NULL), clock()));
376 : }
377 :
378 :
379 4897 : void Random::initialize(const uint32_t seed) {
380 4897 : uint32_t *s = state;
381 : uint32_t *r = state;
382 : int i = 1;
383 4897 : *s++ = seed & 0xffffffffUL;
384 3055728 : for (; i < N; ++i) {
385 3050831 : *s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
386 3050831 : r++;
387 : }
388 4897 : }
389 :
390 46560 : void Random::reload() {
391 46560 : uint32_t *p = state;
392 : int i;
393 10615680 : for (i = N - M; i--; ++p)
394 10569120 : *p = twist(p[M], p[0], p[1]);
395 18484320 : for (i = M; --i; ++p)
396 18437760 : *p = twist(p[M - N], p[0], p[1]);
397 46560 : *p = twist(p[M - N], p[0], state[0]);
398 :
399 46560 : left = N, pNext = state;
400 46560 : }
401 :
402 0 : uint32_t Random::hash(time_t t, clock_t c) {
403 : static uint32_t differ = 0; // guarantee time-based seeds will change
404 :
405 : uint32_t h1 = 0;
406 : unsigned char *p = (unsigned char *) &t;
407 0 : for (size_t i = 0; i < sizeof(t); ++i) {
408 0 : h1 *= std::numeric_limits<unsigned char>::max() + 2U;
409 0 : h1 += p[i];
410 : }
411 : uint32_t h2 = 0;
412 : p = (unsigned char *) &c;
413 0 : for (size_t j = 0; j < sizeof(c); ++j) {
414 0 : h2 *= std::numeric_limits<unsigned char>::max() + 2U;
415 0 : h2 += p[j];
416 : }
417 0 : return (h1 + differ++) ^ h2;
418 : }
419 :
420 0 : void Random::save(uint32_t* saveArray) const {
421 : uint32_t *sa = saveArray;
422 0 : const uint32_t *s = state;
423 : int i = N;
424 0 : for (; i--; *sa++ = *s++) {
425 : }
426 0 : *sa = left;
427 0 : }
428 :
429 2 : const std::vector<uint32_t> &Random::getSeed() const
430 : {
431 2 : return initial_seed;
432 : }
433 :
434 0 : void Random::load(uint32_t * const loadArray) {
435 0 : uint32_t *s = state;
436 : uint32_t *la = loadArray;
437 : int i = N;
438 0 : for (; i--; *s++ = *la++) {
439 : }
440 0 : left = *la;
441 0 : pNext = &state[N - left];
442 0 : }
443 :
444 0 : std::ostream& operator<<(std::ostream& os, const Random& mtrand) {
445 0 : const uint32_t *s = mtrand.state;
446 : int i = mtrand.N;
447 0 : for (; i--; os << *s++ << "\t") {
448 : }
449 0 : return os << mtrand.left;
450 : }
451 :
452 0 : std::istream& operator>>(std::istream& is, Random& mtrand) {
453 0 : uint32_t *s = mtrand.state;
454 : int i = mtrand.N;
455 0 : for (; i--; is >> *s++) {
456 : }
457 0 : is >> mtrand.left;
458 0 : mtrand.pNext = &mtrand.state[mtrand.N - mtrand.left];
459 0 : return is;
460 : }
461 :
462 : #ifdef _OPENMP
463 : #include <omp.h>
464 : #include <stdexcept>
465 :
466 : // see http://stackoverflow.com/questions/8051108/using-the-openmp-threadprivate-directive-on-static-instances-of-c-stl-types
467 : const static int MAX_THREAD = 256;
468 :
469 : struct RANDOM_TLS_ITEM {
470 : Random r;
471 : char padding[(sizeof(Random) / 64 + 1) * 64 - sizeof(Random)];
472 : };
473 :
474 : #ifdef _MSC_VER
475 : __declspec(align(64)) static RANDOM_TLS_ITEM _tls[MAX_THREAD];
476 : #else
477 : __attribute__ ((aligned(64))) static RANDOM_TLS_ITEM _tls[MAX_THREAD];
478 : #endif
479 :
480 2708313 : Random &Random::instance() {
481 2708313 : int i = omp_get_thread_num();
482 2708313 : if (i >= MAX_THREAD)
483 0 : throw std::runtime_error("crpropa::Random: more than MAX_THREAD threads!");
484 2708313 : return _tls[i].r;
485 : }
486 :
487 0 : void Random::seedThreads(const uint32_t oneSeed) {
488 0 : for(size_t i = 0; i < MAX_THREAD; ++i)
489 0 : _tls[i].r.seed(oneSeed + i);
490 0 : }
491 :
492 0 : std::vector< std::vector<uint32_t> > Random::getSeedThreads()
493 : {
494 : std::vector< std::vector<uint32_t> > seeds;
495 0 : for(size_t i = 0; i < omp_get_num_threads(); ++i)
496 0 : seeds.push_back(_tls[i].r.getSeed() );
497 0 : return seeds;
498 0 : }
499 :
500 : #else
501 : static Random _random;
502 : Random &Random::instance() {
503 : return _random;
504 : }
505 : void Random::seedThreads(const uint32_t oneSeed) {
506 : _random.seed(oneSeed);
507 : }
508 : std::vector< std::vector<uint32_t> > Random::getSeedThreads()
509 : {
510 : std::vector< std::vector<uint32_t> > seeds;
511 : seeds.push_back(_random.getSeed() );
512 : return seeds;
513 : }
514 : #endif
515 :
516 0 : const std::string Random::getSeed_base64() const
517 : {
518 0 : return Base64::encode((unsigned char*) &initial_seed[0], sizeof(initial_seed[0]) * initial_seed.size() / sizeof(unsigned char));
519 : }
520 :
521 1 : void Random::seed(const std::string &b64Seed)
522 : {
523 1 : std::string decoded_data = Base64::decode(b64Seed);
524 1 : size_t seedSize = decoded_data.size() * sizeof(decoded_data[0]) / sizeof(uint32_t);
525 1 : seed((uint32_t*)decoded_data.c_str(), seedSize );
526 1 : }
527 :
528 : } // namespace crpropa
529 :
|
