Geant4/externals/clhep/include/CLHEP/Random/MixMaxRng.h

Geant4 LXR

Geant4 Cross Reference

Cross-Referencing Geant4
Geant4/externals/clhep/include/CLHEP/Random/MixMaxRng.h

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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/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 HepRandomEngine
 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[counter++]);
 73   }
 74   // Returns a pseudo random number between 0 and 1
 75   // excluding the zero: in (0,1]
 76   // smallest number which it will give is approximately 10^-19
 77 
 78   void flatArray (const int size, double* vect);
 79   // Fills the array "vect" of specified size with flat random values.
 80 
 81   inline void setSeed(long longSeed, int = 0 /* extraSeed */)
 82   {
 83     seed_spbox(theSeed = longSeed);
 84   }
 85   // Sets the state of the algorithm according to seed.
 86 
 87   void setSeeds(const long * seeds, int seedNum=0);
 88   // 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 platform, only the lower 32-bits are used.
 90   // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
 91   // to be independent and non-colliding for at least the next 10^100 random numbers
 92 
 93   void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
 94   // Saves the the current engine state in the file given, by default MixMaxRngState.conf
 95 
 96   void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
 97   // Reads a valid engine state from a given file, by default MixMaxRngState.conf
 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 static_cast<unsigned int>(get_next()); }
110   // 32-bit flat. clhep_get_next() returns a 64-bit integer, of which
111   // the lower 61 bits are random and upper 3 bits are zero
112 
113   virtual std::ostream & put (std::ostream & os) const;
114   virtual std::istream & get (std::istream & is);
115   static  std::string beginTag ( );
116   virtual std::istream & getState ( std::istream & is );
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> & vec);
123   bool getState (const std::vector<unsigned long> & vec);
124 
125 private:
126 
127   static constexpr long long int SPECIAL   = 0;
128   static constexpr long long int SPECIALMUL= 36;
129   static constexpr int BITS=61;
130   static constexpr myuint_t M61=2305843009213693951ULL;
131   static constexpr double INV_M61=0.43368086899420177360298E-18;
132   static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX
133 
134   inline myuint_t MIXMAX_MOD_MERSENNE(myuint_t k)
135   {
136     return ((((k)) & M61) + (((k)) >> BITS) );
137   }
138 
139   static constexpr int rng_get_N();
140   void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
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, myID_t runID, myID_t  streamID );    // Constructor with four 32-bit seeds
150   inline void seed64(myuint_t seedval)   // seed with one 64-bit seed
151   {
152     seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (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; // will keep a running sum of all new elements
160     tempP = 0;                            // will keep a partial sum of all old elements
161     myuint_t tempPO;
162     tempPO = MULWU(tempP); tempP = modadd(tempP, V[1] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[1]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
163     tempPO = MULWU(tempP); tempP = modadd(tempP, V[2] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[2]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
164     tempPO = MULWU(tempP); tempP = modadd(tempP, V[3] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[3]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
165     tempPO = MULWU(tempP); tempP = modadd(tempP, V[4] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[4]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
166     tempPO = MULWU(tempP); tempP = modadd(tempP, V[5] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[5]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
167     tempPO = MULWU(tempP); tempP = modadd(tempP, V[6] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[6]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
168     tempPO = MULWU(tempP); tempP = modadd(tempP, V[7] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[7]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
169     tempPO = MULWU(tempP); tempP = modadd(tempP, V[8] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[8]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
170     tempPO = MULWU(tempP); tempP = modadd(tempP, V[9] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[9]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
171     tempPO = MULWU(tempP); tempP = modadd(tempP, V[10]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[10] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
172     tempPO = MULWU(tempP); tempP = modadd(tempP, V[11]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[11] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
173     tempPO = MULWU(tempP); tempP = modadd(tempP, V[12]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[12] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
174     tempPO = MULWU(tempP); tempP = modadd(tempP, V[13]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[13] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
175     tempPO = MULWU(tempP); tempP = modadd(tempP, V[14]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[14] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
176     tempPO = MULWU(tempP); tempP = modadd(tempP, V[15]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[15] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
177     tempPO = MULWU(tempP); tempP = modadd(tempP, V[16]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[16] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
178     sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(insumtot) + (ovflow <<3 ));
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) >> (BITS-SPECIALMUL))  );
187   }
188   myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
189   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_t xbar)
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, myuint_t b);
198 #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, myuint_t a);
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

Cross-Referencing Geant4 Geant4/externals/clhep/include/CLHEP/Random/MixMaxRng.h

Cross-Referencing Geant4
Geant4/externals/clhep/include/CLHEP/Random/MixMaxRng.h