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]);
286 template <
class Engine = detail::RandomState<detail::TT800> >
290 mutable double normalCached_;
291 mutable bool normalCachedValid_;
301 : normalCached_(0.0),
302 normalCachedValid_(false)
316 : normalCached_(0.0),
317 normalCachedValid_(false)
319 this->seedImpl(RandomSeed);
327 : normalCached_(0.0),
328 normalCachedValid_(false)
330 this->seedImpl(theSeed);
338 template<
class Iterator>
340 : normalCached_(0.0),
341 normalCachedValid_(false)
343 this->seedImpl(init, length);
361 this->seedImpl(RandomSeed);
362 normalCachedValid_ =
false;
371 this->seedImpl(theSeed);
372 normalCachedValid_ =
false;
380 template<
class Iterator>
383 this->seedImpl(init, length);
384 normalCachedValid_ =
false;
406 #if 0 // difficult implementation necessary if low bits are not sufficiently random
413 UInt32 factor = factorForUniformInt(beyond);
414 UInt32 res = this->
get() / factor;
419 res = this->
get() / factor;
435 UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
436 UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
441 while(res > lastSafeValue)
454 return ( (this->
get() >> 5) * 67108864.0 + (this->
get() >> 6)) * (1.0/9007199254740992.0);
464 return (
double)this->
get() / 4294967295.0;
472 double uniform(
double lower,
double upper)
const
474 vigra_precondition(lower < upper,
475 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
476 return uniform() * (upper-lower) + lower;
490 double normal(
double mean,
double stddev)
const
492 vigra_precondition(stddev > 0.0,
493 "RandomNumberGenerator::normal(): standard deviation must be positive.");
494 return normal()*stddev + mean;
510 return (range > 2147483648U || range == 0)
516 template <
class Engine>
519 if(normalCachedValid_)
521 normalCachedValid_ =
false;
522 return normalCached_;
529 x1 = uniform(-1.0, 1.0);
530 x2 = uniform(-1.0, 1.0);
531 w = x1 * x1 + x2 * x2;
533 while ( w > 1.0 || w == 0.0);
537 normalCached_ = x2 * w;
538 normalCachedValid_ =
true;
568 template <
class Engine>
569 class FunctorTraits<RandomNumberGenerator<Engine> >
572 typedef RandomNumberGenerator<Engine> type;
574 typedef VigraTrueType isInitializer;
576 typedef VigraFalseType isUnaryFunctor;
577 typedef VigraFalseType isBinaryFunctor;
578 typedef VigraFalseType isTernaryFunctor;
580 typedef VigraFalseType isUnaryAnalyser;
581 typedef VigraFalseType isBinaryAnalyser;
582 typedef VigraFalseType isTernaryAnalyser;
597 template <
class Engine = RandomTT800>
600 UInt32 lower_, difference_, factor_;
601 Engine
const & generator_;
615 : lower_(0), difference_(0xffffffff), factor_(1),
616 generator_(generator),
630 Engine
const & generator = Engine::global(),
631 bool useLowBits =
true)
632 : lower_(lower), difference_(upper-lower),
633 factor_(Engine::factorForUniformInt(difference_ + 1)),
634 generator_(generator),
635 useLowBits_(useLowBits)
637 vigra_precondition(lower < upper,
638 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
645 if(difference_ == 0xffffffff)
648 return generator_.uniformInt(difference_+1) + lower_;
651 UInt32 res = generator_() / factor_;
655 while(res > difference_)
656 res = generator_() / factor_;
672 return generator_.uniformInt(beyond);
676 template <
class Engine>
677 class FunctorTraits<UniformIntRandomFunctor<Engine> >
680 typedef UniformIntRandomFunctor<Engine> type;
682 typedef VigraTrueType isInitializer;
684 typedef VigraTrueType isUnaryFunctor;
685 typedef VigraFalseType isBinaryFunctor;
686 typedef VigraFalseType isTernaryFunctor;
688 typedef VigraFalseType isUnaryAnalyser;
689 typedef VigraFalseType isBinaryAnalyser;
690 typedef VigraFalseType isTernaryAnalyser;
703 template <
class Engine = RandomTT800>
706 double offset_, scale_;
707 Engine
const & generator_;
721 generator_(generator)
730 Engine & generator = Engine::global())
732 scale_(upper - lower),
733 generator_(generator)
735 vigra_precondition(lower < upper,
736 "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
743 return generator_.uniform() * scale_ + offset_;
747 template <
class Engine>
748 class FunctorTraits<UniformRandomFunctor<Engine> >
751 typedef UniformRandomFunctor<Engine> type;
753 typedef VigraTrueType isInitializer;
755 typedef VigraFalseType isUnaryFunctor;
756 typedef VigraFalseType isBinaryFunctor;
757 typedef VigraFalseType isTernaryFunctor;
759 typedef VigraFalseType isUnaryAnalyser;
760 typedef VigraFalseType isBinaryAnalyser;
761 typedef VigraFalseType isTernaryAnalyser;
774 template <
class Engine = RandomTT800>
777 double mean_, stddev_;
778 Engine
const & generator_;
792 generator_(generator)
799 Engine & generator = Engine::global())
802 generator_(generator)
804 vigra_precondition(stddev > 0.0,
805 "NormalRandomFunctor(): standard deviation must be positive.");
812 return generator_.normal() * stddev_ + mean_;
817 template <
class Engine>
818 class FunctorTraits<NormalRandomFunctor<Engine> >
821 typedef UniformRandomFunctor<Engine> type;
823 typedef VigraTrueType isInitializer;
825 typedef VigraFalseType isUnaryFunctor;
826 typedef VigraFalseType isBinaryFunctor;
827 typedef VigraFalseType isTernaryFunctor;
829 typedef VigraFalseType isUnaryAnalyser;
830 typedef VigraFalseType isBinaryAnalyser;
831 typedef VigraFalseType isTernaryAnalyser;
838 #endif // VIGRA_RANDOM_HXX