COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/random.h @ 2229:4dbb6dd2dd4b

Last change on this file since 2229:4dbb6dd2dd4b was 2229:4dbb6dd2dd4b, checked in by Balazs Dezso, 13 years ago

Mersenne Twister random number generator

The code is based on the official MT19937 implementation
It is fully rewritten:

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

todo: fixing copyright information

File size: 14.8 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2006
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#ifndef LEMON_RANDOM_H
20#define LEMON_RANDOM_H
21
22#include <algorithm>
23
24#include <ctime>
25#include <cmath>
26#include <cmath>
27
28///\ingroup misc
29///\file
30///\brief Mersenne Twister random number generator
31///
32///\author Balazs Dezso
33
34namespace lemon {
35
36#if WORD_BIT == 32
37
38  /// \ingroup misc
39  ///
40  /// \brief Mersenne Twister random number generator
41  ///
42  /// The Mersenne Twister is a twisted generalized feedback
43  /// shift-register generator of Matsumoto and Nishimura. The period of
44  /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
45  /// 623 dimensions. The time performance of this generator is comparable
46  /// to the commonly used generators.
47  ///
48  /// \author Balazs Dezso
49  class Random {
50
51    static const int length = 624;
52    static const int shift = 397;
53
54  public:
55
56    static const unsigned long min = 0ul;
57    static const unsigned long max = ~0ul;
58 
59    /// \brief Constructor
60    ///
61    /// Constructor with time dependent seeding.
62    Random() { initState(std::time(0)); }
63
64    /// \brief Constructor
65    ///
66    /// Constructor
67    Random(unsigned long seed) { initState(seed); }
68
69    /// \brief Copy constructor
70    ///
71    /// Copy constructor. The generated sequence will be identical to
72    /// the other sequence.
73    Random(const Random& other) {
74      std::copy(other.state, other.state + length, state);
75      current = state + (other.current - other.state);
76    }
77
78    /// \brief Assign operator
79    ///
80    /// Assign operator. The generated sequence will be identical to
81    /// the other sequence.
82    Random& operator=(const Random& other) {
83      if (&other != this) {
84        std::copy(other.state, other.state + length, state);
85        current = state + (other.current - other.state);
86      }
87      return *this;
88    }
89
90    /// \brief Returns an unsigned random number
91    ///
92    /// It returns an unsigned integer random number from the range
93    /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$.
94    unsigned long getUnsigned() {
95      if (current == state) fillState();
96      --current;
97      unsigned long rnd = *current;
98      rnd ^= (rnd >> 11);
99      rnd ^= (rnd << 7) & 0x9D2C5680ul;
100      rnd ^= (rnd << 15) & 0xEFC60000ul;
101      rnd ^= (rnd >> 18);
102      return rnd;
103    }
104
105    /// \brief Returns a random number
106    ///
107    /// It returns an integer random number from the range
108    /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$.
109    long getInt() {
110      return (long)getUnsigned();
111    }
112   
113    /// \brief Returns an unsigned integer random number
114    ///
115    /// It returns an unsigned integer random number from the range
116    /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$.
117    long getNatural() {
118      return (long)(getUnsigned() >> 1);
119    }
120
121    /// \brief Returns a random bool
122    ///
123    /// It returns a random bool.
124    bool getBool() {
125      return (bool)(getUnsigned() & 1);
126    }
127
128    /// \brief Returns a real random number
129    ///
130    /// It returns a real random number from the range
131    /// \f$ [0, 1) \f$. The double will have 32 significant bits.
132    double getReal() {
133      return std::ldexp((double)getUnsigned(), -32);
134    }
135
136    /// \brief Returns a real random number
137    ///
138    /// It returns a real random number from the range
139    /// \f$ [0, 1) \f$. The double will have 53 significant bits,
140    /// but it is slower than the \c getReal().
141    double getPrecReal() {
142      return std::ldexp((double)(getUnsigned() >> 5), -27) +
143        std::ldexp((double)(getUnsigned() >> 6), -53);
144    }
145
146    /// \brief Returns an unsigned random number from a given range
147    ///
148    /// It returns an unsigned integer random number from the range
149    /// \f$ \{0, 1, \dots, n - 1\} \f$.
150    unsigned long getUnsigned(unsigned long n) {
151      unsigned long mask = n - 1, rnd;
152      mask |= mask >> 1;
153      mask |= mask >> 2;
154      mask |= mask >> 4;
155      mask |= mask >> 8;
156      mask |= mask >> 16;
157      do rnd = getUnsigned() & mask; while (rnd >= n);
158      return rnd;
159    }
160
161    /// \brief Returns a random number from a given range
162    ///
163    /// It returns an unsigned integer random number from the range
164    /// \f$ \{0, 1, \dots, n - 1\} \f$.
165    long getInt(long n) {
166      long mask = n - 1, rnd;
167      mask |= mask >> 1;
168      mask |= mask >> 2;
169      mask |= mask >> 4;
170      mask |= mask >> 8;
171      mask |= mask >> 16;
172      do rnd = getUnsigned() & mask; while (rnd >= n);
173      return rnd;
174    }
175
176  private:
177
178    void initState(unsigned long seed) {
179      static const unsigned long mul = 0x6c078965ul;
180
181      current = state;
182
183      unsigned long *curr = state + length - 1;
184      curr[0] = seed; --curr;
185      for (int i = 1; i < length; ++i) {
186        curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i);
187        --curr;
188      }
189    }
190 
191    void fillState() {
192      static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul };
193      static const unsigned long loMask = (1ul << 31) - 1;
194      static const unsigned long hiMask = ~loMask;
195
196      current = state + length;
197
198      register unsigned long *curr = state + length - 1;
199      register long num;
200     
201      num = length - shift;
202      while (num--) {
203        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
204          curr[- shift] ^ mask[curr[-1] & 1ul];
205        --curr;
206      }
207      num = shift - 1;
208      while (num--) {
209        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
210          curr[length - shift] ^ mask[curr[-1] & 1ul];
211        --curr;
212      }
213      curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
214        curr[length - shift] ^ mask[curr[length - 1] & 1ul];
215
216    }
217 
218    unsigned long *current;
219    unsigned long state[length];
220   
221  };
222
223#elif WORD_BIT == 64
224
225  /// \brief Mersenne Twister random number generator
226  ///
227  /// The Mersenne Twister is a twisted generalized feedback
228  /// shift-register generator of Matsumoto and Nishimura. The period of
229  /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
230  /// 311 dimensions. The time performance of this generator is comparable
231  /// to the commonly used generators.
232  class Random {
233
234    static const int length = 312;
235    static const int shift = 156;
236
237  public:
238
239    static const unsigned long min = 0ul;
240    static const unsigned long max = ~0ul;
241 
242    /// \brief Constructor
243    ///
244    /// Constructor with time dependent seeding.
245    Random() { initState(std::time(0)); }
246
247    /// \brief Constructor
248    ///
249    /// Constructor
250    Random(unsigned long seed) { initState(seed); }
251
252    /// \brief Copy constructor
253    ///
254    /// Copy constructor. The generated sequence will be identical to
255    /// the other sequence.
256    Random(const Random& other) {
257      std::copy(other.state, other.state + length, state);
258      current = state + (other.current - other.state);
259    }
260
261    /// \brief Assign operator
262    ///
263    /// Assign operator. The generated sequence will be identical to
264    /// the other sequence.
265    Random& operator=(const Random& other) {
266      if (&other != this) {
267        std::copy(other.state, other.state + length, state);
268        current = state + (other.current - other.state);
269      }
270      return *this;
271    }
272
273    /// \brief Returns an unsigned random number
274    ///
275    /// It returns an unsigned integer random number from the range
276    /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$.
277    unsigned long getUnsigned() {
278      if (current == state) fillState();
279      --current;
280      unsigned long rnd = *current;
281      rnd ^= (rnd >> 29) & 0x5555555555555555ul;
282      rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;
283      rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;
284      rnd ^= (rnd >> 43);
285      return rnd;
286    }
287
288    /// \brief Returns a random number
289    ///
290    /// It returns an integer random number from the range
291    /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$.
292    long getInt() {
293      return (long)getUnsigned();
294    }
295   
296    /// \brief Returns an unsigned integer random number
297    ///
298    /// It returns an unsigned integer random number from the range
299    /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$.
300    long getNatural() {
301      return (long)(getUnsigned() >> 1);
302    }
303
304    /// \brief Returns a random bool
305    ///
306    /// It returns a random bool.
307    bool getBool() {
308      return (bool)(getUnsigned() & 1);
309    }
310
311    /// \brief Returns a real random number
312    ///
313    /// It returns a real random number from the range
314    /// \f$ [0, 1) \f$.
315    double getReal() {
316      return std::ldexp((double)(getUnsigned() >> 11), -53);
317    }
318
319    /// \brief Returns a real random number
320    ///
321    /// It returns a real random number from the range
322    /// \f$ [0, 1) \f$. This function is identical to the \c getReal().
323    double getPrecReal() {
324      return getReal();
325    }
326
327    /// \brief Returns an unsigned random number from a given range
328    ///
329    /// It returns an unsigned integer random number from the range
330    /// \f$ \{0, 1, \dots, n - 1\} \f$.
331    unsigned long getUnsigned(unsigned long n) {
332      unsigned long mask = n - 1, rnd;
333      mask |= mask >> 1;
334      mask |= mask >> 2;
335      mask |= mask >> 4;
336      mask |= mask >> 8;
337      mask |= mask >> 16;
338      mask |= mask >> 32;
339      do rnd = getUnsigned() & mask; while (rnd >= n);
340      return rnd;
341    }
342
343    /// \brief Returns a random number from a given range
344    ///
345    /// It returns an unsigned integer random number from the range
346    /// \f$ \{0, 1, \dots, n - 1\} \f$.
347    long getInt(long n) {
348      long mask = n - 1, rnd;
349      mask |= mask >> 1;
350      mask |= mask >> 2;
351      mask |= mask >> 4;
352      mask |= mask >> 8;
353      mask |= mask >> 16;
354      mask |= mask >> 32;
355      do rnd = getUnsigned() & mask; while (rnd >= n);
356      return rnd;
357    }
358
359  private:
360
361    void initState(unsigned long seed) {
362
363      static const unsigned long mul = 0x5851F42D4C957F2Dul;
364
365      current = state;
366
367      unsigned long *curr = state + length - 1;
368      curr[0] = seed; --curr;
369      for (int i = 1; i < length; ++i) {
370        curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);
371        --curr;
372      }
373    }
374 
375    void fillState() {
376      static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };
377      static const unsigned long loMask = (1ul << 31) - 1;
378      static const unsigned long hiMask = ~loMask;
379
380      current = state + length;
381
382      register unsigned long *curr = state + length - 1;
383      register int num;
384     
385      num = length - shift;
386      while (num--) {
387        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
388          curr[- shift] ^ mask[curr[-1] & 1ul];
389        --curr;
390      }
391      num = shift - 1;
392      while (num--) {
393        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
394          curr[length - shift] ^ mask[curr[-1] & 1ul];
395        --curr;
396      }
397      curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
398        curr[length - shift] ^ mask[curr[length - 1] & 1ul];
399
400    }
401 
402    unsigned long *current;
403    unsigned long state[length];
404   
405  };
406
407#else
408
409  /// \brief Mersenne Twister random number generator
410  ///
411  /// The Mersenne Twister is a twisted generalized feedback
412  /// shift-register generator of Matsumoto and Nishimura. There is
413  /// two different implementation in the Lemon library, one for
414  /// 32-bit processors and one for 64-bit processors. The period of
415  /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated
416  /// sequence of 32-bit random numbers is equi-distributed in 623
417  /// dimensions. The time performance of this generator is comparable
418  /// to the commonly used generators.
419  class Random {
420  public:
421
422    static const unsigned long min = 0ul;
423    static const unsigned long max = ~0ul;
424 
425    /// \brief Constructor
426    ///
427    /// Constructor with time dependent seeding.
428    Random() {}
429
430    /// \brief Constructor
431    ///
432    /// Constructor
433    Random(unsigned long seed) {}
434
435    /// \brief Copy constructor
436    ///
437    /// Copy constructor. The generated sequence will be identical to
438    /// the other sequence.
439    Random(const Random& other) {}
440
441    /// \brief Assign operator
442    ///
443    /// Assign operator. The generated sequence will be identical to
444    /// the other sequence.
445    Random& operator=(const Random& other) { return *this; }
446
447    /// \brief Returns an unsigned random number
448    ///
449    /// It returns an unsigned integer random number from the range
450    /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from
451    /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors.
452    unsigned long getUnsigned() { return 0ul; }
453
454    /// \brief Returns a random number
455    ///
456    /// It returns an integer random number from the range
457    /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from
458    /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
459    long getInt() { return 0l; }
460   
461    /// \brief Returns an unsigned integer random number
462    ///
463    /// It returns an unsigned integer random number from the range
464    /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and
465    /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
466    long getNatural() { return 0l; }
467
468    /// \brief Returns a random bool
469    ///
470    /// It returns a random bool.
471    bool getBool() { return false; }
472
473    /// \brief Returns a real random number
474    ///
475    /// It returns a real random number from the range
476    /// \f$ [0, 1) \f$. For 32-bit processors the generated random
477    /// number will just have 32 significant bits.
478    double getReal() { return 0.0; }
479
480    /// \brief Returns a real random number
481    ///
482    /// It returns a real random number from the range
483    /// \f$ [0, 1) \f$. This function returns random numbers with 53
484    /// significant bits for 32-bit processors. For 64-bit processors
485    /// it is identical to the \c getReal().
486    double getPrecReal() { return 0.0; }
487
488    /// \brief Returns an unsigned random number from a given range
489    ///
490    /// It returns an unsigned integer random number from the range
491    /// \f$ \{0, 1, \dots, n - 1\} \f$.
492    unsigned long getUnsigned(unsigned long n) { return 0; }
493
494    /// \brief Returns a random number from a given range
495    ///
496    /// It returns an unsigned integer random number from the range
497    /// \f$ \{0, 1, \dots, n - 1\} \f$.
498    long getInt(long n) { return 0; }
499
500  };
501
502#endif
503
504  /// \brief Global random number generator instance
505  ///
506  /// A global mersenne twister random number generator instance
507  extern Random rnd;
508
509}
510
511#endif
Note: See TracBrowser for help on using the repository browser.