| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms. 2 // See the file tools.license for terms.
3 3
4 #ifndef tools_rtausmeui 4 #ifndef tools_rtausmeui
5 #define tools_rtausmeui 5 #define tools_rtausmeui
6 6
7 // G.Barrand : not so clear if 0 and tools::ui 7 // G.Barrand : not so clear if 0 and tools::uint32_max() are included.
8 // A simple program shows that "if 8 // A simple program shows that "if(r.shoot()==tools::uint32_max())" shows hits.
9 // (But we did not see hits for "i 9 // (But we did not see hits for "if(r.shoot()==0)"). (See tools/tests/rand.cpp).
10 10
11 // tausme is for Tausworthe maxmally equidistr 11 // tausme is for Tausworthe maxmally equidistributed.
12 12
13 // From logic and code of CERN-ROOT/TRandom2 c 13 // From logic and code of CERN-ROOT/TRandom2 class.
14 14
15 // Random number generator class based on the 15 // Random number generator class based on the maximally quidistributed combined
16 // Tausworthe generator by L'Ecuyer. 16 // Tausworthe generator by L'Ecuyer.
17 // 17 //
18 // The period of the generator is 2**88 (about 18 // The period of the generator is 2**88 (about 10**26) and it uses only 3 words
19 // for the state. 19 // for the state.
20 // 20 //
21 // For more information see: 21 // For more information see:
22 // P. L'Ecuyer, Mathematics of Computation, 65 22 // P. L'Ecuyer, Mathematics of Computation, 65, 213 (1996)
23 // P. L'Ecuyer, Mathematics of Computation, 68 23 // P. L'Ecuyer, Mathematics of Computation, 68, 225 (1999)
24 // 24 //
25 // The publication are available online at 25 // The publication are available online at
26 // http://www.iro.umontreal.ca/~lecuyer/myftp 26 // http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
27 // http://www.iro.umontreal.ca/~lecuyer/myftp 27 // http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
28 28
29 #ifdef TOOLS_MEM 29 #ifdef TOOLS_MEM
30 #include "mem" 30 #include "mem"
31 #include "S_STRING" 31 #include "S_STRING"
32 #endif 32 #endif
33 33
34 namespace tools { 34 namespace tools {
35 35
36 class rtausmeui { 36 class rtausmeui {
37 #ifdef TOOLS_MEM 37 #ifdef TOOLS_MEM
38 TOOLS_SCLASS(tools::rtausmeui) 38 TOOLS_SCLASS(tools::rtausmeui)
39 #endif 39 #endif
40 public: 40 public:
41 rtausmeui(unsigned int a_seed = 1):m_seed(0) 41 rtausmeui(unsigned int a_seed = 1):m_seed(0),m_seed1(0),m_seed2(0){
42 #ifdef TOOLS_MEM 42 #ifdef TOOLS_MEM
43 mem::increment(s_class().c_str()); 43 mem::increment(s_class().c_str());
44 #endif 44 #endif
45 set_seed(a_seed); 45 set_seed(a_seed);
46 } 46 }
47 virtual ~rtausmeui(){ 47 virtual ~rtausmeui(){
48 #ifdef TOOLS_MEM 48 #ifdef TOOLS_MEM
49 mem::decrement(s_class().c_str()); 49 mem::decrement(s_class().c_str());
50 #endif 50 #endif
51 } 51 }
52 public: 52 public:
53 rtausmeui(const rtausmeui& a_from):m_seed(a_ 53 rtausmeui(const rtausmeui& a_from):m_seed(a_from.m_seed),m_seed1(a_from.m_seed1),m_seed2(a_from.m_seed2){
54 #ifdef TOOLS_MEM 54 #ifdef TOOLS_MEM
55 mem::increment(s_class().c_str()); 55 mem::increment(s_class().c_str());
56 #endif 56 #endif
57 } 57 }
58 rtausmeui& operator=(const rtausmeui& a_from 58 rtausmeui& operator=(const rtausmeui& a_from) {
59 m_seed = a_from.m_seed; 59 m_seed = a_from.m_seed;
60 m_seed1 = a_from.m_seed1; 60 m_seed1 = a_from.m_seed1;
61 m_seed2 = a_from.m_seed2; 61 m_seed2 = a_from.m_seed2;
62 return *this; 62 return *this;
63 } 63 }
64 public: 64 public:
65 void set_seed(unsigned int a_seed) { 65 void set_seed(unsigned int a_seed) {
66 m_seed = a_seed?a_seed:1; 66 m_seed = a_seed?a_seed:1;
67 67
68 // Generate m_seed[1,2] needed for the gen 68 // Generate m_seed[1,2] needed for the generator state using
69 // a linear congruential generator 69 // a linear congruential generator
70 // The only condition, stated at the end o 70 // The only condition, stated at the end of the 1999 L'Ecuyer paper is that the seeds
71 // must be greater than 1,7 and 15. 71 // must be greater than 1,7 and 15.
72 72
73 m_seed = LCG(m_seed); 73 m_seed = LCG(m_seed);
74 if (m_seed < 2) m_seed += 2UL; 74 if (m_seed < 2) m_seed += 2UL;
75 m_seed1 = LCG(m_seed); 75 m_seed1 = LCG(m_seed);
76 if (m_seed1 < 8) m_seed1 += 8UL; 76 if (m_seed1 < 8) m_seed1 += 8UL;
77 m_seed2 = LCG(m_seed1); 77 m_seed2 = LCG(m_seed1);
78 if (m_seed2 < 16) m_seed2 += 16UL; 78 if (m_seed2 < 16) m_seed2 += 16UL;
79 79
80 // "warm it up" by calling it 6 times 80 // "warm it up" by calling it 6 times
81 for (unsigned int i = 0; i < 6; ++i) shoot 81 for (unsigned int i = 0; i < 6; ++i) shoot();
82 } 82 }
83 unsigned int seed() const {return m_seed;} 83 unsigned int seed() const {return m_seed;}
84 84
85 unsigned int shoot() { 85 unsigned int shoot() {
86 // TausWorth generator from L'Ecuyer, use 86 // TausWorth generator from L'Ecuyer, uses as seed 3x32bits integers
87 // Use a mask of 0xffffffffUL to make in 87 // Use a mask of 0xffffffffUL to make in work on 64 bit machines
88 // Periodicity of about 10**26 88 // Periodicity of about 10**26
89 89
90 unsigned int y; 90 unsigned int y;
91 91
92 do { 92 do {
93 m_seed = TAUSWORTHE (m_seed, 13, 19, 42 93 m_seed = TAUSWORTHE (m_seed, 13, 19, 4294967294UL, 12);
94 m_seed1 = TAUSWORTHE (m_seed1, 2, 25, 42 94 m_seed1 = TAUSWORTHE (m_seed1, 2, 25, 4294967288UL, 4);
95 m_seed2 = TAUSWORTHE (m_seed2, 3, 11, 42 95 m_seed2 = TAUSWORTHE (m_seed2, 3, 11, 4294967280UL, 17);
96 96
97 y = m_seed ^ m_seed1 ^ m_seed2; 97 y = m_seed ^ m_seed1 ^ m_seed2;
98 } while(!y); 98 } while(!y);
99 99
100 return y; 100 return y;
101 } 101 }
102 protected: 102 protected:
103 static unsigned int LCG(unsigned int a_n) { 103 static unsigned int LCG(unsigned int a_n) {
104 return ((69069 * a_n) & 0xffffffffUL); // 104 return ((69069 * a_n) & 0xffffffffUL); // linear congurential generator
105 } 105 }
106 static unsigned int TAUSWORTHE(unsigned int 106 static unsigned int TAUSWORTHE(unsigned int a_s,unsigned int a_a,unsigned int a_b,unsigned int a_c,unsigned int a_d) {
107 return (((a_s & a_c) << a_d) & 0xffffffffU 107 return (((a_s & a_c) << a_d) & 0xffffffffUL ) ^ ((((a_s << a_a) & 0xffffffffUL )^ a_s) >> a_b);
108 } 108 }
109 protected: 109 protected:
110 unsigned int m_seed; 110 unsigned int m_seed;
111 unsigned int m_seed1; 111 unsigned int m_seed1;
112 unsigned int m_seed2; 112 unsigned int m_seed2;
113 }; 113 };
114 114
115 } 115 }
116 116
117 #endif 117 #endif