| Geant4 Cross Reference |
1 //
<< 1 // $Id:$
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 PseudoRandom Number Generator
11 // proposed by:
11 // proposed by:
12 //
<< 12 // N.Z. Akopov, G.K.Saviddy & N.G. Ter-Arutyunian
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
<< 13 // "Matrix Generator of Pseudorandom Numbers",
14 // On the Monte Carlo simulation of physical << 14 // J.Compt.Phy. 97, 573 (1991)
15 // J.Comput.Phys. 97, 566 (1991);
<< 15 // Preprint: EPI-867(18)-86, Yerevan June 1986.
16 // Preprint EPI-865-16-86, Yerevan, Jan. 198 << 16 // G. Savvidy & N. Savvidy
17 // http://dx.doi.org/10.1016/0021-9991(91)90 << 17 // "On the Monte Carlo Simulation of Physical Systems",
18 //
<< 18 // J.Comput.Phys. 97 (1991) 566
19 // K.Savvidy
<< 19
20 // "The MIXMAX random number generator"
<< 20 // =======================================================================
21 // Comp. Phys. Commun. (2015)
<< 21 // Implementation by Konstantin Savvidy - 2004-2015
22 // http://dx.doi.org/10.1016/j.cpc.2015.06.0 << 22 // Release 0.99 and later: released under the LGPL license version 3.0
23 //
<<
24 // K.Savvidy and G.Savvidy
<<
25 // "Spectrum and Entropy of C-systems. MIXMA <<
26 // Chaos, Solitons & Fractals, Volume 91, (2 <<
27 // http://dx.doi.org/10.1016/j.chaos.2016.05 <<
28 //
<<
29 // =========================================== 23 // =======================================================================
30 // Implementation by Konstantin Savvidy - Copy << 24 // CLHEP interface implemented by
31 // July 2023 - Updated class structure upon su << 25 // J. Apostolakis, G. Cosmo & K. Savvidy - Created: 6th July 2015
32 // September 2023 - fix (re-)initialization fr << 26 // CLHEP interface released under the LGPL license version 3.0
33 // =========================================== 27 // =======================================================================
34
28
35 #ifndef MixMaxRng_h
29 #ifndef MixMaxRng_h
36 #define MixMaxRng_h 1
30 #define MixMaxRng_h 1
37
31
38 #include <array>
<<
39 #include <cstdint>
<<
40 #include "CLHEP/Random/RandomEngine.h"
32 #include "CLHEP/Random/RandomEngine.h"
>> 33 #include "CLHEP/Random/mixmax.h"
41
34
42 namespace CLHEP {
35 namespace CLHEP {
43
<<
44 /**
<<
45 * @author K.Savvidy
<<
46 * @ingroup random
<<
47 */
<<
48
<<
49 using myID_t = std::uint32_t;
<<
50 using myuint_t = std::uint64_t;
<<
51
36
52 class alignas(128) MixMaxRng : public HepRando << 37 /**
53 {
<< 38 * @author K. Savvidy
54
<< 39 * @ingroup random
55 static const int N = 17;
<< 40 */
>> 41 class MixMaxRng: public HepRandomEngine {
56
42
57 public:
43 public:
58
44
59 MixMaxRng(std::istream& is);
45 MixMaxRng(std::istream& is);
60 MixMaxRng();
46 MixMaxRng();
61 MixMaxRng(long seed);
47 MixMaxRng(long seed);
62 ~MixMaxRng();
<< 48 MixMaxRng(int rowIndex, int colIndex);
63 // Constructors and destructor.
<< 49 virtual ~MixMaxRng();
>> 50 // Constructor and destructor.
64
51
65 MixMaxRng(const MixMaxRng& rng);
52 MixMaxRng(const MixMaxRng& rng);
66 MixMaxRng& operator=(const MixMaxRng& rng);
53 MixMaxRng& operator=(const MixMaxRng& rng);
67 // Copy constructor and assignment operator. 54 // Copy constructor and assignment operator.
68
55
69 inline double flat()
<< 56 double flat();
70 {
<<
71 if (counter >= N) iterate();
<<
72 return INV_M61*static_cast<double>(V[count <<
73 }
<<
74 // Returns a pseudo random number between 0 57 // Returns a pseudo random number between 0 and 1
75 // excluding the zero: in (0,1]
<< 58 // (excluding the zero: in (0,1] )
76 // smallest number which it will give is app 59 // smallest number which it will give is approximately 10^-19
77
60
78 void flatArray (const int size, double* vect 61 void flatArray (const int size, double* vect);
79 // Fills the array "vect" of specified size 62 // Fills the array "vect" of specified size with flat random values.
80
63
81 inline void setSeed(long longSeed, int = 0 / << 64 void setSeed(long seed, int dum=0);
82 {
<<
83 seed_spbox(theSeed = longSeed);
<<
84 }
<<
85 // Sets the state of the algorithm according 65 // Sets the state of the algorithm according to seed.
86
66
87 void setSeeds(const long * seeds, int seedNu 67 void setSeeds(const long * seeds, int seedNum=0);
88 // Sets the initial state of the engine acco 68 // 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 69 // 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 70 // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
91 // to be independent and non-colliding for a 71 // to be independent and non-colliding for at least the next 10^100 random numbers
92
72
93 void saveStatus( const char filename[] = "Mi 73 void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
94 // Saves the the current engine state in the 74 // Saves the the current engine state in the file given, by default MixMaxRngState.conf
95
75
96 void restoreStatus( const char filename[] = 76 void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
97 // Reads a valid engine state from a given f 77 // Reads a valid engine state from a given file, by default MixMaxRngState.conf
98 // and restores it.
78 // and restores it.
99
79
100 void showStatus() const;
80 void showStatus() const;
101 // Dumps the engine status on the screen.
81 // Dumps the engine status on the screen.
102
82
103 inline operator double() { return flat(); }
<< 83 operator unsigned int();
104 // Returns same as flat()
<< 84 // 32-bit flat
105
<<
106 inline operator float() { return float( flat <<
107 // less precise flat, faster if possible
<<
108
<<
109 inline operator unsigned int() { return stat <<
110 // 32-bit flat. clhep_get_next() returns a 6 <<
111 // the lower 61 bits are random and upper 3 <<
112
85
113 virtual std::ostream & put (std::ostream & o 86 virtual std::ostream & put (std::ostream & os) const;
114 virtual std::istream & get (std::istream & i 87 virtual std::istream & get (std::istream & is);
115 static std::string beginTag ( );
88 static std::string beginTag ( );
116 virtual std::istream & getState ( std::istre 89 virtual std::istream & getState ( std::istream & is );
117
90
118 std::string name() const { return "MixMaxRng << 91 std::string name() const;
119 static std::string engineName();
<< 92 static std::string engineName() {return "MixMaxRng";}
120
93
121 std::vector<unsigned long> put () const;
94 std::vector<unsigned long> put () const;
122 bool get (const std::vector<unsigned long> & << 95 bool get (const std::vector<unsigned long> & v);
123 bool getState (const std::vector<unsigned lo << 96 bool getState (const std::vector<unsigned long> & v);
124
<<
125 private:
<<
126
<<
127 static constexpr long long int SPECIAL = 0 <<
128 static constexpr long long int SPECIALMUL= 3 <<
129 static constexpr int BITS=61;
<<
130 static constexpr myuint_t M61=23058430092136 <<
131 static constexpr double INV_M61=0.4336808689 <<
132 static constexpr unsigned int VECTOR_STATE_S <<
133
<<
134 inline myuint_t MIXMAX_MOD_MERSENNE(myuint_t <<
135 {
<<
136 return ((((k)) & M61) + (((k)) >> BITS) ); <<
137 }
<<
138
<<
139 static constexpr int rng_get_N();
<<
140 void seed_uniquestream( myID_t clusterID, my <<
141 void seed_spbox(myuint_t seed);
<<
142 void print_state() const;
<<
143 myuint_t precalc();
<<
144 myuint_t get_next();
<<
145
97
146 MixMaxRng Branch();
<< 98 static const unsigned int VECTOR_STATE_SIZE = 516; // 2N+4 for MIXMAX
147 void BranchInplace(int id);
<< 99
>> 100 private:
148
101
149 MixMaxRng(myID_t clusterID, myID_t machineID << 102 // Members defining the current status of the generator.
150 inline void seed64(myuint_t seedval) // se << 103 rng_state_st* fRngState;
151 {
<<
152 seed_uniquestream( 0, 0, (myID_t)(seedval> <<
153 }
<<
154
<<
155 inline void iterate()
<<
156 {
<<
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
<<
180 counter=1;
<<
181 }
<<
182
<<
183 void state_init();
<<
184 inline myuint_t MULWU (myuint_t k)
<<
185 {
<<
186 return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) <<
187 }
<<
188 myuint_t iterate_raw_vec(myuint_t* Y, myuint <<
189 myuint_t apply_bigskip(myuint_t* Vout, myuin <<
190 inline myuint_t modadd(myuint_t xfoo, myuint <<
191 {
<<
192 return MIXMAX_MOD_MERSENNE(xfoo+xbar);
<<
193 }
<<
194
<<
195 #if defined(__x86_64__)
<<
196 myuint_t mod128(__uint128_t s);
<<
197 myuint_t fmodmulM61(myuint_t cum, myuint_t a <<
198 #else // on all other platforms, including 32- <<
199 myuint_t fmodmulM61(myuint_t cum, myuint_t s <<
200 #endif
<<
201
<<
202 // Engine state
<<
203
<<
204 myuint_t V[N] = {0};
<<
205 myuint_t sumtot = 0;
<<
206 int counter = N;
<<
207 };
104 };
208
105
209 } // namespace CLHEP
106 } // namespace CLHEP
210
107
211 #endif
108 #endif
212 109