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>
59 #include <sys/syscall.h>
60 #include <AvailabilityMacros.h>
65 enum RandomSeedTag { RandomSeed };
69 enum RandomEngineTag { TT800, MT19937 };
72 template<RandomEngineTag EngineTag>
75 template <RandomEngineTag EngineTag>
76 void seed(
UInt32 theSeed, RandomState<EngineTag> & engine)
78 engine.state_[0] = theSeed;
79 for(
UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
81 engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
85 template <
class Iterator, RandomEngineTag EngineTag>
86 void seed(Iterator init,
UInt32 key_length, RandomState<EngineTag> & engine)
88 const UInt32 N = RandomState<EngineTag>::N;
89 int k = (int)std::max(N, key_length);
94 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
100 engine.state_[0] = engine.state_[N-1];
112 engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
117 engine.state_[0] = engine.state_[N-1];
122 engine.state_[0] = 0x80000000U;
125 template <RandomEngineTag EngineTag>
126 void seed(RandomSeedTag, RandomState<EngineTag> & engine)
128 static UInt32 globalCount = 0;
129 ArrayVector<UInt32> seedData;
131 seedData.push_back((
UInt32)time(0));
132 seedData.push_back((
UInt32)clock());
133 seedData.push_back(++globalCount);
135 std::size_t ptr((
char*)&engine - (
char*)0);
136 seedData.push_back((
UInt32)(ptr & 0xffffffff));
137 seedData.push_back((
UInt32)(ptr >> 32));
140 seedData.push_back((
UInt32)GetCurrentProcessId());
141 seedData.push_back((
UInt32)GetCurrentThreadId());
145 seedData.push_back((
UInt32)getpid());
147 seedData.push_back((
UInt32)syscall(SYS_gettid));
152 seedData.push_back((
UInt32)getpid());
153 #if defined(SYS_thread_selfid) && (MAC_OS_X_VERSION_MIN_REQUIRED >= MAC_OS_X_VERSION_10_6)
155 seedData.push_back((
UInt32)syscall(SYS_thread_selfid));
159 seed(seedData.begin(), seedData.size(), engine);
164 struct RandomState<TT800>
166 static const UInt32 N = 25, M = 7;
175 0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
176 0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
177 0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
178 0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
179 0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
183 state_[i] = seeds[i];
191 generateNumbers<void>();
193 UInt32 y = state_[current_++];
194 y ^= (y << 7) & 0x2b5b2500;
195 y ^= (y << 15) & 0xdb8b0000;
196 return y ^ (y >> 16);
199 template <
class DUMMY>
200 void generateNumbers()
const;
202 void seedImpl(RandomSeedTag)
204 seed(RandomSeed, *
this);
207 void seedImpl(
UInt32 theSeed)
209 seed(theSeed, *
this);
212 template<
class Iterator>
213 void seedImpl(Iterator init,
UInt32 length)
215 seed(init, length, *
this);
219 template <
class DUMMY>
220 void RandomState<TT800>::generateNumbers()
const
222 UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
224 for(
UInt32 i=0; i<N-M; ++i)
226 state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
228 for (
UInt32 i=N-M; i<N; ++i)
230 state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
237 struct RandomState<MT19937>
239 static const UInt32 N = 624, M = 397;
247 seed(19650218U, *
this);
255 generateNumbers<void>();
257 UInt32 x = state_[current_++];
259 x ^= (x << 7) & 0x9D2C5680U;
260 x ^= (x << 15) & 0xEFC60000U;
261 return x ^ (x >> 18);
264 template <
class DUMMY>
265 void generateNumbers()
const;
269 return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
270 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
273 void seedImpl(RandomSeedTag)
275 seed(RandomSeed, *
this);
276 generateNumbers<void>();
279 void seedImpl(
UInt32 theSeed)
281 seed(theSeed, *
this);
282 generateNumbers<void>();
285 template<
class Iterator>
286 void seedImpl(Iterator init,
UInt32 length)
288 seed(19650218U, *
this);
289 seed(init, length, *
this);
290 generateNumbers<void>();
294 template <
class DUMMY>
295 void RandomState<MT19937>::generateNumbers()
const
297 for (
unsigned int i = 0; i < (N - M); ++i)
299 state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
301 for (
unsigned int i = N - M; i < (N - 1); ++i)
303 state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
305 state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
334 template <
class Engine = detail::RandomState<detail::TT800> >
338 mutable double normalCached_;
339 mutable bool normalCachedValid_;
349 : normalCached_(0.0),
350 normalCachedValid_(false)
364 : normalCached_(0.0),
365 normalCachedValid_(false)
367 this->seedImpl(RandomSeed);
375 : normalCached_(0.0),
376 normalCachedValid_(false)
378 this->seedImpl(theSeed);
386 template<
class Iterator>
388 : normalCached_(0.0),
389 normalCachedValid_(false)
391 this->seedImpl(init, length);
409 this->seedImpl(RandomSeed);
410 normalCachedValid_ =
false;
419 this->seedImpl(theSeed);
420 normalCachedValid_ =
false;
428 template<
class Iterator>
431 this->seedImpl(init, length);
432 normalCachedValid_ =
false;
454 #if 0 // difficult implementation necessary if low bits are not sufficiently random
461 UInt32 factor = factorForUniformInt(beyond);
462 UInt32 res = this->
get() / factor;
467 res = this->
get() / factor;
483 UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
484 UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
489 while(res > lastSafeValue)
502 return ( (this->
get() >> 5) * 67108864.0 + (this->
get() >> 6)) * (1.0/9007199254740992.0);
512 return (
double)this->
get() / 4294967295.0;
520 double uniform(
double lower,
double upper)
const
522 vigra_precondition(lower < upper,
523 "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
524 return uniform() * (upper-lower) + lower;
538 double normal(
double mean,
double stddev)
const
540 vigra_precondition(stddev > 0.0,
541 "RandomNumberGenerator::normal(): standard deviation must be positive.");
542 return normal()*stddev + mean;
558 return (range > 2147483648U || range == 0)
564 template <
class Engine>
567 if(normalCachedValid_)
569 normalCachedValid_ =
false;
570 return normalCached_;
577 x1 = uniform(-1.0, 1.0);
578 x2 = uniform(-1.0, 1.0);
579 w = x1 * x1 + x2 * x2;
581 while ( w > 1.0 || w == 0.0);
585 normalCached_ = x2 * w;
586 normalCachedValid_ =
true;
616 template <
class Engine>
617 class FunctorTraits<RandomNumberGenerator<Engine> >
620 typedef RandomNumberGenerator<Engine> type;
622 typedef VigraTrueType isInitializer;
624 typedef VigraFalseType isUnaryFunctor;
625 typedef VigraFalseType isBinaryFunctor;
626 typedef VigraFalseType isTernaryFunctor;
628 typedef VigraFalseType isUnaryAnalyser;
629 typedef VigraFalseType isBinaryAnalyser;
630 typedef VigraFalseType isTernaryAnalyser;
647 template <
class Engine = RandomTT800>
650 UInt32 lower_, difference_, factor_;
651 Engine
const & generator_;
665 : lower_(0), difference_(0xffffffff), factor_(1),
666 generator_(generator),
680 Engine
const & generator = Engine::global(),
681 bool useLowBits =
true)
682 : lower_(lower), difference_(upper-lower),
683 factor_(Engine::factorForUniformInt(difference_ + 1)),
684 generator_(generator),
685 useLowBits_(useLowBits)
687 vigra_precondition(lower < upper,
688 "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
695 if(difference_ == 0xffffffff)
698 return generator_.uniformInt(difference_+1) + lower_;
701 UInt32 res = generator_() / factor_;
705 while(res > difference_)
706 res = generator_() / factor_;
722 return generator_.uniformInt(beyond);
726 template <
class Engine>
727 class FunctorTraits<UniformIntRandomFunctor<Engine> >
730 typedef UniformIntRandomFunctor<Engine> type;
732 typedef VigraTrueType isInitializer;
734 typedef VigraTrueType isUnaryFunctor;
735 typedef VigraFalseType isBinaryFunctor;
736 typedef VigraFalseType isTernaryFunctor;
738 typedef VigraFalseType isUnaryAnalyser;
739 typedef VigraFalseType isBinaryAnalyser;
740 typedef VigraFalseType isTernaryAnalyser;
754 template <
class Engine = RandomTT800>
757 double offset_, scale_;
758 Engine
const & generator_;
772 generator_(generator)
781 Engine & generator = Engine::global())
783 scale_(upper - lower),
784 generator_(generator)
786 vigra_precondition(lower < upper,
787 "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
794 return generator_.uniform() * scale_ + offset_;
798 template <
class Engine>
799 class FunctorTraits<UniformRandomFunctor<Engine> >
802 typedef UniformRandomFunctor<Engine> type;
804 typedef VigraTrueType isInitializer;
806 typedef VigraFalseType isUnaryFunctor;
807 typedef VigraFalseType isBinaryFunctor;
808 typedef VigraFalseType isTernaryFunctor;
810 typedef VigraFalseType isUnaryAnalyser;
811 typedef VigraFalseType isBinaryAnalyser;
812 typedef VigraFalseType isTernaryAnalyser;
826 template <
class Engine = RandomTT800>
829 double mean_, stddev_;
830 Engine
const & generator_;
844 generator_(generator)
851 Engine & generator = Engine::global())
854 generator_(generator)
856 vigra_precondition(stddev > 0.0,
857 "NormalRandomFunctor(): standard deviation must be positive.");
864 return generator_.normal() * stddev_ + mean_;
869 template <
class Engine>
870 class FunctorTraits<NormalRandomFunctor<Engine> >
873 typedef UniformRandomFunctor<Engine> type;
875 typedef VigraTrueType isInitializer;
877 typedef VigraFalseType isUnaryFunctor;
878 typedef VigraFalseType isBinaryFunctor;
879 typedef VigraFalseType isTernaryFunctor;
881 typedef VigraFalseType isUnaryAnalyser;
882 typedef VigraFalseType isBinaryAnalyser;
883 typedef VigraFalseType isTernaryAnalyser;
890 #endif // VIGRA_RANDOM_HXX