| Geant4 Cross Reference |
1 //
1
2 // -*- C++ -*-
3 //
4 // -------------------------------------------
5 // HEP Random
6 // --- MixMaxRng ---
7 // class header file
8 // -------------------------------------------
9 //
10 // This file interfaces the MixMax PseudoRando
11 // proposed by:
12 //
13 // G.K.Savvidy and N.G.Ter-Arutyunian,
14 // On the Monte Carlo simulation of physical
15 // J.Comput.Phys. 97, 566 (1991);
16 // Preprint EPI-865-16-86, Yerevan, Jan. 198
17 // http://dx.doi.org/10.1016/0021-9991(91)90
18 //
19 // K.Savvidy
20 // "The MIXMAX random number generator"
21 // Comp. Phys. Commun. (2015)
22 // http://dx.doi.org/10.1016/j.cpc.2015.06.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 // ===========================================
30 // Implementation by Konstantin Savvidy - Copy
31 // July 2023 - Updated class structure upon su
32 // September 2023 - fix (re-)initialization fr
33 // ===========================================
34
35 #ifndef MixMaxRng_h
36 #define MixMaxRng_h 1
37
38 #include <array>
39 #include <cstdint>
40 #include "CLHEP/Random/RandomEngine.h"
41
42 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
52 class alignas(128) MixMaxRng : public HepRando
53 {
54
55 static const int N = 17;
56
57 public:
58
59 MixMaxRng(std::istream& is);
60 MixMaxRng();
61 MixMaxRng(long seed);
62 ~MixMaxRng();
63 // Constructors and destructor.
64
65 MixMaxRng(const MixMaxRng& rng);
66 MixMaxRng& operator=(const MixMaxRng& rng);
67 // Copy constructor and assignment operator.
68
69 inline 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
75 // excluding the zero: in (0,1]
76 // smallest number which it will give is app
77
78 void flatArray (const int size, double* vect
79 // Fills the array "vect" of specified size
80
81 inline void setSeed(long longSeed, int = 0 /
82 {
83 seed_spbox(theSeed = longSeed);
84 }
85 // Sets the state of the algorithm according
86
87 void setSeeds(const long * seeds, int seedNu
88 // Sets the initial state of the engine acco
89 // If the size of long is greater on the pla
90 // Streams created from seeds differing by a
91 // to be independent and non-colliding for a
92
93 void saveStatus( const char filename[] = "Mi
94 // Saves the the current engine state in the
95
96 void restoreStatus( const char filename[] =
97 // Reads a valid engine state from a given f
98 // and restores it.
99
100 void showStatus() const;
101 // Dumps the engine status on the screen.
102
103 inline operator double() { return flat(); }
104 // Returns same as 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
113 virtual std::ostream & put (std::ostream & o
114 virtual std::istream & get (std::istream & i
115 static std::string beginTag ( );
116 virtual std::istream & getState ( std::istre
117
118 std::string name() const { return "MixMaxRng
119 static std::string engineName();
120
121 std::vector<unsigned long> put () const;
122 bool get (const std::vector<unsigned long> &
123 bool getState (const std::vector<unsigned lo
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
146 MixMaxRng Branch();
147 void BranchInplace(int id);
148
149 MixMaxRng(myID_t clusterID, myID_t machineID
150 inline void seed64(myuint_t seedval) // se
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 };
208
209 } // namespace CLHEP
210
211 #endif
212