libstdc++-v3/include/parallel/random_number.h

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2007, 2008 Free Software Foundation, Inc.
   4 //
   5 // This file is part of the GNU ISO C++ Library.  This library is free
   6 // software; you can redistribute it and/or modify it under the terms
   7 // of the GNU General Public License as published by the Free Software
   8 // Foundation; either version 2, or (at your option) any later
   9 // version.
  10
  11 // This library is distributed in the hope that it will be useful, but
  12 // WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // General Public License for more details.
  15
  16 // You should have received a copy of the GNU General Public License
  17 // along with this library; see the file COPYING.  If not, write to
  18 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
  19 // MA 02111-1307, USA.
  20
  21 // As a special exception, you may use this file as part of a free
  22 // software library without restriction.  Specifically, if other files
  23 // instantiate templates or use macros or inline functions from this
  24 // file, or you compile this file and link it with other files to
  25 // produce an executable, this file does not by itself cause the
  26 // resulting executable to be covered by the GNU General Public
  27 // License.  This exception does not however invalidate any other
  28 // reasons why the executable file might be covered by the GNU General
  29 // Public License.
  30
  31 /** @file parallel/random_number.h
  32  *  @brief Random number generator based on the Mersenne twister.
  33  *  This file is a GNU parallel extension to the Standard C++ Library.
  34  */
  35
  36 // Written by Johannes Singler.
  37
  38 #ifndef _GLIBCXX_PARALLEL_RANDOM_NUMBER_H
  39 #define _GLIBCXX_PARALLEL_RANDOM_NUMBER_H 1
  40
  41 #include <parallel/types.h>
  42
  43 namespace __gnu_parallel
  44 {
  45   // XXX use tr1 random number.
  46   // http://www.math.keio.ac.jp/matumoto/emt.html
  47   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
  48            int s, UIntType b, int t, UIntType c, int l, UIntType val>
  49     class mersenne_twister
  50     {
  51     public:
  52       typedef UIntType result_type;
  53       static const int word_size = w;
  54       static const int state_size = n;
  55       static const int shift_size = m;
  56       static const int mask_bits = r;
  57       static const UIntType parameter_a = a;
  58       static const int output_u = u;
  59       static const int output_s = s;
  60       static const UIntType output_b = b;
  61       static const int output_t = t;
  62       static const UIntType output_c = c;
  63       static const int output_l = l;
  64
  65       static const bool has_fixed_range = false;
  66
  67       mersenne_twister() { seed(); }
  68
  69 #if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
  70       // Work around overload resolution problem (Gennadiy E. Rozental)
  71       explicit
  72       mersenne_twister(const UIntType& value)
  73 #else
  74       explicit
  75       mersenne_twister(UIntType value)
  76 #endif
  77       { seed(value); }
  78
  79       template<typename It>
  80         mersenne_twister(It& first, It last)
  81         { seed(first,last); }
  82
  83       template<typename Generator>
  84         explicit
  85         mersenne_twister(Generator & gen)
  86         { seed(gen); }
  87
  88       // compiler-generated copy ctor and assignment operator are fine
  89
  90       void
  91       seed()
  92       { seed(UIntType(5489)); }
  93
  94 #if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
  95       // Work around overload resolution problem (Gennadiy E. Rozental)
  96       void
  97       seed(const UIntType& value)
  98 #else
  99       void
 100       seed(UIntType value)
 101 #endif
 102       {
 103         // New seeding algorithm from
 104       // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
 105         // In the previous versions, MSBs of the seed affected only MSBs of the
 106         // state x[].
 107         const UIntType mask = ~0u;
 108         x[0] = value & mask;
 109         for (i = 1; i < n; ++i)
 110           {
 111             // See Knuth "The Art of Computer Programming" Vol. 2,
 112             // 3rd ed., page 106
 113             x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
 114           }
 115       }
 116
 117       // For GCC, moving this function out-of-line prevents inlining, which may
 118       // reduce overall object code size.  However, MSVC does not grok
 119       // out-of-line definitions of member function templates.
 120       template<typename Generator>
 121         void
 122         seed(Generator & gen)
 123         {
 124           // I could have used std::generate_n, but it takes "gen" by value
 125           for (int j = 0; j < n; ++j)
 126             x[j] = gen();
 127           i = n;
 128         }
 129
 130       template<typename It>
 131         void
 132         seed(It& first, It last)
 133         {
 134           int j;
 135           for (j = 0; j < n && first != last; ++j, ++first)
 136             x[j] = *first;
 137           i = n;
 138           /*    if (first == last && j < n)
 139                 throw std::invalid_argument("mersenne_twister::seed");*/
 140         }
 141
 142       result_type
 143       min() const
 144       { return 0; }
 145
 146       result_type
 147       max() const
 148       {
 149         // avoid "left shift count >= with of type" warning
 150         result_type res = 0;
 151         for (int i = 0; i < w; ++i)
 152           res |= (1u << i);
 153         return res;
 154       }
 155
 156       result_type
 157       operator()();
 158
 159       static bool
 160       validation(result_type v)
 161       { return val == v; }
 162
 163 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
 164
 165       friend bool
 166       operator==(const mersenne_twister& x, const mersenne_twister& y)
 167       {
 168         for (int j = 0; j < state_size; ++j)
 169           if (x.compute(j) != y.compute(j))
 170             return false;
 171         return true;
 172       }
 173
 174       friend bool
 175       operator!=(const mersenne_twister& x, const mersenne_twister& y)
 176       { return !(x == y); }
 177 #else
 178       // Use a member function; Streamable concept not supported.
 179       bool
 180       operator==(const mersenne_twister& rhs) const
 181       {
 182         for (int j = 0; j < state_size; ++j)
 183           if (compute(j) != rhs.compute(j))
 184             return false;
 185         return true;
 186       }
 187
 188       bool
 189       operator!=(const mersenne_twister& rhs) const
 190       { return !(*this == rhs); }
 191 #endif
 192
 193     private:
 194       // returns x(i-n+index), where index is in 0..n-1
 195       UIntType
 196       compute(unsigned int index) const
 197       {
 198         // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers
 199         return x[ (i + n + index) % (2*n) ];
 200       }
 201
 202       void
 203       twist(int block);
 204
 205       // state representation: next output is o(x(i))
 206       //  x[0]  ... x[k] x[k+1] ... x[n-1]     x[n]     ... x[2*n-1] represents
 207       //  x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)]
 208       // The goal is to always have x(i-n) ... x(i-1) available for
 209       // operator== and save/restore.
 210
 211       UIntType x[2*n];
 212       int i;
 213   };
 214
 215 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 216   //  A definition is required even for integral static constants
 217   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 218            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 219     const bool
 220     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range;
 221
 222   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 223            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 224     const int
 225     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::state_size;
 226
 227   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 228            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 229     const int
 230     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::shift_size;
 231
 232   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 233            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 234     const int
 235     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::mask_bits;
 236
 237   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 238            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 239     const UIntType
 240     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::parameter_a;
 241
 242   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 243            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 244     const int
 245     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_u;
 246
 247   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 248            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 249     const int
 250     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_s;
 251
 252   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 253            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 254     const UIntType
 255     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_b;
 256
 257   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 258            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 259     const int
 260     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_t;
 261
 262   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 263            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 264     const UIntType
 265     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_c;
 266
 267   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 268            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 269     const int
 270     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_l;
 271 #endif
 272
 273   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 274            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 275   void
 276   mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block)
 277   {
 278     const UIntType upper_mask = (~0u) << r;
 279     const UIntType lower_mask = ~upper_mask;
 280
 281     if (block == 0)
 282       {
 283         for (int j = n; j < 2*n; ++j)
 284           {
 285             UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);
 286             x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
 287           }
 288       }
 289     else if (block == 1)
 290       {
 291         // split loop to avoid costly modulo operations
 292         {  // extra scope for MSVC brokenness w.r.t. for scope
 293           for (int j = 0; j < n-m; ++j)
 294             {
 295               UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
 296               x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);
 297             }
 298         }
 299
 300         for (int j = n-m; j < n-1; ++j)
 301           {
 302             UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
 303             x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
 304           }
 305         // last iteration
 306         UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);
 307         x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);
 308         i = 0;
 309       }
 310   }
 311
 312   template<typename UIntType, int w, int n, int m, int r, UIntType a, int u,
 313            int s, UIntType b, int t, UIntType c, int l, UIntType val>
 314     inline
 315     typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type
 316     mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()()
 317     {
 318       if (i == n)
 319         twist(0);
 320       else if (i >= 2*n)
 321         twist(1);
 322       // Step 4
 323       UIntType z = x[i];
 324       ++i;
 325       z ^= (z >> u);
 326       z ^= ((z << s) & b);
 327       z ^= ((z << t) & c);
 328       z ^= (z >> l);
 329       return z;
 330     }
 331
 332
 333   typedef mersenne_twister<uint32,32,351,175,19,0xccab8ee7,11,
 334                            7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b;
 335
 336   // validation by experiment from mt19937.c
 337   typedef mersenne_twister<uint32,32,624,397,31,0x9908b0df,11,
 338                            7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
 339
 340   /** @brief Random number generator, based on the Mersenne twister. */
 341   class random_number
 342   {
 343   private:
 344     mt19937 mt;
 345     uint64 supremum, RAND_SUP;
 346     double supremum_reciprocal, RAND_SUP_REC;
 347
 348     uint64 cache;  /* assumed to be twice as long as the usual random number */
 349     int bits_left; /* bit results */
 350
 351     static uint32
 352     scale_down(uint64 x,
 353 #if _GLIBCXX_SCALE_DOWN_FPU
 354                uint64 /*supremum*/, double supremum_reciprocal)
 355 #else
 356                uint64 supremum, double /*supremum_reciprocal*/)
 357 #endif
 358         {
 359 #if _GLIBCXX_SCALE_DOWN_FPU
 360           return (uint32)(x * supremum_reciprocal);
 361 #else
 362           return static_cast<uint32>(x % supremum);
 363 #endif
 364         }
 365
 366   public:
 367     /** @brief Default constructor. Seed with 0. */
 368     random_number()
 369     : mt(0), supremum(0x100000000ULL),
 370       RAND_SUP(1ULL << (sizeof(uint32) * 8)),
 371       supremum_reciprocal((double)supremum / (double)RAND_SUP),
 372       RAND_SUP_REC(1.0 / (double)RAND_SUP),
 373       cache(0), bits_left(0) { }
 374
 375     /** @brief Constructor.
 376      *  @param seed Random seed.
 377      *  @param supremum Generate integer random numbers in the
 378      *                  interval @c [0,supremum). */
 379     random_number(uint32 seed, uint64 supremum = 0x100000000ULL)
 380     : mt(seed), supremum(supremum),
 381       RAND_SUP(1ULL << (sizeof(uint32) * 8)),
 382       supremum_reciprocal((double)supremum / (double)RAND_SUP),
 383       RAND_SUP_REC(1.0 / (double)RAND_SUP),
 384       cache(0), bits_left(0) { }
 385
 386     /** @brief Generate unsigned random 32-bit integer. */
 387     uint32
 388     operator()()
 389     { return scale_down(mt(), supremum, supremum_reciprocal); }
 390
 391     /** @brief Generate unsigned random 32-bit integer in the
 392         interval @c [0,local_supremum). */
 393     uint32
 394     operator()(uint64 local_supremum)
 395     {
 396       return scale_down(mt(), local_supremum,
 397                         (double)local_supremum * RAND_SUP_REC);
 398     }
 399
 400     /** @brief Set the random seed.
 401      *  @param seed to set. */
 402     void
 403     set_seed(uint32 seed)
 404     {
 405       mt.seed(seed);
 406       cache = mt();
 407       bits_left = 32;
 408     }
 409
 410     /** @brief Generate a number of random bits, compile-time parameter. */
 411     template<int bits>
 412       unsigned long
 413       genrand_bits()
 414       {
 415         unsigned long res = cache & ((1 << bits) - 1);
 416         cache = cache >> bits;
 417         bits_left -= bits;
 418         if (bits_left < 32)
 419           {
 420             cache |= (((uint64)mt()) << bits_left);
 421             bits_left += 32;
 422           }
 423         return res;
 424       }
 425
 426     /** @brief Generate a number of random bits, run-time parameter.
 427      *  @param bits Number of bits to generate. */
 428     unsigned long
 429     genrand_bits(int bits)
 430     {
 431       unsigned long res = cache & ((1 << bits) - 1);
 432       cache = cache >> bits;
 433       bits_left -= bits;
 434       if (bits_left < 32)
 435         {
 436           cache |= (((uint64)mt()) << bits_left);
 437           bits_left += 32;
 438         }
 439       return res;
 440     }
 441 };
 442
 443 } // namespace __gnu_parallel
 444
 445 #endif