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-2017
31 // =======================================================================
32 
33 #ifndef MixMaxRng_h
34 #define MixMaxRng_h 1
35 
36 #include <array>
38 
39 namespace CLHEP {
40 
41 /**
42  * @author K.Savvidy
43  * @ingroup random
44  */
45 
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(int rowIndex, int colIndex);
59  ~MixMaxRng();
60  // Constructor and destructor.
61 
62  MixMaxRng(const MixMaxRng& rng);
63  MixMaxRng& operator=(const MixMaxRng& rng);
64  // Copy constructor and assignment operator.
65 
66  double flat() { return (S.counter<=(N-1)) ? generate(S.counter):iterate(); }
67  // Returns a pseudo random number between 0 and 1
68  // (excluding the zero: in (0,1] )
69  // smallest number which it will give is approximately 10^-19
70 
71  void flatArray (const int size, double* vect);
72  // Fills the array "vect" of specified size with flat random values.
73 
74  void setSeed(long seed, int dum=0);
75  // Sets the state of the algorithm according to seed.
76 
77  void setSeeds(const long * seeds, int seedNum=0);
78  // Sets the initial state of the engine according to the array of between one and four 32-bit seeds.
79  // If the size of long is greater on the platform, only the lower 32-bits are used.
80  // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
81  // to be independent and non-colliding for at least the next 10^100 random numbers
82 
83  void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
84  // Saves the the current engine state in the file given, by default MixMaxRngState.conf
85 
86  void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
87  // Reads a valid engine state from a given file, by default MixMaxRngState.conf
88  // and restores it.
89 
90  void showStatus() const;
91  // Dumps the engine status on the screen.
92 
93  operator unsigned int();
94  // 32-bit flat
95 
96  virtual std::ostream & put (std::ostream & os) const;
97  virtual std::istream & get (std::istream & is);
98  static std::string beginTag ( );
99  virtual std::istream & getState ( std::istream & is );
100 
101  std::string name() const { return "MixMaxRng"; }
102  static std::string engineName();
103 
104  std::vector<unsigned long> put () const;
105  bool get (const std::vector<unsigned long> & v);
106  bool getState (const std::vector<unsigned long> & v);
107 
108 private:
109 
110  static constexpr long long int SPECIAL = ((N==17)? 0 : ((N==240)? 487013230256099140ULL:0) ); // etc...
111  static constexpr long long int SPECIALMUL= ((N==17)? 36: ((N==240)? 51 :53) ); // etc...
112  // Note the potential for confusion...
113  static constexpr int BITS=61;
114  static constexpr myuint_t M61=2305843009213693951ULL;
115  static constexpr double INV_M61=0.43368086899420177360298E-18;
116  static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX
117 
118  #define MIXMAX_MOD_MERSENNE(k) ((((k)) & M61) + (((k)) >> BITS) )
119 
120  static constexpr int rng_get_N();
121  static constexpr long long int rng_get_SPECIAL();
122  static constexpr int rng_get_SPECIALMUL();
123  void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
124  void seed_spbox(myuint_t);
125  void print_state() const;
126  myuint_t precalc();
127  myuint_t get_next() ;
128  inline double get_next_float() { return get_next_float_packbits(); }
129  // Returns a random double with all 52 bits random, in the range (0,1]
130 
131  MixMaxRng Branch();
132  void BranchInplace(int id);
133 
134  MixMaxRng(myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID ); // Constructor with four 32-bit seeds
135  inline void seed64(myuint_t seedval) { seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (myID_t)seedval ); } // seed with one 64-bit seed
136 
137  double generate(int i);
138  double iterate();
139 
140  double get_next_float_packbits();
141 #if defined __GNUC__
142 #pragma GCC diagnostic push
143 #pragma GCC diagnostic ignored "-Wstrict-aliasing"
144 #endif
145  inline double convert1double(myuint_t u)
146  {
147  const double one = 1;
148  const myuint_t onemask = *(myuint_t*)&one;
149  myuint_t tmp = (u>>9) | onemask; // bits between 52 and 62 dont affect the result!
150  double d = *(double*)&tmp;
151  return d-1.0;
152  }
153 #if defined __GNUC__
154 #pragma GCC diagnostic pop
155 #endif
156  myuint_t MOD_MULSPEC(myuint_t k);
157  myuint_t MULWU(myuint_t k);
158  void seed_vielbein( unsigned int i); // seeds with the i-th unit vector, i = 0..N-1, for testing only
159  myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
160  myuint_t apply_bigskip(myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID );
161  myuint_t modadd(myuint_t foo, myuint_t bar);
162 #if defined(__x86_64__)
163  myuint_t mod128(__uint128_t s);
164  myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b);
165 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
166  myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a);
167 #endif
168 
169 private:
170 
172  {
173  std::array<myuint_t, N> V;
174  myuint_t sumtot;
175  int counter;
176  };
177 
178  typedef struct rng_state_st rng_state_t; // struct alias
180 };
181 
182 } // namespace CLHEP
183 
184 #endif
static constexpr int rng_get_SPECIALMUL()
double convert1double(myuint_t u)
Definition: MixMaxRng.h:145
void BranchInplace(int id)
static constexpr long long int SPECIAL
Definition: MixMaxRng.h:110
unsigned long long int myuint_t
Definition: MixMaxRng.h:47
MixMaxRng Branch()
static std::string beginTag()
unsigned long int myID_t
Definition: MixMaxRng.h:46
myuint_t precalc()
myuint_t apply_bigskip(myuint_t *Vout, myuint_t *Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
MixMaxRng & operator=(const MixMaxRng &rng)
void restoreStatus(const char filename[]="MixMaxRngState.conf")
string filename
Definition: shutoffs.py:106
Float_t tmp
Definition: plot.C:36
std::string name() const
Definition: MixMaxRng.h:101
virtual std::istream & getState(std::istream &is)
static constexpr unsigned int VECTOR_STATE_SIZE
Definition: MixMaxRng.h:116
myuint_t MOD_MULSPEC(myuint_t k)
std::array< myuint_t, N > V
Definition: MixMaxRng.h:173
void seed64(myuint_t seedval)
Definition: MixMaxRng.h:135
Float_t Y
Definition: plot.C:38
static std::string engineName()
const XML_Char * s
Definition: expat.h:262
myuint_t iterate_raw_vec(myuint_t *Y, myuint_t sumtotOld)
void setSeeds(const long *seeds, int seedNum=0)
double generate(int i)
myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a)
unsigned int seed
Definition: runWimpSim.h:102
double get_next_float_packbits()
void seed_vielbein(unsigned int i)
static constexpr double INV_M61
Definition: MixMaxRng.h:115
const double a
void showStatus() const
Float_t d
Definition: plot.C:236
void saveStatus(const char filename[]="MixMaxRngState.conf") const
static constexpr myuint_t M61
Definition: MixMaxRng.h:114
double get_next_float()
Definition: MixMaxRng.h:128
myuint_t MULWU(myuint_t k)
double flat()
Definition: MixMaxRng.h:66
void print_state() const
void flatArray(const int size, double *vect)
static constexpr int rng_get_N()
static constexpr double bar
myuint_t modadd(myuint_t foo, myuint_t bar)
static constexpr long long int rng_get_SPECIAL()
const hit & b
Definition: hits.cxx:21
void seed_uniquestream(myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
static constexpr int BITS
Definition: MixMaxRng.h:113
std::vector< unsigned long > put() const
auto one()
Definition: PMNS.cxx:49
static const int N
Definition: MixMaxRng.h:51
rng_state_t S
Definition: MixMaxRng.h:179
myuint_t get_next()
void setSeed(long seed, int dum=0)
static constexpr long long int SPECIALMUL
Definition: MixMaxRng.h:111
void seed_spbox(myuint_t)
enum BeamMode string