Geant4  10.02.p01
MixMaxRng.cc
Go to the documentation of this file.
1 // $Id:$
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- MixMaxRng ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 //
10 // This file interfaces the PseudoRandom Number Generator
11 // proposed by N.Z. Akopov, G.K.Saviddy & N.G. Ter-Arutyunian
12 // "Matrix Generator of Pseudorandom Numbers"
13 // J. Compt. Phy. 97, 573 (1991)
14 // Preprint: EPI-867(18)-86, Yerevan June 1986.
15 //
16 // Implementation by Konstantin Savvidy - 2004-2015
17 // "The MIXMAX random number generator"
18 // Comp. Phys. Commun. (2015)
19 // http://dx.doi.org/10.1016/j.cpc.2015.06.003
20 //
21 // Release 0.99 and later: released under the LGPL license version 3.0
22 // =======================================================================
23 // CLHEP interface implemented by
24 // J. Apostolakis, G. Cosmo & K. Savvidy - Created: 6th July 2015
25 // CLHEP interface released under the LGPL license version 3.0
26 // =======================================================================
27 
28 #include "CLHEP/Random/Random.h"
29 #include "CLHEP/Random/MixMaxRng.h"
30 #include "CLHEP/Random/engineIDulong.h"
31 #include "CLHEP/Utility/atomic_int.h"
32 
33 #include <string.h> // for strcmp
34 #include <cmath>
35 #include <cstdlib>
36 #include <stdint.h>
37 
38 #include "CLHEP/Random/mixmax.h"
39 
40 const unsigned long MASK32=0xffffffff;
41 
42 namespace CLHEP {
43 
44 namespace {
45  // Number of instances with automatic seed selection
46  CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
47 }
48 
49 static const int MarkerLen = 64; // Enough room to hold a begin or end marker.
50 
51 std::string MixMaxRng::name() const { return "MixMaxRng"; } // N=" + N
52 
53 MixMaxRng::MixMaxRng()
54 : HepRandomEngine()
55 {
56  numberOfEngines++;
57  fRngState= rng_alloc();
58  setSeed(static_cast<long>(numberOfEngines));
59 }
60 
61 MixMaxRng::MixMaxRng(long seed)
62 : HepRandomEngine()
63 {
64  fRngState= rng_alloc();
65  setSeed(seed);
66 }
67 
68 MixMaxRng::MixMaxRng(std::istream& is)
69 : HepRandomEngine()
70 {
71  get(is);
72 }
73 
74 MixMaxRng::~MixMaxRng()
75 {
76  rng_free( fRngState );
77 }
78 
79 MixMaxRng::MixMaxRng(const MixMaxRng& rng)
80  : HepRandomEngine(rng)
81 {
82  fRngState= rng_copy( rng.fRngState->V );
83  fRngState->sumtot= rng.fRngState->sumtot;
84  fRngState->counter= rng.fRngState->counter;
85 }
86 
87 MixMaxRng& MixMaxRng::operator=(const MixMaxRng& rng)
88 {
89  // Check assignment to self
90  //
91  if (this == &rng) { return *this; }
92 
93  // Copy base class data
94  //
95  HepRandomEngine::operator=(rng);
96 
97  // Copy data
98  //
99  rng_free( fRngState );
100  fRngState= rng_copy( rng.fRngState->V );
101  fRngState->sumtot= rng.fRngState->sumtot;
102  fRngState->counter= rng.fRngState->counter;
103 
104  return *this;
105 }
106 
107 void MixMaxRng::saveStatus( const char filename[] ) const
108 {
109  // Create a C file-handle or an object that can accept the same output
110  FILE *fh= fopen(filename, "w");
111  if( fh )
112  {
113  fRngState->fh= fh;
114  print_state(fRngState);
115  fclose(fh);
116  }
117  fRngState->fh= 0;
118 }
119 
120 void MixMaxRng::restoreStatus( const char filename[] )
121 {
122  read_state(fRngState, filename);
123 }
124 
125 void MixMaxRng::showStatus() const
126 {
127  std::cout << std::endl;
128  std::cout << "------- MixMaxRng engine status -------" << std::endl;
129 
130  std::cout << " Current state vector is:" << std::endl;
131  fRngState->fh=stdout;
132  print_state(fRngState);
133  std::cout << "---------------------------------------" << std::endl;
134 }
135 
136 void MixMaxRng::setSeed(long longSeed, int /* extraSeed */)
137 {
138  unsigned long seed0, seed1= 0, seed2= 0, seed3= 0;
139 
140  theSeed = longSeed;
141  if( sizeof(long) > 4) // C standard says long is at least 32-bits
142  seed0= static_cast<unsigned long>(longSeed) & MASK32 ;
143  else
144  seed0= longSeed;
145 
146  seed_uniquestream(this->fRngState, seed3, seed2, seed1, seed0 );
147 }
148 
149 // Preferred Seeding method
150 // the values of 'Seeds' must be valid 32-bit integers
151 // Higher order bits will be ignored!!
152 
153 void MixMaxRng::setSeeds(const long* Seeds, int seedNum)
154 {
155  unsigned long seed0, seed1= 0, seed2= 0, seed3= 0;
156 
157  if( seedNum < 1 ) { // Assuming at least 2 seeds in vector...
158  seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
159  seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
160  }
161  else
162  {
163  if( seedNum < 4 ) {
164  seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
165  if( seedNum > 1){ seed1= static_cast<unsigned long>(Seeds[1]) & MASK32; }
166  if( seedNum > 2){ seed2= static_cast<unsigned long>(Seeds[2]) & MASK32; }
167  }
168  if( seedNum >= 4 ) {
169  seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
170  seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
171  seed2= static_cast<unsigned long>(Seeds[2]) & MASK32;
172  seed3= static_cast<unsigned long>(Seeds[3]) & MASK32;
173  }
174  }
175  theSeed = Seeds[0];
176  theSeeds = Seeds;
177  seed_uniquestream(this->fRngState, seed3, seed2, seed1, seed0 );
178 }
179 
180 double MixMaxRng::flat()
181 {
182  return get_next_float(fRngState);
183 }
184 
185 void MixMaxRng::flatArray(const int size, double* arrayDbl )
186 {
187  fill_array( fRngState, size, arrayDbl );
188 }
189 
190 MixMaxRng::operator unsigned int()
191 {
192  return static_cast<unsigned int>(get_next(fRngState));
193  // get_next returns a 64-bit integer, of which the lower 61 bits
194  // are random and upper 3 bits are zero
195 }
196 
197 std::ostream & MixMaxRng::put ( std::ostream& os ) const
198 {
199  char beginMarker[] = "MixMaxRng-begin";
200  char endMarker[] = "MixMaxRng-end";
201 
202  int pr = os.precision(24);
203  os << beginMarker << " ";
204  os << theSeed << " ";
205  for (int i=0; i<rng_get_N(); ++i) {
206  os << fRngState->V[i] << "\n";
207  }
208  os << fRngState->counter << "\n";
209  os << fRngState->sumtot << "\n";
210  os << endMarker << "\n";
211  os.precision(pr);
212  return os;
213 }
214 
215 std::vector<unsigned long> MixMaxRng::put () const
216 {
217  std::vector<unsigned long> v;
218  v.push_back (engineIDulong<MixMaxRng>());
219  for (int i=0; i<rng_get_N(); ++i) {
220  v.push_back(static_cast<unsigned long>(fRngState->V[i] & MASK32));
221  // little-ended order on all platforms
222  v.push_back(static_cast<unsigned long>(fRngState->V[i] >> 32 ));
223  // pack uint64 into a data structure which is 32-bit on some platforms
224  }
225  v.push_back(static_cast<unsigned long>(fRngState->counter));
226  v.push_back(static_cast<unsigned long>(fRngState->sumtot & MASK32));
227  v.push_back(static_cast<unsigned long>(fRngState->sumtot >> 32));
228  return v;
229 }
230 
231 std::istream & MixMaxRng::get ( std::istream& is)
232 {
233  char beginMarker [MarkerLen];
234  is >> std::ws;
235  is.width(MarkerLen); // causes the next read to the char* to be <=
236  // that many bytes, INCLUDING A TERMINATION \0
237  // (Stroustrup, section 21.3.2)
238  is >> beginMarker;
239  if (strcmp(beginMarker,"MixMaxRng-begin")) {
240  is.clear(std::ios::badbit | is.rdstate());
241  std::cerr << "\nInput stream mispositioned or"
242  << "\nMixMaxRng state description missing or"
243  << "\nwrong engine type found." << std::endl;
244  return is;
245  }
246  return getState(is);
247 }
248 
249 std::string MixMaxRng::beginTag ()
250 {
251  return "MixMaxRng-begin";
252 }
253 
254 std::istream & MixMaxRng::getState ( std::istream& is )
255 {
256  char endMarker[MarkerLen];
257  is >> theSeed;
258  for (int i=0; i<rng_get_N(); ++i) is >> fRngState->V[i];
259  is >> fRngState->counter;
260  myuint checksum;
261  is >> checksum;
262  is >> std::ws;
263  is.width(MarkerLen);
264  is >> endMarker;
265  if (strcmp(endMarker,"MixMaxRng-end")) {
266  is.clear(std::ios::badbit | is.rdstate());
267  std::cerr << "\nMixMaxRng state description incomplete."
268  << "\nInput stream is probably mispositioned now.\n";
269  return is;
270  }
271  if ( fRngState->counter < 0 || fRngState->counter > N-1 ) {
272  std::cerr << "\nMixMaxRng::getState(): "
273  << "vector read wrong value of counter from file!"
274  << "\nInput stream is probably mispositioned now.\n";
275  return is;
276  }
277  precalc(fRngState);
278  if ( checksum != fRngState->sumtot) {
279  std::cerr << "\nMixMaxRng::getState(): "
280  << "checksum disagrees with value stored in file!"
281  << "\nInput stream is probably mispositioned now.\n";
282  return is;
283  }
284  return is;
285 }
286 
287 bool MixMaxRng::get (const std::vector<unsigned long> & v)
288 {
289  if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
290  std::cerr <<
291  "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
292  return false;
293  }
294  return getState(v);
295 }
296 
297 bool MixMaxRng::getState (const std::vector<unsigned long> & v)
298 {
299  if (v.size() != VECTOR_STATE_SIZE ) {
300  std::cerr <<
301  "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
302  return false;
303  }
304  for (int i=1; i<2*rng_get_N() ; i=i+2) {
305  fRngState->V[i/2]= ( (v[i] & MASK32) | ( (myuint)(v[i+1]) << 32 ) );
306  // unpack from a data structure which is 32-bit on some platforms
307  }
308  fRngState->counter = v[2*rng_get_N()+1];
309  precalc(fRngState);
310  if ( ( (v[2*rng_get_N()+2] & MASK32)
311  | ( (myuint)(v[2*rng_get_N()+3]) << 32 ) ) != fRngState->sumtot) {
312  std::cerr << "\nMixMaxRng::getState(): vector has wrong checksum!"
313  << "\nInput vector is probably mispositioned now.\n";
314  return false;
315  }
316  return true;
317 }
318 
319 } // namespace CLHEP
void print_state(rng_state_t *X)
Definition: mixmax.cc:270
static const int MarkerLen
Definition: DualRand.cc:67
int rng_get_N(void)
Definition: mixmax.cc:237
rng_state_t * rng_copy(myuint *Y)
Definition: mixmax.cc:165
myuint get_next(rng_state_t *X)
Definition: mixmax.cc:82
const unsigned long MASK32
Definition: MixMaxRng.cc:40
G4String name
Definition: TRTMaterials.hh:40
int rng_free(rng_state_t *X)
Definition: mixmax.cc:159
double get_next_float(rng_state_t *X)
Definition: mixmax.cc:86
void seed_uniquestream(rng_state_t *Xin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
Definition: mixmax.cc:431
double flat()
Definition: G4AblaRandom.cc:47
rng_state_t * rng_alloc()
Definition: mixmax.cc:151
myuint precalc(rng_state_t *X)
Definition: mixmax.cc:225
void fill_array(rng_state_t *X, unsigned int n, double *array)
Definition: mixmax.cc:90
void read_state(rng_state_t *X, const char filename[])
Definition: mixmax.cc:283
void setSeeds(const SeedVector &sv)
Set the seeds of the current generator.
Definition: G4INCLRandom.cc:85