MixMaxRng.cc 9.49 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
// $Id:$
// -*- C++ -*-
//
// -----------------------------------------------------------------------
//                          HEP Random
//                       --- MixMaxRng ---
//                     class implementation file
// -----------------------------------------------------------------------
//
// This file interfaces the PseudoRandom Number Generator 
// proposed by  N.Z. Akopov, G.K.Saviddy & N.G. Ter-Arutyunian
//   "Matrix Generator of Pseudorandom Numbers"
//       J. Compt. Phy. 97, 573 (1991)
//  Preprint: EPI-867(18)-86, Yerevan June 1986.
//
// Implementation by Konstantin Savvidy - 2004-2015
//        "The MIXMAX random number generator"
//        Comp. Phys. Commun. (2015)
//        http://dx.doi.org/10.1016/j.cpc.2015.06.003
//
//  Release 0.99 and later: released under the LGPL license version 3.0
// =======================================================================
// CLHEP interface implemented by 
//     J. Apostolakis, G. Cosmo & K. Savvidy - Created: 6th July 2015
//     CLHEP interface released under the LGPL license version 3.0
// =======================================================================

#include "CLHEP/Random/defs.h"
#include "CLHEP/Random/Random.h"
#include "CLHEP/Random/MixMaxRng.h"
#include "CLHEP/Random/engineIDulong.h"
#include "CLHEP/Utility/atomic_int.h"

#include <string.h>        // for strcmp
#include <cmath>
#include <cstdlib>
#include <stdint.h>

#include "CLHEP/Random/mixmax.h"

const unsigned long MASK32=0xffffffff;

