| Geant4 Cross Reference |
1 //
1 //
2 // -*- C++ -*-
2 // -*- C++ -*-
3 //
3 //
4 // ------------------------------------------- 4 // -----------------------------------------------------------------------
5 // HEP Random
5 // HEP Random
6 // --- MixMaxRng ---
6 // --- MixMaxRng ---
7 // class header file
7 // class header file
8 // ------------------------------------------- 8 // -----------------------------------------------------------------------
9 //
9 //
10 // This file interfaces the MixMax PseudoRando 10 // This file interfaces the MixMax PseudoRandom Number Generator
11 // proposed by:
11 // proposed by:
12 //
12 //
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
14 // On the Monte Carlo simulation of physical 14 // On the Monte Carlo simulation of physical systems,
15 // J.Comput.Phys. 97, 566 (1991);
15 // J.Comput.Phys. 97, 566 (1991);
16 // Preprint EPI-865-16-86, Yerevan, Jan. 198 16 // Preprint EPI-865-16-86, Yerevan, Jan. 1986
17 // http://dx.doi.org/10.1016/0021-9991(91)90 17 // http://dx.doi.org/10.1016/0021-9991(91)90015-D
18 //
18 //
19 // K.Savvidy
19 // K.Savvidy
20 // "The MIXMAX random number generator"
20 // "The MIXMAX random number generator"
21 // Comp. Phys. Commun. (2015)
21 // Comp. Phys. Commun. (2015)
22 // http://dx.doi.org/10.1016/j.cpc.2015.06.0 22 // http://dx.doi.org/10.1016/j.cpc.2015.06.003
23 //
23 //
24 // K.Savvidy and G.Savvidy
24 // K.Savvidy and G.Savvidy
25 // "Spectrum and Entropy of C-systems. MIXMA 25 // "Spectrum and Entropy of C-systems. MIXMAX random number generator"
26 // Chaos, Solitons & Fractals, Volume 91, (2 26 // Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33-38
27 // http://dx.doi.org/10.1016/j.chaos.2016.05 27 // http://dx.doi.org/10.1016/j.chaos.2016.05.003
28 //
28 //
29 // =========================================== 29 // =======================================================================
30 // Implementation by Konstantin Savvidy - Copy << 30 // Implementation by Konstantin Savvidy - Copyright 2004-2017
31 // July 2023 - Updated class structure upon su <<
32 // September 2023 - fix (re-)initialization fr <<
33 // =========================================== 31 // =======================================================================
34
32
35 #ifndef MixMaxRng_h
33 #ifndef MixMaxRng_h
36 #define MixMaxRng_h 1
34 #define MixMaxRng_h 1
37
35
38 #include <array>
36 #include <array>
39 #include <cstdint>
<<
40 #include "CLHEP/Random/RandomEngine.h"
37 #include "CLHEP/Random/RandomEngine.h"
41
38
42 namespace CLHEP {
39 namespace CLHEP {
43
40
44 /**
41 /**
45 * @author K.Savvidy
42 * @author K.Savvidy
46 * @ingroup random
43 * @ingroup random
47 */
44 */
48
45
49 using myID_t = std::uint32_t;
46 using myID_t = std::uint32_t;
50 using myuint_t = std::uint64_t;
<< 47 using myuint_t = unsigned long long int;
51
48
52 class alignas(128) MixMaxRng : public HepRando << 49 class MixMaxRng: public HepRandomEngine {
53 {
<<
54
50
55 static const int N = 17;
51 static const int N = 17;
56
52
57 public:
53 public:
58
54
59 MixMaxRng(std::istream& is);
55 MixMaxRng(std::istream& is);
60 MixMaxRng();
56 MixMaxRng();
61 MixMaxRng(long seed);
57 MixMaxRng(long seed);
62 ~MixMaxRng();
58 ~MixMaxRng();
63 // Constructors and destructor.
59 // Constructors and destructor.
64
60
65 MixMaxRng(const MixMaxRng& rng);
61 MixMaxRng(const MixMaxRng& rng);
66 MixMaxRng& operator=(const MixMaxRng& rng);
62 MixMaxRng& operator=(const MixMaxRng& rng);
67 // Copy constructor and assignment operator. 63 // Copy constructor and assignment operator.
68
64
69 inline double flat()
<< 65 double flat() { return (S.counter<=(N-1)) ? generate(S.counter):iterate(); }
70 {
<<
71 if (counter >= N) iterate();
<<
72 return INV_M61*static_cast<double>(V[count <<
73 }
<<
74 // Returns a pseudo random number between 0 66 // Returns a pseudo random number between 0 and 1
75 // excluding the zero: in (0,1]
<< 67 // (excluding the zero: in (0,1] )
76 // smallest number which it will give is app 68 // smallest number which it will give is approximately 10^-19
77
69
78 void flatArray (const int size, double* vect 70 void flatArray (const int size, double* vect);
79 // Fills the array "vect" of specified size 71 // Fills the array "vect" of specified size with flat random values.
80
72
81 inline void setSeed(long longSeed, int = 0 / << 73 void setSeed(long seed, int dum=0);
82 {
<<
83 seed_spbox(theSeed = longSeed);
<<
84 }
<<
85 // Sets the state of the algorithm according 74 // Sets the state of the algorithm according to seed.
86
75
87 void setSeeds(const long * seeds, int seedNu 76 void setSeeds(const long * seeds, int seedNum=0);
88 // Sets the initial state of the engine acco 77 // Sets the initial state of the engine according to the array of between one and four 32-bit seeds.
89 // If the size of long is greater on the pla 78 // If the size of long is greater on the platform, only the lower 32-bits are used.
90 // Streams created from seeds differing by a 79 // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
91 // to be independent and non-colliding for a 80 // to be independent and non-colliding for at least the next 10^100 random numbers
92
81
93 void saveStatus( const char filename[] = "Mi 82 void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
94 // Saves the the current engine state in the 83 // Saves the the current engine state in the file given, by default MixMaxRngState.conf
95
84
96 void restoreStatus( const char filename[] = 85 void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
97 // Reads a valid engine state from a given f 86 // Reads a valid engine state from a given file, by default MixMaxRngState.conf
98 // and restores it.
87 // and restores it.
99
88
100 void showStatus() const;
89 void showStatus() const;
101 // Dumps the engine status on the screen.
90 // Dumps the engine status on the screen.
102
91
103 inline operator double() { return flat(); }
<< 92 operator double();
104 // Returns same as flat()
93 // Returns same as flat()
105
<< 94 operator float();
106 inline operator float() { return float( flat <<
107 // less precise flat, faster if possible
95 // less precise flat, faster if possible
108
<< 96 operator unsigned int();
109 inline operator unsigned int() { return stat << 97 // 32-bit flat
110 // 32-bit flat. clhep_get_next() returns a 6 <<
111 // the lower 61 bits are random and upper 3 <<
112
98
113 virtual std::ostream & put (std::ostream & o 99 virtual std::ostream & put (std::ostream & os) const;
114 virtual std::istream & get (std::istream & i 100 virtual std::istream & get (std::istream & is);
115 static std::string beginTag ( );
101 static std::string beginTag ( );
116 virtual std::istream & getState ( std::istre 102 virtual std::istream & getState ( std::istream & is );
117
103
118 std::string name() const { return "MixMaxRng 104 std::string name() const { return "MixMaxRng"; }
119 static std::string engineName();
105 static std::string engineName();
120
106
121 std::vector<unsigned long> put () const;
107 std::vector<unsigned long> put () const;
122 bool get (const std::vector<unsigned long> & << 108 bool get (const std::vector<unsigned long> & v);
123 bool getState (const std::vector<unsigned lo << 109 bool getState (const std::vector<unsigned long> & v);
124
110
125 private:
111 private:
126
112
127 static constexpr long long int SPECIAL = 0 << 113 static constexpr long long int SPECIAL = ((N==17)? 0 : ((N==240)? 487013230256099140ULL:0) ); // etc...
128 static constexpr long long int SPECIALMUL= 3 << 114 static constexpr long long int SPECIALMUL= ((N==17)? 36: ((N==240)? 51 :53) ); // etc...
>> 115 // Note the potential for confusion...
129 static constexpr int BITS=61;
116 static constexpr int BITS=61;
130 static constexpr myuint_t M61=23058430092136 117 static constexpr myuint_t M61=2305843009213693951ULL;
131 static constexpr double INV_M61=0.4336808689 118 static constexpr double INV_M61=0.43368086899420177360298E-18;
132 static constexpr unsigned int VECTOR_STATE_S 119 static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX
133
120
134 inline myuint_t MIXMAX_MOD_MERSENNE(myuint_t << 121 #define MIXMAX_MOD_MERSENNE(k) ((((k)) & M61) + (((k)) >> BITS) )
135 {
<<
136 return ((((k)) & M61) + (((k)) >> BITS) ); <<
137 }
<<
138
122
139 static constexpr int rng_get_N();
123 static constexpr int rng_get_N();
>> 124 static constexpr long long int rng_get_SPECIAL();
>> 125 static constexpr int rng_get_SPECIALMUL();
140 void seed_uniquestream( myID_t clusterID, my 126 void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
141 void seed_spbox(myuint_t seed);
<< 127 void seed_spbox(myuint_t);
142 void print_state() const;
128 void print_state() const;
143 myuint_t precalc();
129 myuint_t precalc();
144 myuint_t get_next();
<< 130 myuint_t get_next() ;
>> 131 inline double get_next_float() { return get_next_float_packbits(); }
>> 132 // Returns a random double with all 52 bits random, in the range (0,1]
145
133
146 MixMaxRng Branch();
134 MixMaxRng Branch();
147 void BranchInplace(int id);
135 void BranchInplace(int id);
148
136
149 MixMaxRng(myID_t clusterID, myID_t machineID 137 MixMaxRng(myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ); // Constructor with four 32-bit seeds
150 inline void seed64(myuint_t seedval) // se << 138 inline void seed64(myuint_t seedval) { seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (myID_t)seedval ); } // seed with one 64-bit seed
151 {
<<
152 seed_uniquestream( 0, 0, (myID_t)(seedval> <<
153 }
<<
154
139
155 inline void iterate()
<< 140 double generate(int i);
156 {
<< 141 double iterate();
157 myuint_t tempP, tempV;
<<
158 V[0] = ( tempV = sumtot );
<<
159 myuint_t insumtot = V[0], ovflow = 0; // w <<
160 tempP = 0; // w <<
161 myuint_t tempPO;
<<
162 tempPO = MULWU(tempP); tempP = modadd(temp <<
163 tempPO = MULWU(tempP); tempP = modadd(temp <<
164 tempPO = MULWU(tempP); tempP = modadd(temp <<
165 tempPO = MULWU(tempP); tempP = modadd(temp <<
166 tempPO = MULWU(tempP); tempP = modadd(temp <<
167 tempPO = MULWU(tempP); tempP = modadd(temp <<
168 tempPO = MULWU(tempP); tempP = modadd(temp <<
169 tempPO = MULWU(tempP); tempP = modadd(temp <<
170 tempPO = MULWU(tempP); tempP = modadd(temp <<
171 tempPO = MULWU(tempP); tempP = modadd(temp <<
172 tempPO = MULWU(tempP); tempP = modadd(temp <<
173 tempPO = MULWU(tempP); tempP = modadd(temp <<
174 tempPO = MULWU(tempP); tempP = modadd(temp <<
175 tempPO = MULWU(tempP); tempP = modadd(temp <<
176 tempPO = MULWU(tempP); tempP = modadd(temp <<
177 tempPO = MULWU(tempP); tempP = modadd(temp <<
178 sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_ME <<
179
142
180 counter=1;
<< 143 double get_next_float_packbits();
181 }
<< 144 #if defined __GNUC__
182
<< 145 #pragma GCC diagnostic push
183 void state_init();
<< 146 #pragma GCC diagnostic ignored "-Wstrict-aliasing"
184 inline myuint_t MULWU (myuint_t k)
<< 147 #pragma GCC diagnostic ignored "-Wuninitialized"
>> 148 #endif
>> 149 inline double convert1double(myuint_t u)
185 {
150 {
186 return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) << 151 const double one = 1;
>> 152 const myuint_t onemask = *(myuint_t*)&one;
>> 153 myuint_t tmp = (u>>9) | onemask; // bits between 52 and 62 dont affect the result!
>> 154 double d = *(double*)&tmp;
>> 155 return d-1.0;
187 }
156 }
>> 157 #if defined __GNUC__
>> 158 #pragma GCC diagnostic pop
>> 159 #endif
>> 160 myuint_t MOD_MULSPEC(myuint_t k);
>> 161 myuint_t MULWU(myuint_t k);
>> 162 void seed_vielbein( unsigned int i); // seeds with the i-th unit vector, i = 0..N-1, for testing only
188 myuint_t iterate_raw_vec(myuint_t* Y, myuint 163 myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
189 myuint_t apply_bigskip(myuint_t* Vout, myuin 164 myuint_t apply_bigskip(myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
190 inline myuint_t modadd(myuint_t xfoo, myuint << 165 myuint_t modadd(myuint_t foo, myuint_t bar);
191 {
<<
192 return MIXMAX_MOD_MERSENNE(xfoo+xbar);
<<
193 }
<<
194
<<
195 #if defined(__x86_64__)
166 #if defined(__x86_64__)
196 myuint_t mod128(__uint128_t s);
167 myuint_t mod128(__uint128_t s);
197 myuint_t fmodmulM61(myuint_t cum, myuint_t a 168 myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b);
198 #else // on all other platforms, including 32- 169 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
199 myuint_t fmodmulM61(myuint_t cum, myuint_t s 170 myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a);
200 #endif
171 #endif
201
172
202 // Engine state
<< 173 private:
>> 174
>> 175 struct rng_state_st
>> 176 {
>> 177 std::array<myuint_t, N> V;
>> 178 myuint_t sumtot;
>> 179 int counter;
>> 180 };
203
181
204 myuint_t V[N] = {0};
<< 182 typedef struct rng_state_st rng_state_t; // struct alias
205 myuint_t sumtot = 0;
<< 183 rng_state_t S;
206 int counter = N;
<<
207 };
184 };
208
185
209 } // namespace CLHEP
186 } // namespace CLHEP
210
187
211 #endif
188 #endif
212 189