37#ifndef VIGRA_RANDOM_HXX
38#define VIGRA_RANDOM_HXX
40#include "mathutil.hxx"
41#include "functortraits.hxx"
42#include "array_vector.hxx"
49 #include <vigra/windows.h>
54 #include <sys/syscall.h>
60 #include <sys/syscall.h>
61 #include <AvailabilityMacros.h>
66enum RandomSeedTag { RandomSeed };
70enum RandomEngineTag { TT800, MT19937 };
73template<RandomEngineTag EngineTag>
76template <RandomEngineTag EngineTag>
77void seed(
UInt32 theSeed, RandomState<EngineTag> & engine)
79 engine.state_[0] = theSeed;
80 for(
UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
82 engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
86template <
class Iterator, RandomEngineTag EngineTag>
87void seed(Iterator init,
UInt32 key_length, RandomState<EngineTag> & engine)
89 const UInt32 N = RandomState<EngineTag>::N;
90 int k =
static_cast<int>(std::max(N, key_length));
95 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
101 engine.state_[0] = engine.state_[N-1];
113 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
118 engine.state_[0] = engine.state_[N-1];
123 engine.state_[0] = 0x80000000U;
126template <RandomEngineTag EngineTag>
127void seed(RandomSeedTag, RandomState<EngineTag> & engine)
129 static UInt32 globalCount = 0;
130 ArrayVector<UInt32> seedData;
132 seedData.push_back(
static_cast<UInt32>(time(0)));
133 seedData.push_back(
static_cast<UInt32>(clock()));
134 seedData.push_back(++globalCount);
136 std::size_t ptr((
char*)&engine - (
char*)0);
137 seedData.push_back(
static_cast<UInt32>((ptr & 0xffffffff)));
138 static const UInt32 shift =
sizeof(ptr) > 4 ? 32 : 16;
139 seedData.push_back(
static_cast<UInt32>((ptr >> shift)));
142 seedData.push_back(
static_cast<UInt32>(GetCurrentProcessId()));
143 seedData.push_back(
static_cast<UInt32>(GetCurrentThreadId()));
147 seedData.push_back(
static_cast<UInt32>(getpid()));
149 seedData.push_back(
static_cast<UInt32>(syscall(SYS_gettid)));
154 seedData.push_back(
static_cast<UInt32>(getpid()));
155 #if __MAC_OS_X_VERSION_MAX_ALLOWED >= MAC_OS_X_VERSION_10_12
157 pthread_threadid_np(NULL, &tid64);
158 seedData.push_back(
static_cast<UInt32>(tid64));
159 #elif defined(SYS_thread_selfid) && (__MAC_OS_X_VERSION_MAX_ALLOWED >= MAC_OS_X_VERSION_10_6)
161 seedData.push_back(
static_cast<UInt32>(syscall(SYS_thread_selfid)));
162 #elif defined(SYS_gettid)
163 seedData.push_back(
static_cast<UInt32>(syscall(SYS_gettid)));
167 seed(seedData.begin(), seedData.size(), engine);
172struct RandomState<TT800>
174 static const UInt32 N = 25, M = 7;
183 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
184 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
185 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
186 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
187 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
191 state_[i] = seeds[i];
199 generateNumbers<void>();
201 UInt32 y = state_[current_++];
202 y ^= (y << 7) & 0x2b5b2500;
203 y ^= (y << 15) & 0xdb8b0000;
204 return y ^ (y >> 16);
207 template <
class DUMMY>
208 void generateNumbers()
const;
210 void seedImpl(RandomSeedTag)
212 seed(RandomSeed, *
this);
215 void seedImpl(
UInt32 theSeed)
217 seed(theSeed, *
this);
220 template<
class Iterator>
221 void seedImpl(Iterator init,
UInt32 length)
223 seed(init, length, *
this);
227template <
class DUMMY>
228void RandomState<TT800>::generateNumbers()
const
230 UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
232 for(
UInt32 i=0; i<N-M; ++i)
234 state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
236 for (
UInt32 i=N-M; i<N; ++i)
238 state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
245struct RandomState<MT19937>
247 static const UInt32 N = 624, M = 397;
255 seed(19650218U, *
this);
263 generateNumbers<void>();
265 UInt32 x = state_[current_++];
267 x ^= (x << 7) & 0x9D2C5680U;
268 x ^= (x << 15) & 0xEFC60000U;
269 return x ^ (x >> 18);
272 template <
class DUMMY>
273 void generateNumbers()
const;
277 return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
278 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
281 void seedImpl(RandomSeedTag)
283 seed(RandomSeed, *
this);
284 generateNumbers<void>();
287 void seedImpl(
UInt32 theSeed)
289 seed(theSeed, *
this);
290 generateNumbers<void>();
293 template<
class Iterator>
294 void seedImpl(Iterator init,
UInt32 length)
296 seed(19650218U, *
this);
297 seed(init, length, *
this);
298 generateNumbers<void>();
302template <
class DUMMY>
303void RandomState<MT19937>::generateNumbers()
const
305 for (
unsigned int i = 0; i < (N - M); ++i)
307 state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
309 for (
unsigned int i = N - M; i < (N - 1); ++i)
311 state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
313 state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
342template <
class Engine = detail::RandomState<detail::MT19937> >
346 mutable double normalCached_;
347 mutable bool normalCachedValid_;
357 : normalCached_(0.0),
358 normalCachedValid_(false)
372 : normalCached_(0.0),
373 normalCachedValid_(false)
375 this->seedImpl(RandomSeed);
388 : normalCached_(0.0),
389 normalCachedValid_(false)
392 this->seedImpl(RandomSeed);
394 this->seedImpl(theSeed);
402 template<
class Iterator>
404 : normalCached_(0.0),
405 normalCachedValid_(false)
407 this->seedImpl(init, length);
425 this->seedImpl(RandomSeed);
426 normalCachedValid_ =
false;
440 this->seedImpl(RandomSeed);
442 this->seedImpl(theSeed);
443 normalCachedValid_ =
false;
451 template<
class Iterator>
454 this->seedImpl(init, length);
455 normalCachedValid_ =
false;
484 UInt32 factor = factorForUniformInt(beyond);
485 UInt32 res = this->get() / factor;
490 res = this->get() / factor;
506 UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
507 UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
512 while(res > lastSafeValue)
525 return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0);
535 return static_cast<double>(this->get()) / 4294967295.0;
543 double uniform(
double lower,
double upper)
const
545 vigra_precondition(lower < upper,
546 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
547 return uniform() * (upper-lower) + lower;
561 double normal(
double mean,
double stddev)
const
563 vigra_precondition(stddev > 0.0,
564 "RandomNumberGenerator::normal(): standard deviation must be positive.");
565 return normal()*stddev + mean;
580 return (range > 2147483648U || range == 0)
588template <
class Engine>
592template <
class Engine>
595 if(normalCachedValid_)
597 normalCachedValid_ =
false;
598 return normalCached_;
607 w = x1 * x1 + x2 * x2;
609 while ( w > 1.0 || w == 0.0);
611 w = std::sqrt( -2.0 * std::log( w ) / w );
613 normalCached_ = x2 * w;
614 normalCachedValid_ =
true;
644template <
class Engine>
645class FunctorTraits<RandomNumberGenerator<Engine> >
648 typedef RandomNumberGenerator<Engine> type;
650 typedef VigraTrueType isInitializer;
652 typedef VigraFalseType isUnaryFunctor;
653 typedef VigraFalseType isBinaryFunctor;
654 typedef VigraFalseType isTernaryFunctor;
656 typedef VigraFalseType isUnaryAnalyser;
657 typedef VigraFalseType isBinaryAnalyser;
658 typedef VigraFalseType isTernaryAnalyser;
675template <
class Engine = MersenneTwister>
678 UInt32 lower_, difference_, factor_;
679 Engine
const & generator_;
693 : lower_(0), difference_(0xffffffff), factor_(1),
694 generator_(generator),
708 Engine
const & generator = Engine::global(),
709 bool useLowBits =
true)
710 : lower_(lower), difference_(upper-lower),
711 factor_(Engine::factorForUniformInt(difference_ + 1)),
712 generator_(generator),
713 useLowBits_(useLowBits)
715 vigra_precondition(lower < upper,
716 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
723 if(difference_ == 0xffffffff)
726 return generator_.uniformInt(difference_+1) + lower_;
729 UInt32 res = generator_() / factor_;
733 while(res > difference_)
734 res = generator_() / factor_;
750 return generator_.uniformInt(beyond);
754template <
class Engine>
755class FunctorTraits<UniformIntRandomFunctor<Engine> >
758 typedef UniformIntRandomFunctor<Engine> type;
760 typedef VigraTrueType isInitializer;
762 typedef VigraTrueType isUnaryFunctor;
763 typedef VigraFalseType isBinaryFunctor;
764 typedef VigraFalseType isTernaryFunctor;
766 typedef VigraFalseType isUnaryAnalyser;
767 typedef VigraFalseType isBinaryAnalyser;
768 typedef VigraFalseType isTernaryAnalyser;
782template <
class Engine = MersenneTwister>
785 double offset_, scale_;
786 Engine
const & generator_;
800 generator_(generator)
809 Engine & generator = Engine::global())
811 scale_(upper - lower),
812 generator_(generator)
814 vigra_precondition(lower < upper,
815 "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
822 return generator_.uniform() * scale_ + offset_;
826template <
class Engine>
827class FunctorTraits<UniformRandomFunctor<Engine> >
830 typedef UniformRandomFunctor<Engine> type;
832 typedef VigraTrueType isInitializer;
834 typedef VigraFalseType isUnaryFunctor;
835 typedef VigraFalseType isBinaryFunctor;
836 typedef VigraFalseType isTernaryFunctor;
838 typedef VigraFalseType isUnaryAnalyser;
839 typedef VigraFalseType isBinaryAnalyser;
840 typedef VigraFalseType isTernaryAnalyser;
854template <
class Engine = MersenneTwister>
857 double mean_, stddev_;
858 Engine
const & generator_;
872 generator_(generator)
879 Engine & generator = Engine::global())
882 generator_(generator)
884 vigra_precondition(stddev > 0.0,
885 "NormalRandomFunctor(): standard deviation must be positive.");
892 return generator_.normal() * stddev_ + mean_;
897template <
class Engine>
898class FunctorTraits<NormalRandomFunctor<Engine> >
901 typedef UniformRandomFunctor<Engine> type;
903 typedef VigraTrueType isInitializer;
905 typedef VigraFalseType isUnaryFunctor;
906 typedef VigraFalseType isBinaryFunctor;
907 typedef VigraFalseType isTernaryFunctor;
909 typedef VigraFalseType isUnaryAnalyser;
910 typedef VigraFalseType isBinaryAnalyser;
911 typedef VigraFalseType isTernaryAnalyser;
NormalRandomFunctor(Engine const &generator=Engine::global())
Definition random.hxx:869
double result_type
STL required functor result type.
Definition random.hxx:862
NormalRandomFunctor(double mean, double stddev, Engine &generator=Engine::global())
Definition random.hxx:878
double operator()() const
Definition random.hxx:890
Definition random.hxx:345
RandomNumberGenerator(Iterator init, UInt32 length)
Definition random.hxx:403
RandomNumberGenerator()
Definition random.hxx:356
double uniform(double lower, double upper) const
Definition random.hxx:543
void seed(Iterator init, UInt32 length)
Definition random.hxx:452
UInt32 operator()() const
Definition random.hxx:462
void seed(RandomSeedTag)
Definition random.hxx:423
UInt32 uniformInt() const
Definition random.hxx:471
UInt32 uniformInt(UInt32 beyond) const
Definition random.hxx:499
void seed(UInt32 theSeed, bool ignoreSeed=false)
Definition random.hxx:437
RandomNumberGenerator(RandomSeedTag)
Definition random.hxx:371
double normal() const
Definition random.hxx:593
double uniform53() const
Definition random.hxx:522
double normal(double mean, double stddev) const
Definition random.hxx:561
RandomNumberGenerator(UInt32 theSeed, bool ignoreSeed=false)
Definition random.hxx:387
static RandomNumberGenerator & global()
Definition random.hxx:573
double uniform() const
Definition random.hxx:533
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:626
RandomNumberGenerator< detail::RandomState< detail::MT19937 > > MersenneTwister
Definition random.hxx:634
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:630
RandomTT800 & randomTT800()
Definition random.hxx:638
RandomNumberGenerator< detail::RandomState< detail::TT800 > > RandomTT800
Definition random.hxx:622
RandomMT19937 & randomMT19937()
Definition random.hxx:642