Line data Source code
1 : // 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 : #ifndef RANDOM_H
61 : #define RANDOM_H
62 :
63 : // Not thread safe (unless auto-initialization is avoided and each thread has
64 : // its own Random object)
65 : #include "crpropa/Vector3.h"
66 :
67 : #include <iostream>
68 : #include <limits>
69 : #include <ctime>
70 : #include <cmath>
71 : #include <vector>
72 : #include <stdexcept>
73 : #include <algorithm>
74 :
75 : #include <stdint.h>
76 : #include <string>
77 :
78 : //necessary for win32
79 : #ifndef M_PI
80 : #define M_PI 3.14159265358979323846
81 : #endif
82 :
83 : namespace crpropa {
84 :
85 : /**
86 : * \addtogroup Core
87 : * @{
88 : */
89 : /**
90 : @class Random
91 : @brief Random number generator.
92 :
93 : Mersenne Twister random number generator -- a C++ class Random
94 : Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
95 : Richard J. Wagner v1.0 15 May 2003 rjwagner\@writeme.com
96 : */
97 15 : class Random {
98 : public:
99 : enum {N = 624}; // length of state vector
100 : enum {SAVE = N + 1}; // length of array for save()
101 :
102 : protected:
103 : enum {M = 397}; // period parameter
104 : uint32_t state[N];// internal state
105 : std::vector<uint32_t> initial_seed;//
106 : uint32_t *pNext;// next value to get from state
107 : int left;// number of values left before reload needed
108 :
109 : //Methods
110 : public:
111 : /// initialize with a simple uint32_t
112 : Random( const uint32_t& oneSeed );
113 : // initialize with an array
114 : Random( uint32_t *const bigSeed, uint32_t const seedLength = N );
115 : /// auto-initialize with /dev/urandom or time() and clock()
116 : /// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
117 : /// values together, otherwise the generator state can be learned after
118 : /// reading 624 consecutive values.
119 : Random();
120 : // Access to 32-bit random numbers
121 : double rand();///< real number in [0,1]
122 : double rand( const double& n );///< real number in [0,n]
123 : double randExc();///< real number in [0,1)
124 : double randExc( const double& n );///< real number in [0,n)
125 : double randDblExc();///< real number in (0,1)
126 : double randDblExc( const double& n );///< real number in (0,n)
127 : // Pull a 32-bit integer from the generator state
128 : // Every other access function simply transforms the numbers extracted here
129 : uint32_t randInt();///< integer in [0,2**32-1]
130 : uint32_t randInt( const uint32_t& n );///< integer in [0,n] for n < 2**32
131 :
132 : uint64_t randInt64(); ///< integer in [0, 2**64 -1]. PROBABLY NOT SECURE TO USE
133 : uint64_t randInt64(const uint64_t &n); ///< integer in [0, n] for n < 2**64 -1. PROBABLY NOT SECURE TO USE
134 :
135 0 : double operator()() {return rand();} ///< same as rand()
136 :
137 : // Access to 53-bit random numbers (capacity of IEEE double precision)
138 : double rand53();///< real number in [0,1) (capacity of IEEE double precision)
139 : ///Exponential distribution in (0,inf)
140 : double randExponential();
141 : /// Normal distributed random number
142 1715954 : double randNorm( const double& mean = 0.0, const double& variance = 1.0 );
143 : /// Uniform distribution in [min, max]
144 : double randUniform(double min, double max);
145 : /// Rayleigh distributed random number
146 : double randRayleigh(double sigma);
147 : /// Fisher distributed random number
148 : double randFisher(double k);
149 :
150 : /// Draw a random bin from a (unnormalized) cumulative distribution function, without leading zero.
151 : size_t randBin(const std::vector<float> &cdf);
152 : size_t randBin(const std::vector<double> &cdf);
153 :
154 : /// Random point on a unit-sphere
155 : Vector3d randVector();
156 : /// Random vector with given angular separation around mean direction
157 : Vector3d randVectorAroundMean(const Vector3d &meanDirection, double angle);
158 : /// Fisher distributed random vector
159 : Vector3d randFisherVector(const Vector3d &meanDirection, double kappa);
160 : /// Uniform distributed random vector inside a cone
161 : Vector3d randConeVector(const Vector3d &meanDirection, double angularRadius);
162 : /// Random lamberts distributed vector with theta distribution: sin(t) * cos(t),
163 : /// aka cosine law (https://en.wikipedia.org/wiki/Lambert%27s_cosine_law),
164 : /// for a surface element with normal vector pointing in positive z-axis (0, 0, 1)
165 : Vector3d randVectorLamberts();
166 : /// Same as above but rotated to the respective normalVector of surface element
167 : Vector3d randVectorLamberts(const Vector3d &normalVector);
168 : ///_Position vector uniformly distributed within propagation step size bin
169 : Vector3d randomInterpolatedPosition(const Vector3d &a, const Vector3d &b);
170 :
171 : /// Power-law distribution of a given differential spectral index
172 : double randPowerLaw(double index, double min, double max);
173 : /// Broken power-law distribution
174 : double randBrokenPowerLaw(double index1, double index2, double breakpoint, double min, double max );
175 :
176 : /// Seed the generator with a simple uint32_t
177 : void seed( const uint32_t oneSeed );
178 : /// Seed the generator with an array of uint32_t's
179 : /// There are 2^19937-1 possible initial states. This function allows
180 : /// all of those to be accessed by providing at least 19937 bits (with a
181 : /// default seed length of N = 624 uint32_t's). Any bits above the lower 32
182 : /// in each element are discarded.
183 : /// Just call seed() if you want to get array from /dev/urandom
184 : void seed( uint32_t *const bigSeed, const uint32_t seedLength = N );
185 : // seed via an b64 encoded string
186 : void seed( const std::string &b64Seed);
187 : /// Seed the generator with an array from /dev/urandom if available
188 : /// Otherwise use a hash of time() and clock() values
189 : void seed();
190 :
191 : // Saving and loading generator state
192 : void save( uint32_t* saveArray ) const;// to array of size SAVE
193 : void load( uint32_t *const loadArray );// from such array
194 : const std::vector<uint32_t> &getSeed() const; // copy the seed to the array
195 : const std::string getSeed_base64() const; // get the base 64 encoded seed
196 :
197 : friend std::ostream& operator<<( std::ostream& os, const Random& mtrand );
198 : friend std::istream& operator>>( std::istream& is, Random& mtrand );
199 :
200 : static Random &instance();
201 : static void seedThreads(const uint32_t oneSeed);
202 : static std::vector< std::vector<uint32_t> > getSeedThreads();
203 :
204 : protected:
205 : /// Initialize generator state with seed
206 : /// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
207 : /// In previous versions, most significant bits (MSBs) of the seed affect
208 : /// only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
209 : void initialize( const uint32_t oneSeed );
210 :
211 : /// Generate N new values in state
212 : /// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
213 : void reload();
214 28746432 : uint32_t hiBit( const uint32_t& u ) const {return u & 0x80000000UL;}
215 28746432 : uint32_t loBit( const uint32_t& u ) const {return u & 0x00000001UL;}
216 28746432 : uint32_t loBits( const uint32_t& u ) const {return u & 0x7fffffffUL;}
217 : uint32_t mixBits( const uint32_t& u, const uint32_t& v ) const
218 28746432 : { return hiBit(u) | loBits(v);}
219 :
220 : #ifdef _MSC_VER
221 : #pragma warning( push )
222 : #pragma warning( disable : 4146 )
223 : #endif
224 : uint32_t twist( const uint32_t& m, const uint32_t& s0, const uint32_t& s1 ) const
225 28700364 : { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL);}
226 :
227 : #ifdef _MSC_VER
228 : #pragma warning( pop )
229 : #endif
230 :
231 : /// Get a uint32_t from t and c
232 : /// Better than uint32_t(x) in case x is floating point in [0,1]
233 : /// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
234 : static uint32_t hash( time_t t, clock_t c );
235 :
236 : };
237 : /** @}*/
238 :
239 : } //namespace crpropa
240 :
241 : #endif // RANDOM_H
|
