37#ifndef VIGRA_RANDOM_HXX
38#define VIGRA_RANDOM_HXX
40#include "mathutil.hxx"
41#include "functortraits.hxx"
42#include "array_vector.hxx"
50 #include <vigra/windows.h>
55 #include <sys/syscall.h>
61 #include <sys/syscall.h>
62 #include <AvailabilityMacros.h>
67enum RandomSeedTag { RandomSeed };
71enum RandomEngineTag { TT800, MT19937 };
74template<RandomEngineTag EngineTag>
77template <RandomEngineTag EngineTag>
78void seed(
UInt32 theSeed, RandomState<EngineTag> & engine)
80 engine.state_[0] = theSeed;
81 for(
UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
83 engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
87template <
class Iterator, RandomEngineTag EngineTag>
88void seed(Iterator init,
UInt32 key_length, RandomState<EngineTag> & engine)
90 const UInt32 N = RandomState<EngineTag>::N;
91 int k =
static_cast<int>(std::max(N, key_length));
96 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
102 engine.state_[0] = engine.state_[N-1];
114 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
119 engine.state_[0] = engine.state_[N-1];
124 engine.state_[0] = 0x80000000U;
127template <RandomEngineTag EngineTag>
128void seed(RandomSeedTag, RandomState<EngineTag> & engine)
130 static UInt32 globalCount = 0;
131 ArrayVector<UInt32> seedData;
133 seedData.push_back(
static_cast<UInt32>(time(0)));
134 seedData.push_back(
static_cast<UInt32>(clock()));
135 seedData.push_back(++globalCount);
137 const uintptr_t ptr =
reinterpret_cast<uintptr_t
>(&engine);
138 seedData.push_back(
static_cast<UInt32>((ptr & 0xffffffff)));
139 static const UInt32 shift =
sizeof(ptr) > 4 ? 32 : 16;
140 seedData.push_back(
static_cast<UInt32>((ptr >> shift)));
143 seedData.push_back(
static_cast<UInt32>(GetCurrentProcessId()));
144 seedData.push_back(
static_cast<UInt32>(GetCurrentThreadId()));
148 seedData.push_back(
static_cast<UInt32>(getpid()));
150 seedData.push_back(
static_cast<UInt32>(syscall(SYS_gettid)));
155 seedData.push_back(
static_cast<UInt32>(getpid()));
156 #if __MAC_OS_X_VERSION_MAX_ALLOWED >= MAC_OS_X_VERSION_10_12
158 pthread_threadid_np(NULL, &tid64);
159 seedData.push_back(
static_cast<UInt32>(tid64));
160 #elif defined(SYS_thread_selfid) && (__MAC_OS_X_VERSION_MAX_ALLOWED >= MAC_OS_X_VERSION_10_6)
162 seedData.push_back(
static_cast<UInt32>(syscall(SYS_thread_selfid)));
163 #elif defined(SYS_gettid)
164 seedData.push_back(
static_cast<UInt32>(syscall(SYS_gettid)));
168 seed(seedData.begin(), seedData.size(), engine);
173struct RandomState<TT800>
175 static const UInt32 N = 25, M = 7;
184 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
185 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
186 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
187 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
188 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
192 state_[i] = seeds[i];
200 generateNumbers<void>();
202 UInt32 y = state_[current_++];
203 y ^= (y << 7) & 0x2b5b2500;
204 y ^= (y << 15) & 0xdb8b0000;
205 return y ^ (y >> 16);
208 template <
class DUMMY>
209 void generateNumbers()
const;
211 void seedImpl(RandomSeedTag)
213 seed(RandomSeed, *
this);
216 void seedImpl(
UInt32 theSeed)
218 seed(theSeed, *
this);
221 template<
class Iterator>
222 void seedImpl(Iterator init,
UInt32 length)
224 seed(init, length, *
this);
228template <
class DUMMY>
229void RandomState<TT800>::generateNumbers()
const
231 UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
233 for(
UInt32 i=0; i<N-M; ++i)
235 state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
237 for (
UInt32 i=N-M; i<N; ++i)
239 state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
246struct RandomState<MT19937>
248 static const UInt32 N = 624, M = 397;
256 seed(19650218U, *
this);
264 generateNumbers<void>();
266 UInt32 x = state_[current_++];
268 x ^= (x << 7) & 0x9D2C5680U;
269 x ^= (x << 15) & 0xEFC60000U;
270 return x ^ (x >> 18);
273 template <
class DUMMY>
274 void generateNumbers()
const;
278 return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
279 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
282 void seedImpl(RandomSeedTag)
284 seed(RandomSeed, *
this);
285 generateNumbers<void>();
288 void seedImpl(
UInt32 theSeed)
290 seed(theSeed, *
this);
291 generateNumbers<void>();
294 template<
class Iterator>
295 void seedImpl(Iterator init,
UInt32 length)
297 seed(19650218U, *
this);
298 seed(init, length, *
this);
299 generateNumbers<void>();
303template <
class DUMMY>
304void RandomState<MT19937>::generateNumbers()
const
306 for (
unsigned int i = 0; i < (N - M); ++i)
308 state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
310 for (
unsigned int i = N - M; i < (N - 1); ++i)
312 state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
314 state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
343template <
class Engine = detail::RandomState<detail::MT19937> >
347 mutable double normalCached_;
348 mutable bool normalCachedValid_;
358 : normalCached_(0.0),
359 normalCachedValid_(
false)
373 : normalCached_(0.0),
374 normalCachedValid_(
false)
376 this->seedImpl(RandomSeed);
389 : normalCached_(0.0),
390 normalCachedValid_(
false)
393 this->seedImpl(RandomSeed);
403 template<
class Iterator>
405 : normalCached_(0.0),
406 normalCachedValid_(
false)
408 this->seedImpl(init, length);
426 this->seedImpl(RandomSeed);
427 normalCachedValid_ =
false;
441 this->seedImpl(RandomSeed);
444 normalCachedValid_ =
false;
452 template<
class Iterator>
455 this->seedImpl(init, length);
456 normalCachedValid_ =
false;
491 res = this->get() /
factor;
526 return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0);
536 return static_cast<double>(this->get()) / 4294967295.0;
544 double uniform(
double lower,
double upper)
const
546 vigra_precondition(lower < upper,
547 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
548 return uniform() * (upper-lower) + lower;
564 vigra_precondition(
stddev > 0.0,
565 "RandomNumberGenerator::normal(): standard deviation must be positive.");
581 return (range > 2147483648U || range == 0)
589template <
class Engine>
590RandomNumberGenerator<Engine> RandomNumberGenerator<Engine>::global_(RandomSeed);
593template <
class Engine>
596 if(normalCachedValid_)
598 normalCachedValid_ =
false;
599 return normalCached_;
606 x1 = uniform(-1.0, 1.0);
607 x2 = uniform(-1.0, 1.0);
610 while ( w > 1.0 || w == 0.0);
612 w = std::sqrt( -2.0 * std::log( w ) / w );
614 normalCached_ =
x2 * w;
615 normalCachedValid_ =
true;
645template <
class Engine>
646class FunctorTraits<RandomNumberGenerator<Engine> >
649 typedef RandomNumberGenerator<Engine> type;
651 typedef VigraTrueType isInitializer;
653 typedef VigraFalseType isUnaryFunctor;
654 typedef VigraFalseType isBinaryFunctor;
655 typedef VigraFalseType isTernaryFunctor;
657 typedef VigraFalseType isUnaryAnalyser;
658 typedef VigraFalseType isBinaryAnalyser;
659 typedef VigraFalseType isTernaryAnalyser;
676template <
class Engine = MersenneTwister>
679 UInt32 lower_, difference_, factor_;
680 Engine const & generator_;
694 : lower_(0), difference_(0xffffffff), factor_(1),
711 : lower_(lower), difference_(upper-lower),
712 factor_(
Engine::factorForUniformInt(difference_ + 1)),
716 vigra_precondition(lower < upper,
717 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
724 if(difference_ == 0xffffffff)
727 return generator_.uniformInt(difference_+1) + lower_;
730 UInt32 res = generator_() / factor_;
734 while(res > difference_)
735 res = generator_() / factor_;
751 return generator_.uniformInt(
beyond);
755template <
class Engine>
756class FunctorTraits<UniformIntRandomFunctor<Engine> >
759 typedef UniformIntRandomFunctor<Engine> type;
761 typedef VigraTrueType isInitializer;
763 typedef VigraTrueType isUnaryFunctor;
764 typedef VigraFalseType isBinaryFunctor;
765 typedef VigraFalseType isTernaryFunctor;
767 typedef VigraFalseType isUnaryAnalyser;
768 typedef VigraFalseType isBinaryAnalyser;
769 typedef VigraFalseType isTernaryAnalyser;
783template <
class Engine = MersenneTwister>
786 double offset_, scale_;
787 Engine const & generator_;
812 scale_(upper - lower),
815 vigra_precondition(lower < upper,
816 "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
823 return generator_.uniform() * scale_ + offset_;
827template <
class Engine>
828class FunctorTraits<UniformRandomFunctor<Engine> >
831 typedef UniformRandomFunctor<Engine> type;
833 typedef VigraTrueType isInitializer;
835 typedef VigraFalseType isUnaryFunctor;
836 typedef VigraFalseType isBinaryFunctor;
837 typedef VigraFalseType isTernaryFunctor;
839 typedef VigraFalseType isUnaryAnalyser;
840 typedef VigraFalseType isBinaryAnalyser;
841 typedef VigraFalseType isTernaryAnalyser;
855template <
class Engine = MersenneTwister>
858 double mean_, stddev_;
859 Engine const & generator_;
885 vigra_precondition(
stddev > 0.0,
886 "NormalRandomFunctor(): standard deviation must be positive.");
893 return generator_.normal() * stddev_ + mean_;
898template <
class Engine>
899class FunctorTraits<NormalRandomFunctor<Engine> >
902 typedef UniformRandomFunctor<Engine> type;
904 typedef VigraTrueType isInitializer;
906 typedef VigraFalseType isUnaryFunctor;
907 typedef VigraFalseType isBinaryFunctor;
908 typedef VigraFalseType isTernaryFunctor;
910 typedef VigraFalseType isUnaryAnalyser;
911 typedef VigraFalseType isBinaryAnalyser;
912 typedef VigraFalseType isTernaryAnalyser;
Definition random.hxx:857
NormalRandomFunctor(Engine const &generator=Engine::global())
Definition random.hxx:870
double result_type
STL required functor result type.
Definition random.hxx:863
NormalRandomFunctor(double mean, double stddev, Engine &generator=Engine::global())
Definition random.hxx:879
double operator()() const
Definition random.hxx:891
Class for a single RGB value.
Definition rgbvalue.hxx:128
RGBValue()
Definition rgbvalue.hxx:209
Definition random.hxx:346
RandomNumberGenerator(Iterator init, UInt32 length)
Definition random.hxx:404
RandomNumberGenerator()
Definition random.hxx:357
double uniform(double lower, double upper) const
Definition random.hxx:544
void seed(Iterator init, UInt32 length)
Definition random.hxx:453
UInt32 operator()() const
Definition random.hxx:463
void seed(RandomSeedTag)
Definition random.hxx:424
UInt32 uniformInt() const
Definition random.hxx:472
UInt32 uniformInt(UInt32 beyond) const
Definition random.hxx:500
void seed(UInt32 theSeed, bool ignoreSeed=false)
Definition random.hxx:438
RandomNumberGenerator(RandomSeedTag)
Definition random.hxx:372
double normal() const
Definition random.hxx:594
double uniform53() const
Definition random.hxx:523
double normal(double mean, double stddev) const
Definition random.hxx:562
RandomNumberGenerator(UInt32 theSeed, bool ignoreSeed=false)
Definition random.hxx:388
static RandomNumberGenerator & global()
Definition random.hxx:574
double uniform() const
Definition random.hxx:534
LookupTag< TAG, A >::result_type get(A const &a)
Definition accumulator.hxx:2942
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition sized_int.hxx:183
RandomNumberGenerator< detail::RandomState< detail::TT800 > > TemperedTwister
Definition random.hxx:627
RandomNumberGenerator< detail::RandomState< detail::MT19937 > > MersenneTwister
Definition random.hxx:635
UInt32 ceilPower2(UInt32 x)
Round up to the nearest power of 2.
Definition mathutil.hxx:294
RandomNumberGenerator< detail::RandomState< detail::MT19937 > > RandomMT19937
Definition random.hxx:631
RandomTT800 & randomTT800()
Definition random.hxx:639
RandomNumberGenerator< detail::RandomState< detail::TT800 > > RandomTT800
Definition random.hxx:623
RandomMT19937 & randomMT19937()
Definition random.hxx:643