Geant4  10.03
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  int numEngines = ++numberOfEngines;
57  fRngState= rng_alloc();
58  setSeed(static_cast<long>(numEngines));
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;
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_spbox(fRngState, 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(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* vect )
186 {
187  // fill_array( fRngState, size, arrayDbl );
188  for (int i=0; i<size; ++i) { vect[i] = flat(); }
189 }
190 
191 MixMaxRng::operator unsigned int()
192 {
193  return static_cast<unsigned int>(get_next(fRngState));
194  // get_next returns a 64-bit integer, of which the lower 61 bits
195  // are random and upper 3 bits are zero
196 }
197 
198 std::ostream & MixMaxRng::put ( std::ostream& os ) const
199 {
200  char beginMarker[] = "MixMaxRng-begin";
201  char endMarker[] = "MixMaxRng-end";
202 
203  int pr = os.precision(24);
204  os << beginMarker << " ";
205  os << theSeed << " ";
206  for (int i=0; i<rng_get_N(); ++i) {
207  os << fRngState->V[i] << "\n";
208  }
209  os << fRngState->counter << "\n";
210  os << fRngState->sumtot << "\n";
211  os << endMarker << "\n";
212  os.precision(pr);
213  return os;
214 }
215 
216 std::vector<unsigned long> MixMaxRng::put () const
217 {
218  std::vector<unsigned long> v;
219  v.push_back (engineIDulong<MixMaxRng>());
220  for (int i=0; i<rng_get_N(); ++i) {
221  v.push_back(static_cast<unsigned long>(fRngState->V[i] & MASK32));
222  // little-ended order on all platforms
223  v.push_back(static_cast<unsigned long>(fRngState->V[i] >> 32 ));
224  // pack uint64 into a data structure which is 32-bit on some platforms
225  }
226  v.push_back(static_cast<unsigned long>(fRngState->counter));
227  v.push_back(static_cast<unsigned long>(fRngState->sumtot & MASK32));
228  v.push_back(static_cast<unsigned long>(fRngState->sumtot >> 32));
229  return v;
230 }
231 
232 std::istream & MixMaxRng::get ( std::istream& is)
233 {
234  char beginMarker [MarkerLen];
235  is >> std::ws;
236  is.width(MarkerLen); // causes the next read to the char* to be <=
237  // that many bytes, INCLUDING A TERMINATION \0
238  // (Stroustrup, section 21.3.2)
239  is >> beginMarker;
240  if (strcmp(beginMarker,"MixMaxRng-begin")) {
241  is.clear(std::ios::badbit | is.rdstate());
242  std::cerr << "\nInput stream mispositioned or"
243  << "\nMixMaxRng state description missing or"
244  << "\nwrong engine type found." << std::endl;
245  return is;
246  }
247  return getState(is);
248 }
249 
250 std::string MixMaxRng::beginTag ()
251 {
252  return "MixMaxRng-begin";
253 }
254 
255 std::istream & MixMaxRng::getState ( std::istream& is )
256 {
257  char endMarker[MarkerLen];
258  is >> theSeed;
259  for (int i=0; i<rng_get_N(); ++i) is >> fRngState->V[i];
260  is >> fRngState->counter;
261  myuint checksum;
262  is >> checksum;
263  is >> std::ws;
264  is.width(MarkerLen);
265  is >> endMarker;
266  if (strcmp(endMarker,"MixMaxRng-end")) {
267  is.clear(std::ios::badbit | is.rdstate());
268  std::cerr << "\nMixMaxRng state description incomplete."
269  << "\nInput stream is probably mispositioned now.\n";
270  return is;
271  }
272  if ( fRngState->counter < 0 || fRngState->counter > rng_get_N() ) {
273  std::cerr << "\nMixMaxRng::getState(): "
274  << "vector read wrong value of counter from file!"
275  << "\nInput stream is probably mispositioned now.\n";
276  return is;
277  }
278  precalc(fRngState);
279  if ( checksum != fRngState->sumtot) {
280  std::cerr << "\nMixMaxRng::getState(): "
281  << "checksum disagrees with value stored in file!"
282  << "\nInput stream is probably mispositioned now.\n";
283  return is;
284  }
285  return is;
286 }
287 
288 bool MixMaxRng::get (const std::vector<unsigned long> & v)
289 {
290  if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
291  std::cerr <<
292  "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
293  return false;
294  }
295  return getState(v);
296 }
297 
298 bool MixMaxRng::getState (const std::vector<unsigned long> & v)
299 {
300  if (v.size() != VECTOR_STATE_SIZE ) {
301  std::cerr <<
302  "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
303  return false;
304  }
305  for (int i=1; i<2*rng_get_N() ; i=i+2) {
306  fRngState->V[i/2]= ( (v[i] & MASK32) | ( (myuint)(v[i+1]) << 32 ) );
307  // unpack from a data structure which is 32-bit on some platforms
308  }
309  fRngState->counter = v[2*rng_get_N()+1];
310  precalc(fRngState);
311  if ( ( (v[2*rng_get_N()+2] & MASK32)
312  | ( (myuint)(v[2*rng_get_N()+3]) << 32 ) ) != fRngState->sumtot) {
313  std::cerr << "\nMixMaxRng::getState(): vector has wrong checksum!"
314  << "\nInput vector is probably mispositioned now.\n";
315  return false;
316  }
317  return true;
318 }
319 
320 } // 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
int rng_free(rng_state_t *X)
Definition: mixmax.cc:159
const char * name(G4int ptype)
void seed_spbox(rng_state_t *X, myuint seed)
Definition: mixmax.cc:202
double get_next_float(rng_state_t *X)
Definition: mixmax.cc:86
long seed
Definition: chem4.cc:68
void seed_uniquestream(rng_state_t *Xin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t streamID)
Definition: mixmax.cc:435
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 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