ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MersenneTwisterEngine.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: L. Moneta 8/2015
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2015 , ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 // implementation file of Mersenne-Twister engine
12 //
13 //
14 // Created by: Lorenzo Moneta : Tue 4 Aug 2015
15 //
16 //
18 
19 
20 namespace ROOT {
21 namespace Math {
22 
23  /// set the seed x
24  void MersenneTwisterEngine::SetSeed(unsigned int seed) {
25  fCount624 = 624;
26  if (seed > 0) {
27  fMt[0] = seed;
28 
29  // use multipliers from Knuth's "Art of Computer Programming" Vol. 2, 3rd Ed. p.106
30  for(int i=1; i<624; i++) {
31  fMt[i] = (1812433253 * ( fMt[i-1] ^ ( fMt[i-1] >> 30)) + i );
32  }
33  }
34  }
35 
36  /// generate a random double number
38 
39 
40  uint32_t y;
41 
42  const int kM = 397;
43  const int kN = 624;
44  const uint32_t kTemperingMaskB = 0x9d2c5680;
45  const uint32_t kTemperingMaskC = 0xefc60000;
46  const uint32_t kUpperMask = 0x80000000;
47  const uint32_t kLowerMask = 0x7fffffff;
48  const uint32_t kMatrixA = 0x9908b0df;
49 
50  if (fCount624 >= kN) {
51  int i;
52 
53  for (i=0; i < kN-kM; i++) {
54  y = (fMt[i] & kUpperMask) | (fMt[i+1] & kLowerMask);
55  fMt[i] = fMt[i+kM] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
56  }
57 
58  for ( ; i < kN-1 ; i++) {
59  y = (fMt[i] & kUpperMask) | (fMt[i+1] & kLowerMask);
60  fMt[i] = fMt[i+kM-kN] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
61  }
62 
63  y = (fMt[kN-1] & kUpperMask) | (fMt[0] & kLowerMask);
64  fMt[kN-1] = fMt[kM-1] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
65  fCount624 = 0;
66  }
67 
68  y = fMt[fCount624++];
69  y ^= (y >> 11);
70  y ^= ((y << 7 ) & kTemperingMaskB );
71  y ^= ((y << 15) & kTemperingMaskC );
72  y ^= (y >> 18);
73 
74  // 2.3283064365386963e-10 == 1./(max<UINt_t>+1) -> then returned value cannot be = 1.0
75  if (y) return ( (double) y * 2.3283064365386963e-10); // * Power(2,-32)
76  return Rndm_impl();
77 
78  }
79 
80  } // namespace Math
81 } // namespace ROOT
double Rndm_impl()
generate a random double number
void SetSeed(unsigned int seed)
set the seed x
static const double x1[5]
Double_t y[n]
Definition: legend1.C:17