Gabriele Cosmo
101 lines
3.0 KiB
Plaintext
101 lines
3.0 KiB
Plaintext
// Copyright (C) 2010, Guy Barrand. All rights reserved.
|
|
// See the file tools.license for terms.
|
|
|
|
#ifndef tools_rtausmeui
|
|
#define tools_rtausmeui
|
|
|
|
// G.Barrand : not so clear if 0 and tools::uint32_max() are included.
|
|
// A simple program shows that "if(r.shoot()==tools::uint32_max())" shows hits.
|
|
// (But we did not see hits for "if(r.shoot()==0)"). (See tools/tests/rand.cpp).
|
|
|
|
// tausme is for Tausworthe maxmally equidistributed.
|
|
|
|
// From logic and code of CERN-ROOT/TRandom2 class.
|
|
|
|
// Random number generator class based on the maximally quidistributed combined
|
|
// Tausworthe generator by L'Ecuyer.
|
|
//
|
|
// The period of the generator is 2**88 (about 10**26) and it uses only 3 words
|
|
// for the state.
|
|
//
|
|
// For more information see:
|
|
// P. L'Ecuyer, Mathematics of Computation, 65, 213 (1996)
|
|
// P. L'Ecuyer, Mathematics of Computation, 68, 225 (1999)
|
|
//
|
|
// The publication are available online at
|
|
// http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
|
|
// http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
|
|
|
|
namespace tools {
|
|
|
|
class rtausmeui {
|
|
public:
|
|
rtausmeui(unsigned int a_seed = 1):m_seed(0),m_seed1(0),m_seed2(0){
|
|
set_seed(a_seed);
|
|
}
|
|
virtual ~rtausmeui(){
|
|
}
|
|
public:
|
|
rtausmeui(const rtausmeui& a_from):m_seed(a_from.m_seed),m_seed1(a_from.m_seed1),m_seed2(a_from.m_seed2){
|
|
}
|
|
rtausmeui& operator=(const rtausmeui& a_from) {
|
|
m_seed = a_from.m_seed;
|
|
m_seed1 = a_from.m_seed1;
|
|
m_seed2 = a_from.m_seed2;
|
|
return *this;
|
|
}
|
|
public:
|
|
void set_seed(unsigned int a_seed) {
|
|
m_seed = a_seed?a_seed:1;
|
|
|
|
// Generate m_seed[1,2] needed for the generator state using
|
|
// a linear congruential generator
|
|
// The only condition, stated at the end of the 1999 L'Ecuyer paper is that the seeds
|
|
// must be greater than 1,7 and 15.
|
|
|
|
m_seed = LCG(m_seed);
|
|
if (m_seed < 2) m_seed += 2UL;
|
|
m_seed1 = LCG(m_seed);
|
|
if (m_seed1 < 8) m_seed1 += 8UL;
|
|
m_seed2 = LCG(m_seed1);
|
|
if (m_seed2 < 16) m_seed2 += 16UL;
|
|
|
|
// "warm it up" by calling it 6 times
|
|
for (unsigned int i = 0; i < 6; ++i) shoot();
|
|
}
|
|
unsigned int seed() const {return m_seed;}
|
|
|
|
unsigned int shoot() {
|
|
// TausWorth generator from L'Ecuyer, uses as seed 3x32bits integers
|
|
// Use a mask of 0xffffffffUL to make in work on 64 bit machines
|
|
// Periodicity of about 10**26
|
|
|
|
unsigned int y;
|
|
|
|
do {
|
|
m_seed = TAUSWORTHE (m_seed, 13, 19, 4294967294UL, 12);
|
|
m_seed1 = TAUSWORTHE (m_seed1, 2, 25, 4294967288UL, 4);
|
|
m_seed2 = TAUSWORTHE (m_seed2, 3, 11, 4294967280UL, 17);
|
|
|
|
y = m_seed ^ m_seed1 ^ m_seed2;
|
|
} while(!y);
|
|
|
|
return y;
|
|
}
|
|
protected:
|
|
static unsigned int LCG(unsigned int a_n) {
|
|
return ((69069 * a_n) & 0xffffffffUL); // linear congurential generator
|
|
}
|
|
static unsigned int TAUSWORTHE(unsigned int a_s,unsigned int a_a,unsigned int a_b,unsigned int a_c,unsigned int a_d) {
|
|
return (((a_s & a_c) << a_d) & 0xffffffffUL ) ^ ((((a_s << a_a) & 0xffffffffUL )^ a_s) >> a_b);
|
|
}
|
|
protected:
|
|
unsigned int m_seed;
|
|
unsigned int m_seed1;
|
|
unsigned int m_seed2;
|
|
};
|
|
|
|
}
|
|
|
|
#endif
|