ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MixMaxRng.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MixMaxRng.h
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-2017
31 // =======================================================================
32 
33 #ifndef MixMaxRng_h
34 #define MixMaxRng_h 1
35 
36 #include <array>
37 #include "CLHEP/Random/RandomEngine.h"
38 
39 namespace CLHEP {
40 
46 typedef unsigned long int myID_t;
47 typedef unsigned long long int myuint_t;
48 
49 class MixMaxRng: public HepRandomEngine {
50 
51  static const int N = 17;
52 
53 public:
54 
55  MixMaxRng(std::istream& is);
56  MixMaxRng();
57  MixMaxRng(long seed);
58  ~MixMaxRng();
59  // Constructors and destructor.
60 
61  MixMaxRng(const MixMaxRng& rng);
62  MixMaxRng& operator=(const MixMaxRng& rng);
63  // Copy constructor and assignment operator.
64 
65  double flat() { return (S.counter<=(N-1)) ? generate(S.counter):iterate(); }
66  // Returns a pseudo random number between 0 and 1
67  // (excluding the zero: in (0,1] )
68  // smallest number which it will give is approximately 10^-19
69 
70  void flatArray (const int size, double* vect);
71  // Fills the array "vect" of specified size with flat random values.
72 
73  void setSeed(long seed, int dum=0);
74  // Sets the state of the algorithm according to seed.
75 
76  void setSeeds(const long * seeds, int seedNum=0);
77  // Sets the initial state of the engine according to the array of between one and four 32-bit seeds.
78  // If the size of long is greater on the platform, only the lower 32-bits are used.
79  // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
80  // to be independent and non-colliding for at least the next 10^100 random numbers
81 
82  void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
83  // Saves the the current engine state in the file given, by default MixMaxRngState.conf
84 
85  void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
86  // Reads a valid engine state from a given file, by default MixMaxRngState.conf
87  // and restores it.
88 
89  void showStatus() const;
90  // Dumps the engine status on the screen.
91 
92  operator double();
93  // Returns same as flat()
94  operator float();
95  // less precise flat, faster if possible
96  operator unsigned int();
97  // 32-bit flat
98 
99  virtual std::ostream & put (std::ostream & os) const;
100  virtual std::istream & get (std::istream & is);
101  static std::string beginTag ( );
102  virtual std::istream & getState ( std::istream & is );
103 
104  std::string name() const { return "MixMaxRng"; }
105  static std::string engineName();
106 
107  std::vector<unsigned long> put () const;
108  bool get (const std::vector<unsigned long> & v);
109  bool getState (const std::vector<unsigned long> & v);
110 
111 private:
112 
113  static constexpr long long int SPECIAL = ((N==17)? 0 : ((N==240)? 487013230256099140ULL:0) ); // etc...
114  static constexpr long long int SPECIALMUL= ((N==17)? 36: ((N==240)? 51 :53) ); // etc...
115  // Note the potential for confusion...
116  static constexpr int BITS=61;
117  static constexpr myuint_t M61=2305843009213693951ULL;
118  static constexpr double INV_M61=0.43368086899420177360298E-18;
119  static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX
120 
121  #define MIXMAX_MOD_MERSENNE(k) ((((k)) & M61) + (((k)) >> BITS) )
122 
123  static constexpr int rng_get_N();
124  static constexpr long long int rng_get_SPECIAL();
125  static constexpr int rng_get_SPECIALMUL();
126  void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
127  void seed_spbox(myuint_t);
128  void print_state() const;
129  myuint_t precalc();
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]
133 
134  MixMaxRng Branch();
135  void BranchInplace(int id);
136 
137  MixMaxRng(myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ); // Constructor with four 32-bit seeds
138  inline void seed64(myuint_t seedval) { seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (myID_t)seedval ); } // seed with one 64-bit seed
139 
140  double generate(int i);
141  double iterate();
142 
143  double get_next_float_packbits();
144 #if defined __GNUC__
145 #pragma GCC diagnostic push
146 #pragma GCC diagnostic ignored "-Wstrict-aliasing"
147 #endif
148  inline double convert1double(myuint_t u)
149  {
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;
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
162  myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
163  myuint_t apply_bigskip(myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
164  myuint_t modadd(myuint_t foo, myuint_t bar);
165 #if defined(__x86_64__)
166  myuint_t mod128(__uint128_t s);
167  myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b);
168 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
169  myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a);
170 #endif
171 
172 private:
173 
174  struct rng_state_st
175  {
176  std::array<myuint_t, N> V;
177  myuint_t sumtot;
178  int counter;
179  };
180 
181  typedef struct rng_state_st rng_state_t; // struct alias
182  rng_state_t S;
183 };
184 
185 } // namespace CLHEP
186 
187 #endif