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