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 typedef unsigned long int myID_t; 50 using myuint_t = std::uint64_t; << 47 typedef unsigned long long int myuint_t; 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 #endif >> 148 inline double convert1double(myuint_t u) 185 { 149 { 186 return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) << 150 const double one = 1; >> 151 const myuint_t onemask = *(myuint_t*)&one; >> 152 myuint_t tmp = (u>>9) | onemask; // bits between 52 and 62 dont affect the result! >> 153 double d = *(double*)&tmp; >> 154 return d-1.0; 187 } 155 } >> 156 #if defined __GNUC__ >> 157 #pragma GCC diagnostic pop >> 158 #endif >> 159 myuint_t MOD_MULSPEC(myuint_t k); >> 160 myuint_t MULWU(myuint_t k); >> 161 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 162 myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld); 189 myuint_t apply_bigskip(myuint_t* Vout, myuin 163 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 << 164 myuint_t modadd(myuint_t foo, myuint_t bar); 191 { << 192 return MIXMAX_MOD_MERSENNE(xfoo+xbar); << 193 } << 194 << 195 #if defined(__x86_64__) 165 #if defined(__x86_64__) 196 myuint_t mod128(__uint128_t s); 166 myuint_t mod128(__uint128_t s); 197 myuint_t fmodmulM61(myuint_t cum, myuint_t a 167 myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b); 198 #else // on all other platforms, including 32- 168 #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 169 myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a); 200 #endif 170 #endif 201 171 202 // Engine state << 172 private: >> 173 >> 174 struct rng_state_st >> 175 { >> 176 std::array<myuint_t, N> V; >> 177 myuint_t sumtot; >> 178 int counter; >> 179 }; 203 180 204 myuint_t V[N] = {0}; << 181 typedef struct rng_state_st rng_state_t; // struct alias 205 myuint_t sumtot = 0; << 182 rng_state_t S; 206 int counter = N; << 207 }; 183 }; 208 184 209 } // namespace CLHEP 185 } // namespace CLHEP 210 186 211 #endif 187 #endif 212 188