Functions/rng/PkRandom.h

   1
   2 //Mersenne Twister pseudorandom number generator of 32 and 64-bit.
   3 //2005, Diego Park <diegopark@gmail.com>
   4 //
   5 //A C-program for MT19937, with initialization improved 2002/1/26.
   6 //Coded by Takuji Nishimura and Makoto Matsumoto.
   7 //
   8 //Before using, initialize the state by using init_genrand(seed)
   9 //or init_by_array(init_key, key_length).
  10 //
  11 //Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  12 //All rights reserved.
  13 //
  14 //A C-program for MT19937-64 (2004/9/29 version).
  15 //Coded by Takuji Nishimura and Makoto Matsumoto.
  16 //
  17 //This is a 64-bit version of Mersenne Twister pseudorandom number
  18 //generator.
  19 //
  20 //Before using, initialize the state by using init_genrand64(seed)
  21 //or init_by_array64(init_key, key_length).
  22 //
  23 //Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
  24 //All rights reserved.
  25 //
  26 //Redistribution and use in source and binary forms, with or without
  27 //modification, are permitted provided that the following conditions
  28 //are met:
  29 //
  30 //1. Redistributions of source code must retain the above copyright
  31 //notice, this list of conditions and the following disclaimer.
  32 //
  33 //2. Redistributions in binary form must reproduce the above copyright
  34 //notice, this list of conditions and the following disclaimer in the
  35 //documentation and/or other materials provided with the distribution.
  36 //
  37 //3. The names of its contributors may not be used to endorse or promote
  38 //products derived from this software without specific prior written
  39 //permission.
  40 //
  41 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  42 //"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  43 //LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  44 //A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
  45 //CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  46 //EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  47 //PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  48 //PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  49 //LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  50 //NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  51 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  52 //
  53 //References:
  54 //T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
  55 //ACM Transactions on Modeling and
  56 //Computer Simulation 10. (2000) 348--357.
  57 //M. Matsumoto and T. Nishimura,
  58 //``Mersenne Twister: a 623-dimensionally equidistributed
  59 //uniform pseudorandom number generator''
  60 //ACM Transactions on Modeling and
  61 //Computer Simulation 8. (Jan. 1998) 3--30.
  62 //
  63 //Any feedback is very welcome.
  64 //http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html
  65 //email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)
  66
  67 #ifndef PKRANDOM_H
  68 #define PKRANDOM_H
  69
  70 #include <cstring>
  71
  72 #include <boost/random/mersenne_twister.hpp>
  73
  74 typedef boost::mt19937 rng_uint32_type;
  75
  76 // Mersenne Twister
  77
  78 namespace Pk {
  79
  80 template <class Integer>
  81 class Random {
  82 public:
  83         Random(Integer aSeed = 5489);
  84
  85         Integer next();
  86
  87 private:
  88         static const size_t MaxState = 624 * sizeof(unsigned long) / sizeof(Integer);
  89         static const size_t Period;
  90
  91         void seed();
  92
  93         static Integer hashFunction(Integer anInteger);
  94         static Integer twist(Integer a, Integer b, Integer c);
  95
  96         Integer state[MaxState];
  97         Integer index;
  98 };
  99
 100 typedef Random<unsigned long> Random32;
 101 typedef Random<unsigned long long> Random64;
 102
 103 #include "PkRandom.inl"
 104
 105 // Added and implemented by MaRu.
 106 // This is the algorithm used in Java Standard Library, but I'm not convinced if it's the correct way to do it.
 107 // IMHO The correct way should be based on IEEE 754 floating point number encoding standard.
 108 inline double nextDouble(rng_uint32_type &rng) {
 109     unsigned long long l = ((((unsigned long long)rng()) & 0x03FFFFFFLL) << 27) | (((unsigned long long)rng()) & 0x7FFFFFF);
 110     return l / (double)(1LL << 53);
 111 }
 112
 113 inline double nextDouble(Random64& rng) {
 114     unsigned long long l = (rng.next() & 0x1FFFFFFFFFFFFFLL);
 115     return l / (double)(1LL << 53);
 116 }
 117
 118 }
 119
 120 #endif