|
1 /* -*- C++ -*- |
|
2 * |
|
3 * This file is a part of LEMON, a generic C++ optimization library |
|
4 * |
|
5 * Copyright (C) 2003-2007 |
|
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
7 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
8 * |
|
9 * Permission to use, modify and distribute this software is granted |
|
10 * provided that this copyright notice appears in all copies. For |
|
11 * precise terms see the accompanying LICENSE file. |
|
12 * |
|
13 * This software is provided "AS IS" with no warranty of any kind, |
|
14 * express or implied, and with no claim as to its suitability for any |
|
15 * purpose. |
|
16 * |
|
17 */ |
|
18 |
|
19 /* |
|
20 * This file contains the reimplemented version of the Mersenne Twister |
|
21 * Generator of Matsumoto and Nishimura. |
|
22 * |
|
23 * See the appropriate copyright notice below. |
|
24 * |
|
25 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, |
|
26 * All rights reserved. |
|
27 * |
|
28 * Redistribution and use in source and binary forms, with or without |
|
29 * modification, are permitted provided that the following conditions |
|
30 * are met: |
|
31 * |
|
32 * 1. Redistributions of source code must retain the above copyright |
|
33 * notice, this list of conditions and the following disclaimer. |
|
34 * |
|
35 * 2. Redistributions in binary form must reproduce the above copyright |
|
36 * notice, this list of conditions and the following disclaimer in the |
|
37 * documentation and/or other materials provided with the distribution. |
|
38 * |
|
39 * 3. The names of its contributors may not be used to endorse or promote |
|
40 * products derived from this software without specific prior written |
|
41 * permission. |
|
42 * |
|
43 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
|
44 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
|
45 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
|
46 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE |
|
47 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, |
|
48 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
|
49 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR |
|
50 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
|
51 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, |
|
52 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
|
53 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED |
|
54 * OF THE POSSIBILITY OF SUCH DAMAGE. |
|
55 * |
|
56 * |
|
57 * Any feedback is very welcome. |
|
58 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html |
|
59 * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) |
|
60 */ |
|
61 |
|
62 #ifndef LEMON_RANDOM_H |
|
63 #define LEMON_RANDOM_H |
|
64 |
|
65 #include <algorithm> |
|
66 #include <iterator> |
|
67 #include <vector> |
|
68 |
|
69 #include <ctime> |
|
70 #include <cmath> |
|
71 |
|
72 #include <lemon/dim2.h> |
|
73 ///\ingroup misc |
|
74 ///\file |
|
75 ///\brief Mersenne Twister random number generator |
|
76 /// |
|
77 ///\author Balazs Dezso |
|
78 |
|
79 namespace lemon { |
|
80 |
|
81 namespace _random_bits { |
|
82 |
|
83 template <typename _Word, int _bits = std::numeric_limits<_Word>::digits> |
|
84 struct RandomTraits {}; |
|
85 |
|
86 template <typename _Word> |
|
87 struct RandomTraits<_Word, 32> { |
|
88 |
|
89 typedef _Word Word; |
|
90 static const int bits = 32; |
|
91 |
|
92 static const int length = 624; |
|
93 static const int shift = 397; |
|
94 |
|
95 static const Word mul = 0x6c078965u; |
|
96 static const Word arrayInit = 0x012BD6AAu; |
|
97 static const Word arrayMul1 = 0x0019660Du; |
|
98 static const Word arrayMul2 = 0x5D588B65u; |
|
99 |
|
100 static const Word mask = 0x9908B0DFu; |
|
101 static const Word loMask = (1u << 31) - 1; |
|
102 static const Word hiMask = ~loMask; |
|
103 |
|
104 |
|
105 static Word tempering(Word rnd) { |
|
106 rnd ^= (rnd >> 11); |
|
107 rnd ^= (rnd << 7) & 0x9D2C5680u; |
|
108 rnd ^= (rnd << 15) & 0xEFC60000u; |
|
109 rnd ^= (rnd >> 18); |
|
110 return rnd; |
|
111 } |
|
112 |
|
113 }; |
|
114 |
|
115 template <typename _Word> |
|
116 struct RandomTraits<_Word, 64> { |
|
117 |
|
118 typedef _Word Word; |
|
119 static const int bits = 64; |
|
120 |
|
121 static const int length = 312; |
|
122 static const int shift = 156; |
|
123 |
|
124 static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du); |
|
125 static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu); |
|
126 static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u); |
|
127 static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu); |
|
128 |
|
129 static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u); |
|
130 static const Word loMask = (Word(1u) << 31) - 1; |
|
131 static const Word hiMask = ~loMask; |
|
132 |
|
133 static Word tempering(Word rnd) { |
|
134 rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u)); |
|
135 rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u)); |
|
136 rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u)); |
|
137 rnd ^= (rnd >> 43); |
|
138 return rnd; |
|
139 } |
|
140 |
|
141 }; |
|
142 |
|
143 template <typename _Word> |
|
144 class RandomCore { |
|
145 public: |
|
146 |
|
147 typedef _Word Word; |
|
148 |
|
149 private: |
|
150 |
|
151 static const int bits = RandomTraits<Word>::bits; |
|
152 |
|
153 static const int length = RandomTraits<Word>::length; |
|
154 static const int shift = RandomTraits<Word>::shift; |
|
155 |
|
156 public: |
|
157 |
|
158 void initState() { |
|
159 static const Word seedArray[4] = { |
|
160 0x12345u, 0x23456u, 0x34567u, 0x45678u |
|
161 }; |
|
162 |
|
163 initState(seedArray, seedArray + 4); |
|
164 } |
|
165 |
|
166 void initState(Word seed) { |
|
167 |
|
168 static const Word mul = RandomTraits<Word>::mul; |
|
169 |
|
170 current = state; |
|
171 |
|
172 Word *curr = state + length - 1; |
|
173 curr[0] = seed; --curr; |
|
174 for (int i = 1; i < length; ++i) { |
|
175 curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i); |
|
176 --curr; |
|
177 } |
|
178 } |
|
179 |
|
180 template <typename Iterator> |
|
181 void initState(Iterator begin, Iterator end) { |
|
182 |
|
183 static const Word init = RandomTraits<Word>::arrayInit; |
|
184 static const Word mul1 = RandomTraits<Word>::arrayMul1; |
|
185 static const Word mul2 = RandomTraits<Word>::arrayMul2; |
|
186 |
|
187 |
|
188 Word *curr = state + length - 1; --curr; |
|
189 Iterator it = begin; int cnt = 0; |
|
190 int num; |
|
191 |
|
192 initState(init); |
|
193 |
|
194 num = length > end - begin ? length : end - begin; |
|
195 while (num--) { |
|
196 curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) |
|
197 + *it + cnt; |
|
198 ++it; ++cnt; |
|
199 if (it == end) { |
|
200 it = begin; cnt = 0; |
|
201 } |
|
202 if (curr == state) { |
|
203 curr = state + length - 1; curr[0] = state[0]; |
|
204 } |
|
205 --curr; |
|
206 } |
|
207 |
|
208 num = length - 1; cnt = length - (curr - state) - 1; |
|
209 while (num--) { |
|
210 curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2)) |
|
211 - cnt; |
|
212 --curr; ++cnt; |
|
213 if (curr == state) { |
|
214 curr = state + length - 1; curr[0] = state[0]; --curr; |
|
215 cnt = 1; |
|
216 } |
|
217 } |
|
218 |
|
219 state[length - 1] = Word(1) << (bits - 1); |
|
220 } |
|
221 |
|
222 void copyState(const RandomCore& other) { |
|
223 std::copy(other.state, other.state + length, state); |
|
224 current = state + (other.current - other.state); |
|
225 } |
|
226 |
|
227 Word operator()() { |
|
228 if (current == state) fillState(); |
|
229 --current; |
|
230 Word rnd = *current; |
|
231 return RandomTraits<Word>::tempering(rnd); |
|
232 } |
|
233 |
|
234 private: |
|
235 |
|
236 |
|
237 void fillState() { |
|
238 static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask }; |
|
239 static const Word loMask = RandomTraits<Word>::loMask; |
|
240 static const Word hiMask = RandomTraits<Word>::hiMask; |
|
241 |
|
242 current = state + length; |
|
243 |
|
244 register Word *curr = state + length - 1; |
|
245 register long num; |
|
246 |
|
247 num = length - shift; |
|
248 while (num--) { |
|
249 curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ |
|
250 curr[- shift] ^ mask[curr[-1] & 1ul]; |
|
251 --curr; |
|
252 } |
|
253 num = shift - 1; |
|
254 while (num--) { |
|
255 curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ |
|
256 curr[length - shift] ^ mask[curr[-1] & 1ul]; |
|
257 --curr; |
|
258 } |
|
259 curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ |
|
260 curr[length - shift] ^ mask[curr[length - 1] & 1ul]; |
|
261 |
|
262 } |
|
263 |
|
264 |
|
265 Word *current; |
|
266 Word state[length]; |
|
267 |
|
268 }; |
|
269 |
|
270 |
|
271 template <typename Result, |
|
272 int shift = (std::numeric_limits<Result>::digits + 1) / 2> |
|
273 struct Masker { |
|
274 static Result mask(const Result& result) { |
|
275 return Masker<Result, (shift + 1) / 2>:: |
|
276 mask(static_cast<Result>(result | (result >> shift))); |
|
277 } |
|
278 }; |
|
279 |
|
280 template <typename Result> |
|
281 struct Masker<Result, 1> { |
|
282 static Result mask(const Result& result) { |
|
283 return static_cast<Result>(result | (result >> 1)); |
|
284 } |
|
285 }; |
|
286 |
|
287 template <typename Result, typename Word, |
|
288 int rest = std::numeric_limits<Result>::digits, int shift = 0, |
|
289 bool last = rest <= std::numeric_limits<Word>::digits> |
|
290 struct IntConversion { |
|
291 static const int bits = std::numeric_limits<Word>::digits; |
|
292 |
|
293 static Result convert(RandomCore<Word>& rnd) { |
|
294 return static_cast<Result>(rnd() >> (bits - rest)) << shift; |
|
295 } |
|
296 |
|
297 }; |
|
298 |
|
299 template <typename Result, typename Word, int rest, int shift> |
|
300 struct IntConversion<Result, Word, rest, shift, false> { |
|
301 static const int bits = std::numeric_limits<Word>::digits; |
|
302 |
|
303 static Result convert(RandomCore<Word>& rnd) { |
|
304 return (static_cast<Result>(rnd()) << shift) | |
|
305 IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd); |
|
306 } |
|
307 }; |
|
308 |
|
309 |
|
310 template <typename Result, typename Word, |
|
311 bool one_word = (std::numeric_limits<Word>::digits < |
|
312 std::numeric_limits<Result>::digits) > |
|
313 struct Mapping { |
|
314 static Result map(RandomCore<Word>& rnd, const Result& bound) { |
|
315 Word max = Word(bound - 1); |
|
316 Result mask = Masker<Result>::mask(bound - 1); |
|
317 Result num; |
|
318 do { |
|
319 num = IntConversion<Result, Word>::convert(rnd) & mask; |
|
320 } while (num > max); |
|
321 return num; |
|
322 } |
|
323 }; |
|
324 |
|
325 template <typename Result, typename Word> |
|
326 struct Mapping<Result, Word, false> { |
|
327 static Result map(RandomCore<Word>& rnd, const Result& bound) { |
|
328 Word max = Word(bound - 1); |
|
329 Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2> |
|
330 ::mask(max); |
|
331 Word num; |
|
332 do { |
|
333 num = rnd() & mask; |
|
334 } while (num > max); |
|
335 return num; |
|
336 } |
|
337 }; |
|
338 |
|
339 template <typename Result, int exp, bool pos = (exp >= 0)> |
|
340 struct ShiftMultiplier { |
|
341 static const Result multiplier() { |
|
342 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); |
|
343 res *= res; |
|
344 if ((exp & 1) == 1) res *= static_cast<Result>(2.0); |
|
345 return res; |
|
346 } |
|
347 }; |
|
348 |
|
349 template <typename Result, int exp> |
|
350 struct ShiftMultiplier<Result, exp, false> { |
|
351 static const Result multiplier() { |
|
352 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); |
|
353 res *= res; |
|
354 if ((exp & 1) == 1) res *= static_cast<Result>(0.5); |
|
355 return res; |
|
356 } |
|
357 }; |
|
358 |
|
359 template <typename Result> |
|
360 struct ShiftMultiplier<Result, 0, true> { |
|
361 static const Result multiplier() { |
|
362 return static_cast<Result>(1.0); |
|
363 } |
|
364 }; |
|
365 |
|
366 template <typename Result> |
|
367 struct ShiftMultiplier<Result, -20, true> { |
|
368 static const Result multiplier() { |
|
369 return static_cast<Result>(1.0/1048576.0); |
|
370 } |
|
371 }; |
|
372 |
|
373 template <typename Result> |
|
374 struct ShiftMultiplier<Result, -32, true> { |
|
375 static const Result multiplier() { |
|
376 return static_cast<Result>(1.0/424967296.0); |
|
377 } |
|
378 }; |
|
379 |
|
380 template <typename Result> |
|
381 struct ShiftMultiplier<Result, -53, true> { |
|
382 static const Result multiplier() { |
|
383 return static_cast<Result>(1.0/9007199254740992.0); |
|
384 } |
|
385 }; |
|
386 |
|
387 template <typename Result> |
|
388 struct ShiftMultiplier<Result, -64, true> { |
|
389 static const Result multiplier() { |
|
390 return static_cast<Result>(1.0/18446744073709551616.0); |
|
391 } |
|
392 }; |
|
393 |
|
394 template <typename Result, int exp> |
|
395 struct Shifting { |
|
396 static Result shift(const Result& result) { |
|
397 return result * ShiftMultiplier<Result, exp>::multiplier(); |
|
398 } |
|
399 }; |
|
400 |
|
401 template <typename Result, typename Word, |
|
402 int rest = std::numeric_limits<Result>::digits, int shift = 0, |
|
403 bool last = rest <= std::numeric_limits<Word>::digits> |
|
404 struct RealConversion{ |
|
405 static const int bits = std::numeric_limits<Word>::digits; |
|
406 |
|
407 static Result convert(RandomCore<Word>& rnd) { |
|
408 return Shifting<Result, - shift - rest>:: |
|
409 shift(static_cast<Result>(rnd() >> (bits - rest))); |
|
410 } |
|
411 }; |
|
412 |
|
413 template <typename Result, typename Word, int rest, int shift> |
|
414 struct RealConversion<Result, Word, rest, shift, false> { |
|
415 static const int bits = std::numeric_limits<Word>::digits; |
|
416 |
|
417 static Result convert(RandomCore<Word>& rnd) { |
|
418 return Shifting<Result, - shift - bits>:: |
|
419 shift(static_cast<Result>(rnd())) + |
|
420 RealConversion<Result, Word, rest-bits, shift + bits>:: |
|
421 convert(rnd); |
|
422 } |
|
423 }; |
|
424 |
|
425 template <typename Result, typename Word> |
|
426 struct Initializer { |
|
427 |
|
428 template <typename Iterator> |
|
429 static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) { |
|
430 std::vector<Word> ws; |
|
431 for (Iterator it = begin; it != end; ++it) { |
|
432 ws.push_back(Word(*it)); |
|
433 } |
|
434 rnd.initState(ws.begin(), ws.end()); |
|
435 } |
|
436 |
|
437 static void init(RandomCore<Word>& rnd, Result seed) { |
|
438 rnd.initState(seed); |
|
439 } |
|
440 }; |
|
441 |
|
442 template <typename Word> |
|
443 struct BoolConversion { |
|
444 static bool convert(RandomCore<Word>& rnd) { |
|
445 return (rnd() & 1) == 1; |
|
446 } |
|
447 }; |
|
448 |
|
449 template <typename Word> |
|
450 struct BoolProducer { |
|
451 Word buffer; |
|
452 int num; |
|
453 |
|
454 BoolProducer() : num(0) {} |
|
455 |
|
456 bool convert(RandomCore<Word>& rnd) { |
|
457 if (num == 0) { |
|
458 buffer = rnd(); |
|
459 num = RandomTraits<Word>::bits; |
|
460 } |
|
461 bool r = (buffer & 1); |
|
462 buffer >>= 1; |
|
463 --num; |
|
464 return r; |
|
465 } |
|
466 }; |
|
467 |
|
468 } |
|
469 |
|
470 /// \ingroup misc |
|
471 /// |
|
472 /// \brief Mersenne Twister random number generator |
|
473 /// |
|
474 /// The Mersenne Twister is a twisted generalized feedback |
|
475 /// shift-register generator of Matsumoto and Nishimura. The period |
|
476 /// of this generator is \f$ 2^{19937} - 1 \f$ and it is |
|
477 /// equi-distributed in 623 dimensions for 32-bit numbers. The time |
|
478 /// performance of this generator is comparable to the commonly used |
|
479 /// generators. |
|
480 /// |
|
481 /// This implementation is specialized for both 32-bit and 64-bit |
|
482 /// architectures. The generators differ sligthly in the |
|
483 /// initialization and generation phase so they produce two |
|
484 /// completly different sequences. |
|
485 /// |
|
486 /// The generator gives back random numbers of serveral types. To |
|
487 /// get a random number from a range of a floating point type you |
|
488 /// can use one form of the \c operator() or the \c real() member |
|
489 /// function. If you want to get random number from the {0, 1, ..., |
|
490 /// n-1} integer range use the \c operator[] or the \c integer() |
|
491 /// method. And to get random number from the whole range of an |
|
492 /// integer type you can use the argumentless \c integer() or \c |
|
493 /// uinteger() functions. After all you can get random bool with |
|
494 /// equal chance of true and false or given probability of true |
|
495 /// result with the \c boolean() member functions. |
|
496 /// |
|
497 ///\code |
|
498 /// // The commented code is identical to the other |
|
499 /// double a = rnd(); // [0.0, 1.0) |
|
500 /// // double a = rnd.real(); // [0.0, 1.0) |
|
501 /// double b = rnd(100.0); // [0.0, 100.0) |
|
502 /// // double b = rnd.real(100.0); // [0.0, 100.0) |
|
503 /// double c = rnd(1.0, 2.0); // [1.0, 2.0) |
|
504 /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) |
|
505 /// int d = rnd[100000]; // 0..99999 |
|
506 /// // int d = rnd.integer(100000); // 0..99999 |
|
507 /// int e = rnd[6] + 1; // 1..6 |
|
508 /// // int e = rnd.integer(1, 1 + 6); // 1..6 |
|
509 /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1 |
|
510 /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1 |
|
511 /// bool g = rnd.boolean(); // P(g = true) = 0.5 |
|
512 /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 |
|
513 ///\endcode |
|
514 /// |
|
515 /// The lemon provides a global instance of the random number |
|
516 /// generator which name is \ref lemon::rnd "rnd". Usually it is a |
|
517 /// good programming convenience to use this global generator to get |
|
518 /// random numbers. |
|
519 /// |
|
520 /// \author Balazs Dezso |
|
521 class Random { |
|
522 private: |
|
523 |
|
524 // architecture word |
|
525 typedef unsigned long Word; |
|
526 |
|
527 _random_bits::RandomCore<Word> core; |
|
528 _random_bits::BoolProducer<Word> bool_producer; |
|
529 |
|
530 |
|
531 public: |
|
532 |
|
533 /// \brief Constructor |
|
534 /// |
|
535 /// Constructor with constant seeding. |
|
536 Random() { core.initState(); } |
|
537 |
|
538 /// \brief Constructor |
|
539 /// |
|
540 /// Constructor with seed. The current number type will be converted |
|
541 /// to the architecture word type. |
|
542 template <typename Number> |
|
543 Random(Number seed) { |
|
544 _random_bits::Initializer<Number, Word>::init(core, seed); |
|
545 } |
|
546 |
|
547 /// \brief Constructor |
|
548 /// |
|
549 /// Constructor with array seeding. The given range should contain |
|
550 /// any number type and the numbers will be converted to the |
|
551 /// architecture word type. |
|
552 template <typename Iterator> |
|
553 Random(Iterator begin, Iterator end) { |
|
554 typedef typename std::iterator_traits<Iterator>::value_type Number; |
|
555 _random_bits::Initializer<Number, Word>::init(core, begin, end); |
|
556 } |
|
557 |
|
558 /// \brief Copy constructor |
|
559 /// |
|
560 /// Copy constructor. The generated sequence will be identical to |
|
561 /// the other sequence. It can be used to save the current state |
|
562 /// of the generator and later use it to generate the same |
|
563 /// sequence. |
|
564 Random(const Random& other) { |
|
565 core.copyState(other.core); |
|
566 } |
|
567 |
|
568 /// \brief Assign operator |
|
569 /// |
|
570 /// Assign operator. The generated sequence will be identical to |
|
571 /// the other sequence. It can be used to save the current state |
|
572 /// of the generator and later use it to generate the same |
|
573 /// sequence. |
|
574 Random& operator=(const Random& other) { |
|
575 if (&other != this) { |
|
576 core.copyState(other.core); |
|
577 } |
|
578 return *this; |
|
579 } |
|
580 |
|
581 /// \brief Returns a random real number from the range [0, 1) |
|
582 /// |
|
583 /// It returns a random real number from the range [0, 1). The |
|
584 /// default Number type is double. |
|
585 template <typename Number> |
|
586 Number real() { |
|
587 return _random_bits::RealConversion<Number, Word>::convert(core); |
|
588 } |
|
589 |
|
590 double real() { |
|
591 return real<double>(); |
|
592 } |
|
593 |
|
594 /// \brief Returns a random real number the range [0, b) |
|
595 /// |
|
596 /// It returns a random real number from the range [0, b). |
|
597 template <typename Number> |
|
598 Number real(Number b) { |
|
599 return real<Number>() * b; |
|
600 } |
|
601 |
|
602 /// \brief Returns a random real number from the range [a, b) |
|
603 /// |
|
604 /// It returns a random real number from the range [a, b). |
|
605 template <typename Number> |
|
606 Number real(Number a, Number b) { |
|
607 return real<Number>() * (b - a) + a; |
|
608 } |
|
609 |
|
610 /// \brief Returns a random real number from the range [0, 1) |
|
611 /// |
|
612 /// It returns a random double from the range [0, 1). |
|
613 double operator()() { |
|
614 return real<double>(); |
|
615 } |
|
616 |
|
617 /// \brief Returns a random real number from the range [0, b) |
|
618 /// |
|
619 /// It returns a random real number from the range [0, b). |
|
620 template <typename Number> |
|
621 Number operator()(Number b) { |
|
622 return real<Number>() * b; |
|
623 } |
|
624 |
|
625 /// \brief Returns a random real number from the range [a, b) |
|
626 /// |
|
627 /// It returns a random real number from the range [a, b). |
|
628 template <typename Number> |
|
629 Number operator()(Number a, Number b) { |
|
630 return real<Number>() * (b - a) + a; |
|
631 } |
|
632 |
|
633 /// \brief Returns a random integer from a range |
|
634 /// |
|
635 /// It returns a random integer from the range {0, 1, ..., b - 1}. |
|
636 template <typename Number> |
|
637 Number integer(Number b) { |
|
638 return _random_bits::Mapping<Number, Word>::map(core, b); |
|
639 } |
|
640 |
|
641 /// \brief Returns a random integer from a range |
|
642 /// |
|
643 /// It returns a random integer from the range {a, a + 1, ..., b - 1}. |
|
644 template <typename Number> |
|
645 Number integer(Number a, Number b) { |
|
646 return _random_bits::Mapping<Number, Word>::map(core, b - a) + a; |
|
647 } |
|
648 |
|
649 /// \brief Returns a random integer from a range |
|
650 /// |
|
651 /// It returns a random integer from the range {0, 1, ..., b - 1}. |
|
652 template <typename Number> |
|
653 Number operator[](Number b) { |
|
654 return _random_bits::Mapping<Number, Word>::map(core, b); |
|
655 } |
|
656 |
|
657 /// \brief Returns a random non-negative integer |
|
658 /// |
|
659 /// It returns a random non-negative integer uniformly from the |
|
660 /// whole range of the current \c Number type. The default result |
|
661 /// type of this function is unsigned int. |
|
662 template <typename Number> |
|
663 Number uinteger() { |
|
664 return _random_bits::IntConversion<Number, Word>::convert(core); |
|
665 } |
|
666 |
|
667 unsigned int uinteger() { |
|
668 return uinteger<unsigned int>(); |
|
669 } |
|
670 |
|
671 /// \brief Returns a random integer |
|
672 /// |
|
673 /// It returns a random integer uniformly from the whole range of |
|
674 /// the current \c Number type. The default result type of this |
|
675 /// function is int. |
|
676 template <typename Number> |
|
677 Number integer() { |
|
678 static const int nb = std::numeric_limits<Number>::digits + |
|
679 (std::numeric_limits<Number>::is_signed ? 1 : 0); |
|
680 return _random_bits::IntConversion<Number, Word, nb>::convert(core); |
|
681 } |
|
682 |
|
683 int integer() { |
|
684 return integer<int>(); |
|
685 } |
|
686 |
|
687 /// \brief Returns a random bool |
|
688 /// |
|
689 /// It returns a random bool. The generator holds a buffer for |
|
690 /// random bits. Every time when it become empty the generator makes |
|
691 /// a new random word and fill the buffer up. |
|
692 bool boolean() { |
|
693 return bool_producer.convert(core); |
|
694 } |
|
695 |
|
696 ///\name Nonuniform distributions |
|
697 /// |
|
698 |
|
699 ///@{ |
|
700 |
|
701 /// \brief Returns a random bool |
|
702 /// |
|
703 /// It returns a random bool with given probability of true result |
|
704 bool boolean(double p) { |
|
705 return operator()() < p; |
|
706 } |
|
707 |
|
708 /// Standard Gauss distribution |
|
709 |
|
710 /// Standard Gauss distribution. |
|
711 /// \note The Cartesian form of the Box-Muller |
|
712 /// transformation is used to generate a random normal distribution. |
|
713 /// \todo Consider using the "ziggurat" method instead. |
|
714 double gauss() |
|
715 { |
|
716 double V1,V2,S; |
|
717 do { |
|
718 V1=2*real<double>()-1; |
|
719 V2=2*real<double>()-1; |
|
720 S=V1*V1+V2*V2; |
|
721 } while(S>=1); |
|
722 return std::sqrt(-2*std::log(S)/S)*V1; |
|
723 } |
|
724 /// Gauss distribution with given mean and standard deviation |
|
725 |
|
726 /// \sa gauss() |
|
727 /// |
|
728 double gauss(double mean,double std_dev) |
|
729 { |
|
730 return gauss()*std_dev+mean; |
|
731 } |
|
732 |
|
733 /// Exponential distribution with given mean |
|
734 |
|
735 /// This function generates an exponential distribution random number |
|
736 /// with mean <tt>1/lambda</tt>. |
|
737 /// |
|
738 double exponential(double lambda=1.0) |
|
739 { |
|
740 return -std::log(real<double>())/lambda; |
|
741 } |
|
742 |
|
743 /// Gamma distribution with given integer shape |
|
744 |
|
745 /// This function generates a gamma distribution random number. |
|
746 /// |
|
747 ///\param k shape parameter (<tt>k>0</tt> integer) |
|
748 double gamma(int k) |
|
749 { |
|
750 double s = 0; |
|
751 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); |
|
752 return s; |
|
753 } |
|
754 |
|
755 /// Gamma distribution with given shape and scale parameter |
|
756 |
|
757 /// This function generates a gamma distribution random number. |
|
758 /// |
|
759 ///\param k shape parameter (<tt>k>0</tt>) |
|
760 ///\param theta scale parameter |
|
761 /// |
|
762 double gamma(double k,double theta=1.0) |
|
763 { |
|
764 double xi,nu; |
|
765 const double delta = k-std::floor(k); |
|
766 const double v0=M_E/(M_E-delta); |
|
767 do { |
|
768 double V0=1.0-real<double>(); |
|
769 double V1=1.0-real<double>(); |
|
770 double V2=1.0-real<double>(); |
|
771 if(V2<=v0) |
|
772 { |
|
773 xi=std::pow(V1,1.0/delta); |
|
774 nu=V0*std::pow(xi,delta-1.0); |
|
775 } |
|
776 else |
|
777 { |
|
778 xi=1.0-std::log(V1); |
|
779 nu=V0*std::exp(-xi); |
|
780 } |
|
781 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); |
|
782 return theta*(xi-gamma(int(std::floor(k)))); |
|
783 } |
|
784 |
|
785 |
|
786 ///@} |
|
787 |
|
788 ///\name Two dimensional distributions |
|
789 /// |
|
790 |
|
791 ///@{ |
|
792 |
|
793 /// Uniform distribution on the full unit circle. |
|
794 dim2::Point<double> disc() |
|
795 { |
|
796 double V1,V2; |
|
797 do { |
|
798 V1=2*real<double>()-1; |
|
799 V2=2*real<double>()-1; |
|
800 |
|
801 } while(V1*V1+V2*V2>=1); |
|
802 return dim2::Point<double>(V1,V2); |
|
803 } |
|
804 /// A kind of two dimensional Gauss distribution |
|
805 |
|
806 /// This function provides a turning symmetric two-dimensional distribution. |
|
807 /// Both coordinates are of standard normal distribution, but they are not |
|
808 /// independent. |
|
809 /// |
|
810 /// \note The coordinates are the two random variables provided by |
|
811 /// the Box-Muller method. |
|
812 dim2::Point<double> gauss2() |
|
813 { |
|
814 double V1,V2,S; |
|
815 do { |
|
816 V1=2*real<double>()-1; |
|
817 V2=2*real<double>()-1; |
|
818 S=V1*V1+V2*V2; |
|
819 } while(S>=1); |
|
820 double W=std::sqrt(-2*std::log(S)/S); |
|
821 return dim2::Point<double>(W*V1,W*V2); |
|
822 } |
|
823 /// A kind of two dimensional exponential distribution |
|
824 |
|
825 /// This function provides a turning symmetric two-dimensional distribution. |
|
826 /// The x-coordinate is of conditionally exponential distribution |
|
827 /// with the condition that x is positive and y=0. If x is negative and |
|
828 /// y=0 then, -x is of exponential distribution. The same is true for the |
|
829 /// y-coordinate. |
|
830 dim2::Point<double> exponential2() |
|
831 { |
|
832 double V1,V2,S; |
|
833 do { |
|
834 V1=2*real<double>()-1; |
|
835 V2=2*real<double>()-1; |
|
836 S=V1*V1+V2*V2; |
|
837 } while(S>=1); |
|
838 double W=-std::log(S)/S; |
|
839 return dim2::Point<double>(W*V1,W*V2); |
|
840 } |
|
841 |
|
842 ///@} |
|
843 }; |
|
844 |
|
845 |
|
846 extern Random rnd; |
|
847 |
|
848 } |
|
849 |
|
850 #endif |