37 #ifndef VIGRA_RANDOM_HXX
38 #define VIGRA_RANDOM_HXX
40 #include "mathutil.hxx"
41 #include "functortraits.hxx"
47 enum RandomSeedTag { RandomSeed };
51 enum RandomEngineTag { TT800, MT19937 };
53 template<RandomEngineTag EngineTag>
56 template <RandomEngineTag EngineTag>
57 void seed(
UInt32 theSeed, RandomState<EngineTag> & engine)
59 engine.state_[0] = theSeed;
60 for(
UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
62 engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
66 template <
class Iterator, RandomEngineTag EngineTag>
67 void seed(Iterator init,
UInt32 key_length, RandomState<EngineTag> & engine)
69 const UInt32 N = RandomState<EngineTag>::N;
70 int k = (int)
std::max(N, key_length);
75 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
81 engine.state_[0] = engine.state_[N-1];
93 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
98 engine.state_[0] = engine.state_[N-1];
103 engine.state_[0] = 0x80000000U;
106 template <RandomEngineTag EngineTag>
107 void seed(RandomSeedTag, RandomState<EngineTag> & engine)
109 static UInt32 globalCount = 0;
111 seed(init, 3, engine);
117 struct RandomState<TT800>
119 static const UInt32 N = 25, M = 7;
128 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
129 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
130 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
131 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
132 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
136 state_[i] = seeds[i];
144 generateNumbers<void>();
146 UInt32 y = state_[current_++];
147 y ^= (y << 7) & 0x2b5b2500;
148 y ^= (y << 15) & 0xdb8b0000;
149 return y ^ (y >> 16);
152 template <
class DUMMY>
153 void generateNumbers()
const;
155 void seedImpl(RandomSeedTag)
157 seed(RandomSeed, *
this);
160 void seedImpl(
UInt32 theSeed)
162 seed(theSeed, *
this);
165 template<
class Iterator>
166 void seedImpl(Iterator init,
UInt32 length)
168 seed(init, length, *
this);
172 template <
class DUMMY>
173 void RandomState<TT800>::generateNumbers()
const
175 UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
177 for(
UInt32 i=0; i<N-M; ++i)
179 state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
181 for (
UInt32 i=N-M; i<N; ++i)
183 state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
190 struct RandomState<MT19937>
192 static const UInt32 N = 624, M = 397;
200 seed(19650218U, *
this);
208 generateNumbers<void>();
210 UInt32 x = state_[current_++];
212 x ^= (x << 7) & 0x9D2C5680U;
213 x ^= (x << 15) & 0xEFC60000U;
214 return x ^ (x >> 18);
217 template <
class DUMMY>
218 void generateNumbers()
const;
222 return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
223 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
226 void seedImpl(RandomSeedTag)
228 seed(RandomSeed, *
this);
229 generateNumbers<void>();
232 void seedImpl(
UInt32 theSeed)
234 seed(theSeed, *
this);
235 generateNumbers<void>();
238 template<
class Iterator>
239 void seedImpl(Iterator init,
UInt32 length)
241 seed(19650218U, *
this);
242 seed(init, length, *
this);
243 generateNumbers<void>();
247 template <
class DUMMY>
248 void RandomState<MT19937>::generateNumbers()
const
250 for (
unsigned int i = 0; i < (N - M); ++i)
252 state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
254 for (
unsigned int i = N - M; i < (N - 1); ++i)
256 state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
258 state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
287 template <
class Engine = detail::RandomState<detail::TT800> >
291 mutable double normalCached_;
292 mutable bool normalCachedValid_;
302 : normalCached_(0.0),
303 normalCachedValid_(false)
317 : normalCached_(0.0),
318 normalCachedValid_(false)
320 this->seedImpl(RandomSeed);
328 : normalCached_(0.0),
329 normalCachedValid_(false)
331 this->seedImpl(theSeed);
339 template<
class Iterator>
341 : normalCached_(0.0),
342 normalCachedValid_(false)
344 this->seedImpl(init, length);
362 this->seedImpl(RandomSeed);
363 normalCachedValid_ =
false;
372 this->seedImpl(theSeed);
373 normalCachedValid_ =
false;
381 template<
class Iterator>
384 this->seedImpl(init, length);
385 normalCachedValid_ =
false;
407 #if 0 // difficult implementation necessary if low bits are not sufficiently random
414 UInt32 factor = factorForUniformInt(beyond);
415 UInt32 res = this->
get() / factor;
420 res = this->
get() / factor;
442 while(res > lastSafeValue)
455 return ( (this->
get() >> 5) * 67108864.0 + (this->
get() >> 6)) * (1.0/9007199254740992.0);
465 return (
double)this->
get() / 4294967295.0;
473 double uniform(
double lower,
double upper)
const
475 vigra_precondition(lower < upper,
476 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
477 return uniform() * (upper-lower) + lower;
491 double normal(
double mean,
double stddev)
const
493 vigra_precondition(stddev > 0.0,
494 "RandomNumberGenerator::normal(): standard deviation must be positive.");
495 return normal()*stddev + mean;
511 return (range > 2147483648U || range == 0)
517 template <
class Engine>
520 if(normalCachedValid_)
522 normalCachedValid_ =
false;
523 return normalCached_;
530 x1 = uniform(-1.0, 1.0);
531 x2 = uniform(-1.0, 1.0);
532 w = x1 * x1 + x2 * x2;
534 while ( w > 1.0 || w == 0.0);
538 normalCached_ = x2 * w;
539 normalCachedValid_ =
true;
569 template <
class Engine>
570 class FunctorTraits<RandomNumberGenerator<Engine> >
573 typedef RandomNumberGenerator<Engine> type;
575 typedef VigraTrueType isInitializer;
577 typedef VigraFalseType isUnaryFunctor;
578 typedef VigraFalseType isBinaryFunctor;
579 typedef VigraFalseType isTernaryFunctor;
581 typedef VigraFalseType isUnaryAnalyser;
582 typedef VigraFalseType isBinaryAnalyser;
583 typedef VigraFalseType isTernaryAnalyser;
600 template <
class Engine = RandomTT800>
603 UInt32 lower_, difference_, factor_;
604 Engine
const & generator_;
618 : lower_(0), difference_(0xffffffff), factor_(1),
619 generator_(generator),
633 Engine
const & generator = Engine::global(),
634 bool useLowBits =
true)
635 : lower_(lower), difference_(upper-lower),
636 factor_(Engine::factorForUniformInt(difference_ + 1)),
637 generator_(generator),
638 useLowBits_(useLowBits)
640 vigra_precondition(lower < upper,
641 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
648 if(difference_ == 0xffffffff)
651 return generator_.uniformInt(difference_+1) + lower_;
654 UInt32 res = generator_() / factor_;
658 while(res > difference_)
659 res = generator_() / factor_;
675 return generator_.uniformInt(beyond);
679 template <
class Engine>
680 class FunctorTraits<UniformIntRandomFunctor<Engine> >
683 typedef UniformIntRandomFunctor<Engine> type;
685 typedef VigraTrueType isInitializer;
687 typedef VigraTrueType isUnaryFunctor;
688 typedef VigraFalseType isBinaryFunctor;
689 typedef VigraFalseType isTernaryFunctor;
691 typedef VigraFalseType isUnaryAnalyser;
692 typedef VigraFalseType isBinaryAnalyser;
693 typedef VigraFalseType isTernaryAnalyser;
707 template <
class Engine = RandomTT800>
710 double offset_, scale_;
711 Engine
const & generator_;
725 generator_(generator)
734 Engine & generator = Engine::global())
736 scale_(upper - lower),
737 generator_(generator)
739 vigra_precondition(lower < upper,
740 "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
747 return generator_.uniform() * scale_ + offset_;
751 template <
class Engine>
752 class FunctorTraits<UniformRandomFunctor<Engine> >
755 typedef UniformRandomFunctor<Engine> type;
757 typedef VigraTrueType isInitializer;
759 typedef VigraFalseType isUnaryFunctor;
760 typedef VigraFalseType isBinaryFunctor;
761 typedef VigraFalseType isTernaryFunctor;
763 typedef VigraFalseType isUnaryAnalyser;
764 typedef VigraFalseType isBinaryAnalyser;
765 typedef VigraFalseType isTernaryAnalyser;
779 template <
class Engine = RandomTT800>
782 double mean_, stddev_;
783 Engine
const & generator_;
797 generator_(generator)
804 Engine & generator = Engine::global())
807 generator_(generator)
809 vigra_precondition(stddev > 0.0,
810 "NormalRandomFunctor(): standard deviation must be positive.");
817 return generator_.normal() * stddev_ + mean_;
822 template <
class Engine>
823 class FunctorTraits<NormalRandomFunctor<Engine> >
826 typedef UniformRandomFunctor<Engine> type;
828 typedef VigraTrueType isInitializer;
830 typedef VigraFalseType isUnaryFunctor;
831 typedef VigraFalseType isBinaryFunctor;
832 typedef VigraFalseType isTernaryFunctor;
834 typedef VigraFalseType isUnaryAnalyser;
835 typedef VigraFalseType isBinaryAnalyser;
836 typedef VigraFalseType isTernaryAnalyser;
843 #endif // VIGRA_RANDOM_HXX