namespace CLHEP {

namespace {
  // Number of instances with automatic seed selection
  CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
}

static const int MarkerLen = 64; // Enough room to hold a begin or end marker. 

std::string MixMaxRng::name() const { return "MixMaxRng"; } // N=" + N

MixMaxRng::MixMaxRng()
: HepRandomEngine()
{
57
   int numEngines = ++numberOfEngines;
58
   fRngState= rng_alloc();
59
   setSeed(static_cast<long>(numEngines));
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
}

MixMaxRng::MixMaxRng(long seed)
: HepRandomEngine()
{
   fRngState= rng_alloc();
   setSeed(seed);
}

MixMaxRng::MixMaxRng(std::istream& is)
: HepRandomEngine()
{
   get(is);
}

MixMaxRng::~MixMaxRng() 
{
   rng_free( fRngState ); 
}

MixMaxRng::MixMaxRng(const MixMaxRng& rng)
  : HepRandomEngine(rng)
{
   fRngState= rng_copy( rng.fRngState->V );
   fRngState->sumtot= rng.fRngState->sumtot;
   fRngState->counter= rng.fRngState->counter;
}

MixMaxRng& MixMaxRng::operator=(const MixMaxRng& rng)
{
   // Check assignment to self
   //
   if (this == &rng)  { return *this; }

   // Copy base class data
   //
   HepRandomEngine::operator=(rng);

   // Copy data
   //
   rng_free( fRngState );
   fRngState= rng_copy( rng.fRngState->V );
   fRngState->sumtot= rng.fRngState->sumtot;
   fRngState->counter= rng.fRngState->counter;

   return *this;
}

void MixMaxRng::saveStatus( const char filename[] ) const
{
   // Create a C file-handle or an object that can accept the same output
   FILE *fh= fopen(filename, "w");
   if( fh )
   {
     fRngState->fh= fh;
     print_state(fRngState);
     fclose(fh);
   }
   fRngState->fh= 0;
}

void MixMaxRng::restoreStatus( const char filename[] )
{
   read_state(fRngState, filename);
}

void MixMaxRng::showStatus() const
{
   std::cout << std::endl;
   std::cout << "------- MixMaxRng engine status -------" << std::endl;

   std::cout << " Current state vector is:" << std::endl;
   fRngState->fh=stdout;
   print_state(fRngState);
   std::cout << "---------------------------------------" << std::endl;
}

void MixMaxRng::setSeed(long longSeed, int /* extraSeed */)
{
Lynn Garren's avatar
Lynn Garren committed
139
   unsigned long seed0;
140
141
142
143
144
145
146

   theSeed = longSeed;
   if( sizeof(long) > 4)  // C standard says long is at least 32-bits
     seed0= static_cast<unsigned long>(longSeed) & MASK32 ;
   else
     seed0= longSeed;

Lynn Garren's avatar
Lynn Garren committed
147
   seed_spbox(fRngState, seed0);
148
149
150
151
152
153
154
155
156
157
}

//  Preferred Seeding method
//  the values of 'Seeds' must be valid 32-bit integers 
//  Higher order bits will be ignored!!

void MixMaxRng::setSeeds(const long* Seeds, int seedNum)
{
   unsigned long seed0, seed1= 0, seed2= 0, seed3= 0;

Lynn Garren's avatar
Lynn Garren committed
158
159
160
   if( seedNum < 1 ) {  // Assuming at least 2 seeds in vector...
       seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
       seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
   } 
   else
   {
     if( seedNum < 4 ) {
       seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
       if( seedNum > 1){ seed1= static_cast<unsigned long>(Seeds[1]) & MASK32; }
       if( seedNum > 2){ seed2= static_cast<unsigned long>(Seeds[2]) & MASK32; }
     }
     if( seedNum >= 4 ) {
       seed0= static_cast<unsigned long>(Seeds[0]) & MASK32;
       seed1= static_cast<unsigned long>(Seeds[1]) & MASK32;
       seed2= static_cast<unsigned long>(Seeds[2]) & MASK32;
       seed3= static_cast<unsigned long>(Seeds[3]) & MASK32;
     }
   }
   theSeed = Seeds[0];
   theSeeds = Seeds;
Lynn Garren's avatar
Lynn Garren committed
178
   seed_uniquestream(fRngState, seed3, seed2, seed1, seed0);
179
180
181
182
}

double MixMaxRng::flat()
{
183
   return clhep_get_next_float(fRngState);
184
185
}

Lynn Garren's avatar
Lynn Garren committed
186
void MixMaxRng::flatArray(const int size, double* vect )
187
{
Lynn Garren's avatar
Lynn Garren committed
188
189
   // fill_array( fRngState, size, arrayDbl );
   for (int i=0; i<size; ++i) { vect[i] = flat(); }
190
191
192
193
}

MixMaxRng::operator unsigned int()
{
194
195
   return static_cast<unsigned int>(clhep_get_next(fRngState));
   // clhep_get_next returns a 64-bit integer, of which the lower 61 bits
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
   // are random and upper 3 bits are zero
}

std::ostream & MixMaxRng::put ( std::ostream& os ) const
{
   char beginMarker[] = "MixMaxRng-begin";
   char endMarker[]   = "MixMaxRng-end";

   int pr = os.precision(24);
   os << beginMarker << " ";
   os << theSeed << " ";
   for (int i=0; i<rng_get_N(); ++i) {
      os <<  fRngState->V[i] << "\n";
   }
   os << fRngState->counter << "\n";
   os << fRngState->sumtot << "\n";
   os << endMarker << "\n";
   os.precision(pr);
   return os;  
}

std::vector<unsigned long> MixMaxRng::put () const
{
   std::vector<unsigned long> v;
   v.push_back (engineIDulong<MixMaxRng>());
   for (int i=0; i<rng_get_N(); ++i) {
     v.push_back(static_cast<unsigned long>(fRngState->V[i] & MASK32));
       // little-ended order on all platforms
     v.push_back(static_cast<unsigned long>(fRngState->V[i] >> 32  ));
       // pack uint64 into a data structure which is 32-bit on some platforms
   }
   v.push_back(static_cast<unsigned long>(fRngState->counter));
   v.push_back(static_cast<unsigned long>(fRngState->sumtot & MASK32));
   v.push_back(static_cast<unsigned long>(fRngState->sumtot >> 32));
   return v;
}

std::istream & MixMaxRng::get  ( std::istream& is)
{
   char beginMarker [MarkerLen];
   is >> std::ws;
   is.width(MarkerLen);  // causes the next read to the char* to be <=
                         // that many bytes, INCLUDING A TERMINATION \0 
                         // (Stroustrup, section 21.3.2)
   is >> beginMarker;
   if (strcmp(beginMarker,"MixMaxRng-begin")) {
      is.clear(std::ios::badbit | is.rdstate());
      std::cerr << "\nInput stream mispositioned or"
                << "\nMixMaxRng state description missing or"
                << "\nwrong engine type found." << std::endl;
      return is;
   }
   return getState(is);
}

std::string MixMaxRng::beginTag ()
{ 
   return "MixMaxRng-begin"; 
}

std::istream &  MixMaxRng::getState ( std::istream& is )
{
   char endMarker[MarkerLen];
   is >> theSeed;
   for (int i=0; i<rng_get_N(); ++i)  is >> fRngState->V[i];
   is >> fRngState->counter;
   myuint checksum;
   is >> checksum;
   is >> std::ws;
   is.width(MarkerLen);
   is >> endMarker;
   if (strcmp(endMarker,"MixMaxRng-end")) {
       is.clear(std::ios::badbit | is.rdstate());
       std::cerr << "\nMixMaxRng state description incomplete."
                 << "\nInput stream is probably mispositioned now.\n";
       return is;
   }
Lynn Garren's avatar
Lynn Garren committed
273
   if ( fRngState->counter < 0 || fRngState->counter > rng_get_N() ) {
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
       std::cerr << "\nMixMaxRng::getState(): "
                 << "vector read wrong value of counter from file!"
                 << "\nInput stream is probably mispositioned now.\n";
       return is;
   }
   precalc(fRngState);
   if ( checksum != fRngState->sumtot) {
       std::cerr << "\nMixMaxRng::getState(): "
                 << "checksum disagrees with value stored in file!"
                 << "\nInput stream is probably mispositioned now.\n";
       return is;
   }
   return is;
}

bool MixMaxRng::get (const std::vector<unsigned long> & v)
{
   if ((v[0] & 0xffffffffUL) != engineIDulong<MixMaxRng>()) {
     std::cerr << 
        "\nMixMaxRng::get(): vector has wrong ID word - state unchanged\n";
     return false;
   }
   return getState(v);
}

bool MixMaxRng::getState (const std::vector<unsigned long> & v)
{
   if (v.size() != VECTOR_STATE_SIZE ) {
     std::cerr <<
        "\nMixMaxRng::getState(): vector has wrong length - state unchanged\n";
     return false;
   }
   for (int i=1; i<2*rng_get_N() ; i=i+2) {
     fRngState->V[i/2]= ( (v[i] & MASK32) | ( (myuint)(v[i+1]) << 32 ) );
     // unpack from a data structure which is 32-bit on some platforms
   }
   fRngState->counter = v[2*rng_get_N()+1];
   precalc(fRngState);
   if ( ( (v[2*rng_get_N()+2] & MASK32)
        | ( (myuint)(v[2*rng_get_N()+3]) << 32 ) ) != fRngState->sumtot) {
     std::cerr << "\nMixMaxRng::getState(): vector has wrong checksum!"
               << "\nInput vector is probably mispositioned now.\n";
     return false;
   }
   return true;
}

}  // namespace CLHEP