CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
MixMaxRng.h
Go to the documentation of this file.
1//
2// -*- C++ -*-
3//
4// -----------------------------------------------------------------------
5// HEP Random
6// --- MixMaxRng ---
7// class header file
8// -----------------------------------------------------------------------
9//
10// This file interfaces the MixMax PseudoRandom Number Generator
11// proposed by:
12//
13// G.K.Savvidy and N.G.Ter-Arutyunian,
14// On the Monte Carlo simulation of physical systems,
15// J.Comput.Phys. 97, 566 (1991);
16// Preprint EPI-865-16-86, Yerevan, Jan. 1986
17// http://dx.doi.org/10.1016/0021-9991(91)90015-D
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.003
23//
24// K.Savvidy and G.Savvidy
25// "Spectrum and Entropy of C-systems. MIXMAX random number generator"
26// Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33-38
27// http://dx.doi.org/10.1016/j.chaos.2016.05.003
28//
29// =======================================================================
30// Implementation by Konstantin Savvidy - Copyright 2004-2023
31// July 2023 - Updated class structure upon suggestions from Marco Barbone
32// September 2023 - fix (re-)initialization from Gabriele Cosmo
33// =======================================================================
34
35#ifndef MixMaxRng_h
36#define MixMaxRng_h 1
37
38#include <array>
39#include <cstdint>
40#include "CLHEP/Random/defs.h"
42
43namespace CLHEP {
44
49
50using myID_t = std::uint32_t;
51using myuint_t = std::uint64_t;
52
53class
54#if (__cplusplus >= 201703L && __cpp_aligned_new >= 201606L)
55 alignas(128)
56#endif
58{
59
60 static const int N = 17;
61
62public:
63
64 MixMaxRng(std::istream& is);
66 MixMaxRng(long seed);
68 // Constructors and destructor.
69
70 MixMaxRng(const MixMaxRng& rng);
72 // Copy constructor and assignment operator.
73
74 inline double flat()
75 {
76 if (counter >= N) iterate();
77 return INV_M61*static_cast<double>(V[counter++]);
78 }
79 // Returns a pseudo random number between 0 and 1
80 // excluding the zero: in (0,1]
81 // smallest number which it will give is approximately 10^-19
82
83 void flatArray (const int size, double* vect);
84 // Fills the array "vect" of specified size with flat random values.
85
86 inline void setSeed(long longSeed, int = 0 /* extraSeed */)
87 {
88 seed_spbox(theSeed = longSeed);
89 }
90 // Sets the state of the algorithm according to seed.
91
92 void setSeeds(const long * seeds, int seedNum=0);
93 // Sets the initial state of the engine according to the array of between one and four 32-bit seeds.
94 // If the size of long is greater on the platform, only the lower 32-bits are used.
95 // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
96 // to be independent and non-colliding for at least the next 10^100 random numbers
97
98 void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
99 // Saves the the current engine state in the file given, by default MixMaxRngState.conf
100
101 void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
102 // Reads a valid engine state from a given file, by default MixMaxRngState.conf
103 // and restores it.
104
105 void showStatus() const;
106 // Dumps the engine status on the screen.
107
108 inline operator double() { return flat(); }
109 // Returns same as flat()
110
111 inline operator float() { return float( flat() ); }
112 // less precise flat, faster if possible
113
114 inline operator unsigned int() { return static_cast<unsigned int>(get_next()); }
115 // 32-bit flat. clhep_get_next() returns a 64-bit integer, of which
116 // the lower 61 bits are random and upper 3 bits are zero
117
118 virtual std::ostream & put (std::ostream & os) const;
119 virtual std::istream & get (std::istream & is);
120 static std::string beginTag ( );
121 virtual std::istream & getState ( std::istream & is );
122
123 std::string name() const { return "MixMaxRng"; }
124 static std::string engineName();
125
126 std::vector<unsigned long> put () const;
127 bool get (const std::vector<unsigned long> & vec);
128 bool getState (const std::vector<unsigned long> & vec);
129
130private:
131
132 static constexpr long long int SPECIAL = 0;
133 static constexpr long long int SPECIALMUL= 36;
134 static constexpr int BITS=61;
135 static constexpr myuint_t M61=2305843009213693951ULL;
136 static constexpr double INV_M61=0.43368086899420177360298E-18;
137 static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX
138
139 inline myuint_t MIXMAX_MOD_MERSENNE(myuint_t k)
140 {
141 return ((((k)) & M61) + (((k)) >> BITS) );
142 }
143
144 static constexpr int rng_get_N();
145 void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
146 void seed_spbox(myuint_t seed);
147 void print_state() const;
148 myuint_t precalc();
149 myuint_t get_next();
150
151 MixMaxRng Branch();
152 void BranchInplace(int id);
153
154 MixMaxRng(myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ); // Constructor with four 32-bit seeds
155 inline void seed64(myuint_t seedval) // seed with one 64-bit seed
156 {
157 seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (myID_t)seedval );
158 }
159
160 inline void iterate()
161 {
162 myuint_t tempP, tempV;
163 V[0] = ( tempV = sumtot );
164 myuint_t insumtot = V[0], ovflow = 0; // will keep a running sum of all new elements
165 tempP = 0; // will keep a partial sum of all old elements
166 myuint_t tempPO;
167 tempPO = MULWU(tempP); tempP = modadd(tempP, V[1] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[1] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
168 tempPO = MULWU(tempP); tempP = modadd(tempP, V[2] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[2] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
169 tempPO = MULWU(tempP); tempP = modadd(tempP, V[3] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[3] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
170 tempPO = MULWU(tempP); tempP = modadd(tempP, V[4] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[4] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
171 tempPO = MULWU(tempP); tempP = modadd(tempP, V[5] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[5] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
172 tempPO = MULWU(tempP); tempP = modadd(tempP, V[6] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[6] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
173 tempPO = MULWU(tempP); tempP = modadd(tempP, V[7] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[7] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
174 tempPO = MULWU(tempP); tempP = modadd(tempP, V[8] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[8] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
175 tempPO = MULWU(tempP); tempP = modadd(tempP, V[9] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[9] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
176 tempPO = MULWU(tempP); tempP = modadd(tempP, V[10]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[10] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
177 tempPO = MULWU(tempP); tempP = modadd(tempP, V[11]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[11] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
178 tempPO = MULWU(tempP); tempP = modadd(tempP, V[12]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[12] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
179 tempPO = MULWU(tempP); tempP = modadd(tempP, V[13]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[13] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
180 tempPO = MULWU(tempP); tempP = modadd(tempP, V[14]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[14] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
181 tempPO = MULWU(tempP); tempP = modadd(tempP, V[15]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[15] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
182 tempPO = MULWU(tempP); tempP = modadd(tempP, V[16]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[16] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
183 sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(insumtot) + (ovflow <<3 ));
184
185 counter=1;
186 }
187
188 void state_init();
189 inline myuint_t MULWU (myuint_t k)
190 {
191 return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) >> (BITS-SPECIALMUL)) );
192 }
193 myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
194 myuint_t apply_bigskip(myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
195 inline myuint_t modadd(myuint_t xfoo, myuint_t xbar)
196 {
197 return MIXMAX_MOD_MERSENNE(xfoo+xbar);
198 }
199
200#if defined(__x86_64__)
201 myuint_t mod128(__uint128_t s);
202 myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b);
203#else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
204 myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a);
205#endif
206
207 // Engine state
208
209 myuint_t V[N] = {0};
210 myuint_t sumtot = 0;
211 int counter = N;
212};
213
214} // namespace CLHEP
215
216#endif
void restoreStatus(const char filename[]="MixMaxRngState.conf")
void showStatus() const
static std::string engineName()
void flatArray(const int size, double *vect)
MixMaxRng(const MixMaxRng &rng)
std::string name() const
Definition MixMaxRng.h:123
void saveStatus(const char filename[]="MixMaxRngState.conf") const
virtual std::istream & get(std::istream &is)
MixMaxRng & operator=(const MixMaxRng &rng)
std::vector< unsigned long > put() const
bool getState(const std::vector< unsigned long > &vec)
double flat()
Definition MixMaxRng.h:74
MixMaxRng(std::istream &is)
bool get(const std::vector< unsigned long > &vec)
virtual std::istream & getState(std::istream &is)
void setSeeds(const long *seeds, int seedNum=0)
virtual std::ostream & put(std::ostream &os) const
void setSeed(long longSeed, int=0)
Definition MixMaxRng.h:86
static std::string beginTag()
MixMaxRng(long seed)
std::uint64_t myuint_t
Definition MixMaxRng.h:51
std::uint32_t myID_t
Definition MixMaxRng.h:50