Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/MixMaxRng.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/clhep/src/MixMaxRng.cc (Version 11.3.0) and /externals/clhep/src/MixMaxRng.cc (Version 9.1.p1)


  1 //
                                                 1 
  2 // -*- C++ -*-
                                   
  3 //
                                               
  4 // -------------------------------------------    
  5 //                          HEP Random
           
  6 //                       --- MixMaxRng ---
       
  7 //                     class implementation fi    
  8 // -------------------------------------------    
  9 //
                                               
 10 // This file interfaces the MixMax PseudoRando    
 11 // proposed by:
                                  
 12 //
                                               
 13 // G.K.Savvidy and N.G.Ter-Arutyunian,
           
 14 //   On the Monte Carlo simulation of physical    
 15 //   J.Comput.Phys. 97, 566 (1991);
              
 16 //   Preprint EPI-865-16-86, Yerevan, Jan. 198    
 17 //   http://dx.doi.org/10.1016/0021-9991(91)90    
 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.0    
 23 //
                                               
 24 // K.Savvidy and G.Savvidy
                       
 25 //   "Spectrum and Entropy of C-systems. MIXMA    
 26 //   Chaos, Solitons & Fractals, Volume 91, (2    
 27 //   http://dx.doi.org/10.1016/j.chaos.2016.05    
 28 //
                                               
 29 // ===========================================    
 30 // Implementation by Konstantin Savvidy - Copy    
 31 // July 2023 - Updated class structure upon su    
 32 // September 2023 - fix (re-)initialization fr    
 33 // ===========================================    
 34 
                                                 
 35 #include "CLHEP/Random/Random.h"
                 
 36 #include "CLHEP/Random/MixMaxRng.h"
              
 37 #include "CLHEP/Random/engineIDulong.h"
          
 38 #include "CLHEP/Utility/atomic_int.h"
            
 39 
                                                 
 40 #include <atomic>
                                
 41 #include <cmath>
                                 
 42 #include <algorithm>
                             
 43 #include <iostream>
                              
 44 #include <string.h>        // for strcmp
         
 45 #include <vector>
                                
 46 
                                                 
 47 const unsigned long MASK32=0xffffffff;
           
 48 
                                                 
 49 namespace CLHEP {
                                
 50 
                                                 
 51 namespace {
                                      
 52   // Number of instances with automatic seed s    
 53   CLHEP_ATOMIC_INT_TYPE numberOfEngines(0);
      
 54 }
                                                
 55 
                                                 
 56 static const int MarkerLen = 64; // Enough roo    
 57 
                                                 
 58 MixMaxRng::MixMaxRng()
                           
 59 : HepRandomEngine()
                              
 60 {
                                                
 61    int numEngines = ++numberOfEngines;
           
 62    setSeed(static_cast<long>(numEngines));
       
 63 }
                                                
 64 
                                                 
 65 MixMaxRng::MixMaxRng(long seed)
                  
 66 : HepRandomEngine()
                              
 67 {
                                                
 68    theSeed=seed;
                                 
 69    setSeed(seed);
                                
 70 }
                                                
 71 
                                                 
 72 MixMaxRng::MixMaxRng(std::istream& is)
           
 73 : HepRandomEngine()
                              
 74 {
                                                
 75    get(is);
                                      
 76 }
                                                
 77 
                                                 
 78 MixMaxRng::~MixMaxRng() 
                         
 79 {
                                                
 80 }
                                                
 81 
                                                 
 82 MixMaxRng::MixMaxRng(const MixMaxRng& rng)
       
 83   : HepRandomEngine(rng)
                         
 84 {
                                                
 85    std::copy(rng.V, rng.V+N, V);
                 
 86    sumtot= rng.sumtot;
                           
 87    counter= rng.counter;
                         
 88 }
                                                
 89 
                                                 
 90 MixMaxRng& MixMaxRng::operator=(const MixMaxRn    
 91 {
                                                
 92    // Check assignment to self
                   
 93    //
                                            
 94    if (this == &rng)  { return *this; }
          
 95 
                                                 
 96    // Copy base class data
                       
 97    //
                                            
 98    HepRandomEngine::operator=(rng);
              
 99 
                                                 
100    std::copy(rng.V, rng.V+N, V);
                 
101    sumtot= rng.sumtot;
                           
102    counter= rng.counter;
                         
103 
                                                 
104    return *this;
                                 
105 }
                                                
106 
                                                 
107 void MixMaxRng::saveStatus( const char filenam    
108 {
                                                
109    // Create a C file-handle or an object that    
110    FILE *fh= fopen(filename, "w");
               
111    if( fh )
                                      
112    {
                                             
113      int j;
                                      
114      fprintf(fh, "mixmax state, file version 1    
115      fprintf(fh, "N=%u; V[N]={", rng_get_N() )    
116      for (j=0; (j< (rng_get_N()-1) ); ++j) {
     
117          fprintf(fh, "%llu, ", (unsigned long     
118      }
                                           
119      fprintf(fh, "%llu", (unsigned long long)V    
120      fprintf(fh, "}; " );
                        
121      fprintf(fh, "counter=%u; ", counter );
      
122      fprintf(fh, "sumtot=%llu;\n", (unsigned l    
123      fclose(fh);
                                 
124    }
                                             
125 }
                                                
126 
                                                 
127 void MixMaxRng::restoreStatus( const char file    
128 {
                                                
129    // a function for reading the state from a     
130    FILE* fin;
                                    
131    if(  ( fin = fopen(filename, "r") ) )
         
132    {
                                             
133       char l=0;
                                  
134       while ( l != '{' ) { // 0x7B = "{"
         
135         l=fgetc(fin); // proceed until hitting    
136       }
                                          
137       ungetc(' ', fin);
                          
138    }
                                             
139    else
                                          
140    {
                                             
141       fprintf(stderr, "mixmax -> read_state: e    
142       throw std::runtime_error("Error in readi    
143    }
                                             
144     
                                             
145    myuint_t vecVal;
                              
146    //printf("mixmax -> read_state: starting to    
147    if (!fscanf(fin, "%llu", (unsigned long lon    
148    {
                                             
149      fprintf(stderr, "mixmax -> read_state: er    
150      throw std::runtime_error("Error in readin    
151    }
                                             
152 
                                                 
153    for (int i = 1; i < rng_get_N(); ++i)
         
154    {
                                             
155      if (!fscanf(fin, ", %llu", (unsigned long    
156      {
                                           
157        fprintf(stderr, "mixmax -> read_state:     
158        throw std::runtime_error("Error in read    
159      }
                                           
160      if( vecVal <= MixMaxRng::M61 )
              
161      {
                                           
162        V[i] = vecVal;
                            
163      }
                                           
164      else
                                        
165      {
                                           
166        fprintf(stderr, "mixmax -> read_state:     
167                " ( must be less than %llu ) "
    
168                " obtained from reading file %s    
169                , (unsigned long long)vecVal, (    
170      }
                                           
171    }
                                             
172     
                                             
173    int incounter;
                                
174    if (!fscanf( fin, "}; counter=%i; ", &incou    
175    {
                                             
176      fprintf(stderr, "mixmax -> read_state: er    
177      throw std::runtime_error("Error in readin    
178    }
                                             
179    if( incounter <= rng_get_N() )
                
180    {
                                             
181      counter = incounter;
                        
182    }
                                             
183    else
                                          
184    {
                                             
185      fprintf(stderr, "mixmax -> read_state: In    
186              "  Must be 0 <= counter < %u\n" ,    
187      print_state();
                              
188      throw std::runtime_error("Error in readin    
189    }
                                             
190    precalc();
                                    
191    myuint_t insumtot;
                            
192    if (!fscanf( fin, "sumtot=%llu\n", (unsigne    
193    {
                                             
194      fprintf(stderr, "mixmax -> read_state: er    
195      throw std::runtime_error("Error in readin    
196    }
                                             
197 
                                                 
198    if (sumtot != insumtot)
                       
199    {
                                             
200      fprintf(stderr, "mixmax -> checksum error    
201      throw std::runtime_error("Error in readin    
202    }
                                             
203    fclose(fin);
                                  
204 }
                                                
205 
                                                 
206 void MixMaxRng::showStatus() const
               
207 {
                                                
208    std::cout << std::endl;
                       
209    std::cout << "------- MixMaxRng engine stat    
210 
                                                 
211    std::cout << " Current state vector is:" <<    
212    print_state();
                                
213    std::cout << "-----------------------------    
214 }
                                                
215 
                                                 
216 //  Preferred Seeding method
                     
217 //  the values of 'Seeds' must be valid 32-bit    
218 //  Higher order bits will be ignored!!
          
219 
                                                 
220 void MixMaxRng::setSeeds(const long* Seeds, in    
221 {
                                                
222    myID_t seed0, seed1= 0, seed2= 0, seed3= 0;    
223 
                                                 
224    if( seedNum < 1 ) {  // Assuming at least 2    
225        seed0= static_cast<myID_t>(Seeds[0]) &     
226        seed1= static_cast<myID_t>(Seeds[1]) &     
227    } 
                                            
228    else
                                          
229    {
                                             
230      if( seedNum < 4 ) {
                         
231        seed0= static_cast<myID_t>(Seeds[0]) &     
232        if( seedNum > 1){ seed1= static_cast<my    
233        if( seedNum > 2){ seed2= static_cast<my    
234      }
                                           
235      if( seedNum >= 4 ) {
                        
236        seed0= static_cast<myID_t>(Seeds[0]) &     
237        seed1= static_cast<myID_t>(Seeds[1]) &     
238        seed2= static_cast<myID_t>(Seeds[2]) &     
239        seed3= static_cast<myID_t>(Seeds[3]) &     
240      }
                                           
241    }
                                             
242    theSeed = Seeds[0];
                           
243    theSeeds = Seeds;
                             
244    seed_uniquestream(seed3, seed2, seed1, seed    
245 }
                                                
246 
                                                 
247 std::string MixMaxRng::engineName()
              
248 {
                                                
249    return "MixMaxRng";
                           
250 }
                                                
251 
                                                 
252 constexpr int MixMaxRng::rng_get_N()
             
253 {
                                                
254    return N;
                                     
255 }
                                                
256 
                                                 
257 void MixMaxRng::flatArray(const int size, doub    
258 {
                                                
259    for (int i=0; i<size; ++i) { vect[i] = flat    
260 }
                                                
261 
                                                 
262 std::ostream & MixMaxRng::put ( std::ostream&     
263 {
                                                
264    char beginMarker[] = "MixMaxRng-begin";
       
265    char endMarker[]   = "MixMaxRng-end";
         
266 
                                                 
267    long pr = os.precision(24);
                   
268    os << beginMarker << " ";
                     
269    os << theSeed << "\n";
                        
270    for (int i=0; i<rng_get_N(); ++i) {
           
271       os <<  V[i] << "\n";
                       
272    }
                                             
273    os << counter << "\n";
                        
274    os << sumtot << "\n";
                         
275    os << endMarker << "\n";
                      
276    os.precision(pr);
                             
277    return os;  
                                  
278 }
                                                
279 
                                                 
280 std::vector<unsigned long> MixMaxRng::put () c    
281 {
                                                
282    std::vector<unsigned long> vec;
               
283    vec.push_back (engineIDulong<MixMaxRng>());    
284    for (int i=0; i<rng_get_N(); ++i)
             
285    {
                                             
286      vec.push_back(static_cast<unsigned long>(    
287        // little-ended order on all platforms
    
288      vec.push_back(static_cast<unsigned long>(    
289        // pack uint64 into a data structure wh    
290    }
                                             
291    vec.push_back(static_cast<unsigned long>(co    
292    vec.push_back(static_cast<unsigned long>(su    
293    vec.push_back(static_cast<unsigned long>(su    
294    return vec;
                                   
295 }
                                                
296 
                                                 
297 std::istream & MixMaxRng::get  ( std::istream&    
298 {
                                                
299    char beginMarker [MarkerLen];
                 
300    is >> std::ws;
                                
301    is.width(MarkerLen);  // causes the next re    
302                          // that many bytes, I    
303                          // (Stroustrup, secti    
304    is >> beginMarker;
                            
305    if (strcmp(beginMarker,"MixMaxRng-begin"))     
306       is.clear(std::ios::badbit | is.rdstate()    
307       std::cerr << "\nInput stream misposition    
308                 << "\nMixMaxRng state descript    
309                 << "\nwrong engine type found.    
310       return is;
                                 
311    }
                                             
312    return getState(is);
                          
313 }
                                                
314 
                                                 
315 std::string MixMaxRng::beginTag ()
               
316 { 
                                               
317    return "MixMaxRng-begin"; 
                    
318 }
                                                
319 
                                                 
320 std::istream &  MixMaxRng::getState ( std::ist    
321 {
                                                
322    char endMarker[MarkerLen];
                    
323    is >> theSeed;
                                
324    for (int i=0; i<rng_get_N(); ++i)  is >> V[    
325    is >> counter;
                                
326    myuint_t checksum;
                            
327    is >> checksum;
                               
328    is >> std::ws;
                                
329    is.width(MarkerLen);
                          
330    is >> endMarker;
                              
331    if (strcmp(endMarker,"MixMaxRng-end")) {
      
332        is.clear(std::ios::badbit | is.rdstate(    
333        std::cerr << "\nMixMaxRng state descrip    
334                  << "\nInput stream is probabl    
335        return is;
                                
336    }
                                             
337    if ( counter < 0 || counter > rng_get_N() )    
338        std::cerr << "\nMixMaxRng::getState():     
339                  << "vector read wrong value o    
340                  << "\nInput stream is probabl    
341        return is;
                                
342    }
                                             
343    precalc();
                                    
344    if ( checksum != sumtot) {
                    
345        std::cerr << "\nMixMaxRng::getState():     
346                  << "checksum disagrees with v    
347                  << "\nInput stream is probabl    
348        return is;
                                
349    }
                                             
350    return is;
                                    
351 }
                                                
352 
                                                 
353 bool MixMaxRng::get (const std::vector<unsigne    
354 {
                                                
355    if ((vec[0] & 0xffffffffUL) != engineIDulon    
356    {
                                             
357      std::cerr << 
                               
358         "\nMixMaxRng::get(): vector has wrong     
359      return false;
                               
360    }
                                             
361    return getState(vec);
                         
362 }
                                                
363 
                                                 
364 bool MixMaxRng::getState (const std::vector<un    
365 {
                                                
366    if (vec.size() != VECTOR_STATE_SIZE ) {
       
367      std::cerr <<
                                
368         "\nMixMaxRng::getState(): vector has w    
369      return false;
                               
370    }
                                             
371    for (int i=1; i<2*rng_get_N() ; i=i+2) {
      
372      V[i/2]= ( (vec[i] & MASK32) | ( (myuint_t    
373      // unpack from a data structure which is     
374    }
                                             
375    counter = (int)vec[2*rng_get_N()+1];
          
376    precalc();
                                    
377    if ( ( (vec[2*rng_get_N()+2] & MASK32)
        
378         | ( (myuint_t)(vec[2*rng_get_N()+3]) <    
379      std::cerr << "\nMixMaxRng::getState(): ve    
380                << "\nInput vector is probably     
381      return false;
                               
382    }
                                             
383    return true;
                                  
384 }
                                                
385 
                                                 
386 myuint_t MixMaxRng::iterate_raw_vec(myuint_t*     
387 {
                                                
388    // operates with a raw vector, uses known s    
389             
                                     
390    myuint_t  tempP, tempV;
                       
391    Y[0] = ( tempV = sumtotOld);
                  
392    myuint_t insumtot = Y[0], ovflow = 0; // wi    
393    tempP = 0;              // will keep a part    
394    for (int i=1; (i<N); ++i)
                     
395    {
                                             
396      myuint_t tempPO = MULWU(tempP);
             
397      tempP = modadd(tempP, Y[i]);
                
398      tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+t    
399      Y[i] = tempV;
                               
400      insumtot += tempV; if (insumtot < tempV)     
401    }
                                             
402    return MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSE    
403 }
                                                
404         
                                         
405 myuint_t MixMaxRng::get_next()
                   
406 {
                                                
407    int i = counter;
                              
408             
                                     
409    if ((i<=(N-1)) )
                              
410    {
                                             
411      ++counter;
                                  
412      return V[i];
                                
413    }
                                             
414    else
                                          
415    {
                                             
416      sumtot = iterate_raw_vec(V, sumtot);
        
417      counter=2;
                                  
418      return V[1];
                                
419    }
                                             
420 }
                                                
421 
                                                 
422 myuint_t MixMaxRng::precalc()
                    
423 {
                                                
424    myuint_t temp = 0;
                            
425    for (int i=0; i < N; ++i) { temp = MIXMAX_M    
426    sumtot = temp;
                                
427    return temp;
                                  
428 }
                                                
429         
                                         
430 void MixMaxRng::state_init()
                     
431 {
                                                
432    for (int i=1; i < N; ++i) { V[i] = 0; }
       
433    V[0] = 1;
                                     
434    counter = N;  // set the counter to N if it    
435    sumtot = 1;
                                   
436 }
                                                
437 
                                                 
438 void MixMaxRng::seed_spbox(myuint_t seed)
        
439 {
                                                
440    // a 64-bit LCG from Knuth line 26, in comb    
441 
                                                 
442    if (seed == 0)
                                
443      throw std::runtime_error("try seeding wit    
444 
                                                 
445    const myuint_t MULT64=6364136223846793005UL    
446    sumtot = 0;
                                   
447             
                                     
448    myuint_t l = seed;
                            
449             
                                     
450    for (int i=0; i < N; ++i)
                     
451    {
                                             
452      l*=MULT64; l = (l << 32) ^ (l>>32);
         
453      V[i] = l & M61;
                             
454      sumtot = MIXMAX_MOD_MERSENNE(sumtot + V[(    
455    }
                                             
456    counter = N;  // set the counter to N if it    
457 }
                                                
458 
                                                 
459 void MixMaxRng::seed_uniquestream( myID_t clus    
460 {
                                                
461    state_init();
                                 
462    sumtot = apply_bigskip(V, V, clusterID, mac    
463    counter = 1;
                                  
464 }
                                                
465         
                                         
466 myuint_t MixMaxRng::apply_bigskip( myuint_t* V    
467 {
                                                
468   /*
                                             
469     makes a derived state vector, Vout, from t    
470     by skipping a large number of steps, deter    
471              
                                    
472     it is mathematically guaranteed that the s    
473     1) at least one bit of ID is different
       
474     2) less than 10^100 numbers are drawn from    
475     (this is good enough : a single CPU will n    
476     even if it had a clock cycle of Planch tim    
477              
                                    
478     Caution: never apply this to a derived vec    
479     and use it in all your runs, just change r    
480              
                                    
481     clusterID and machineID are provided for t    
482     which is running in parallel on a large nu    
483              
                                    
484     did i repeat it enough times? the non-coll    
485 
                                                 
486    */
                                            
487 
                                                 
488    const myuint_t skipMat17[128][17] =
           
489    #include "CLHEP/Random/mixmax_skip_N17.icc"    
490    ;
                                             
491             
                                     
492    const myuint_t* skipMat[128];
                 
493    for (int i=0; i<128; ++i) { skipMat[i] = sk    
494             
                                     
495    myID_t IDvec[4] = {streamID, runID, machine    
496    int r,i,j,  IDindex;
                          
497    myID_t id;
                                    
498    myuint_t Y[N], cum[N];
                        
499    myuint_t coeff;
                               
500    myuint_t* rowPtr;
                             
501    myuint_t insumtot=0;
                          
502             
                                     
503    for (i=0; i<N; ++i) { Y[i] = Vin[i]; insumt    
504    for (IDindex=0; IDindex<4; ++IDindex)
         
505    { // go from lower order to higher order ID    
506      id=IDvec[IDindex];
                          
507      //printf("now doing ID at level %d, with     
508      r = 0;
                                      
509      while (id)
                                  
510      {
                                           
511        if (id & 1)
                               
512        {
                                         
513          rowPtr = (myuint_t*)skipMat[r + IDind    
514          for (i=0; i<N; ++i){ cum[i] = 0; }
      
515          for (j=0; j<N; ++j)
                     
516          { // j is lag, enumerates terms of th    
517            // for zero lag Y is already given
    
518            coeff = rowPtr[j]; // same coeff fo    
519            for (i =0; i<N; ++i){
                 
520              cum[i] =  fmodmulM61( cum[i], coe    
521            }
                                     
522            insumtot = iterate_raw_vec(Y, insum    
523          }
                                       
524          insumtot=0;
                             
525          for (i=0; i<N; ++i){ Y[i] = cum[i]; i    
526        }
                                         
527        id = (id >> 1); ++r; // bring up the r-    
528      }
                                           
529    }
                                             
530    insumtot=0;
                                   
531    for (i=0; i<N; ++i){ Vout[i] = Y[i]; insumt    
532    // returns sumtot, and copy the vector over    
533    return insumtot;
                              
534 }
                                                
535         
                                         
536 #if defined(__x86_64__)
                          
537   myuint_t MixMaxRng::mod128(__uint128_t s)
      
538   {
                                              
539     myuint_t s1;
                                 
540     s1 = ( (  ((myuint_t)s)&M61 ) + (  ((myuin    
541     return  MIXMAX_MOD_MERSENNE(s1);
             
542   }
                                              
543   myuint_t MixMaxRng::fmodmulM61(myuint_t cum,    
544   {
                                              
545     __uint128_t temp;
                            
546     temp = (__uint128_t)a*(__uint128_t)b + cum    
547     return mod128(temp);
                         
548   }
                                              
549 #else // on all other platforms, including 32-    
550   myuint_t MixMaxRng::fmodmulM61(myuint_t cum,    
551   {
                                              
552     const myuint_t MASK32=0xFFFFFFFFULL;
         
553     myuint_t o,ph,pl,ah,al;
                      
554     o=(s)*a;
                                     
555     ph = ((s)>>32);
                              
556     pl = (s) & MASK32;
                           
557     ah = a>>32;
                                  
558     al = a & MASK32;
                             
559     o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*    
560     o += cum;
                                    
561     o = (o & M61) + ((o>>61));
                   
562     return o;
                                    
563   }
                                              
564 #endif
                                           
565 
                                                 
566 void MixMaxRng::print_state() const
              
567 {
                                                
568    std::cout << "mixmax state, file version 1.    
569    std::cout << "N=" << rng_get_N() << "; V[N]    
570    for (int j=0; (j< (rng_get_N()-1) ); ++j) {    
571    std::cout << V[rng_get_N()-1];
                
572    std::cout << "}; ";
                           
573    std::cout << "counter= " << counter;
          
574    std::cout << "sumtot= " << sumtot << "\n";
    
575 }
                                                
576 
                                                 
577 MixMaxRng MixMaxRng::Branch()
                    
578 {
                                                
579    sumtot = iterate_raw_vec(V, sumtot); counte    
580    MixMaxRng tmp=*this;
                          
581    tmp.BranchInplace(0); // daughter id
          
582    return tmp;
                                   
583 }
                                                
584     
                                             
585 void MixMaxRng::BranchInplace(int id)
            
586 {
                                                
587    // Dont forget to iterate the mother, when     
588    // a 64-bit LCG from Knuth line 26, is used    
589    constexpr myuint_t MULT64=63641362238467930    
590    myuint_t tmp=V[id];
                           
591    V[1] *= MULT64; V[id] &= M61;
                 
592    sumtot = MIXMAX_MOD_MERSENNE( sumtot + V[id    
593    sumtot = iterate_raw_vec(V, sumtot);
          
594    counter = 1;
                                  
595 }
                                                
596 
                                                 
597 }  // namespace CLHEP
                            
598