libstdc++-v3/include/tr1/random.tcc

   1 // random number generation (out of line) -*- C++ -*-
   2
   3 // Copyright (C) 2006 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
   7 // terms of the GNU General Public License as published by the
   8 // Free Software Foundation; either version 2, or (at your option)
   9 // any later version.
  10
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 // GNU General Public License for more details.
  15
  16 // You should have received a copy of the GNU General Public License along
  17 // with this library; see the file COPYING.  If not, write to the Free
  18 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
  19 // USA.
  20
  21 // As a special exception, you may use this file as part of a free software
  22 // library without restriction.  Specifically, if other files instantiate
  23 // templates or use macros or inline functions from this file, or you compile
  24 // this file and link it with other files to produce an executable, this
  25 // file does not by itself cause the resulting executable to be covered by
  26 // the GNU General Public License.  This exception does not however
  27 // invalidate any other reasons why the executable file might be covered by
  28 // the GNU General Public License.
  29
  30 #include <limits>
  31
  32 namespace std
  33 {
  34 _GLIBCXX_BEGIN_NAMESPACE(tr1)
  35
  36   /*
  37    * Implementation-space details.
  38    */
  39   namespace _Private
  40   {
  41     // General case for x = (ax + c) mod m -- use Schrage's algorithm to avoid
  42     // integer overflow.
  43     //
  44     // Because a and c are compile-time integral constants the compiler kindly
  45     // elides any unreachable paths.
  46     //
  47     // Preconditions:  a > 0, m > 0.
  48     //
  49     template<typename _Tp, _Tp a, _Tp c, _Tp m, bool _m_is_zero>
  50       struct Mod
  51       {
  52         static _Tp
  53         calc(_Tp x)
  54         {
  55           if (a == 1)
  56             x %= m;
  57           else
  58             {
  59               static const _Tp q = m / a;
  60               static const _Tp r = m % a;
  61
  62               _Tp t1 = a * (x % q);
  63               _Tp t2 = r * (x / q);
  64               if (t1 >= t2)
  65                 x = t1 - t2;
  66               else
  67                 x = m - t2 + t1;
  68             }
  69
  70           if (c != 0)
  71             {
  72               const _Tp d = m - x;
  73               if (d > c)
  74                 x += c;
  75               else
  76                 x = c - d;
  77             }
  78           return x;
  79         }
  80       };
  81
  82     // Special case for m==0 -- use unsigned integer overflow as modulo
  83     // operator.
  84     template<typename _Tp, _Tp a, _Tp c, _Tp m>
  85       struct Mod<_Tp, a, c, m, true>
  86       {
  87         static _Tp
  88         calc(_Tp x)
  89         { return a * x + c; }
  90       };
  91
  92     // Dispatch based on modulus value to prevent divide-by-zero compile-time
  93     // errors when m == 0.
  94     template<typename _Tp, _Tp a, _Tp c, _Tp m>
  95       inline _Tp
  96       mod(_Tp x)
  97       { return Mod<_Tp, a, c, m, m == 0>::calc(x); }
  98
  99     // Like the above, for a==1, c==0, in terms of w.
 100     template<typename _Tp, _Tp w, bool>
 101       struct Mod_w
 102       {
 103         static _Tp
 104         calc(_Tp x)
 105         { return x % (_Tp(1) << w); }
 106       };
 107
 108     template<typename _Tp, _Tp w>
 109       struct Mod_w<_Tp, w, true>
 110       {
 111         static _Tp
 112         calc(_Tp x)
 113         { return x; }
 114       };
 115
 116     template<typename _Tp, _Tp w>
 117       inline _Tp
 118       mod_w(_Tp x)
 119       { return Mod_w<_Tp, w, w == std::numeric_limits<_Tp>::digits>::calc(x); }
 120
 121     // Selector to return the maximum value possible that will fit in
 122     // @p w bits of @p _Tp.
 123     template<typename _Tp, _Tp w, bool>
 124       struct Max_w
 125       {
 126         static _Tp
 127         value()
 128         { return (_Tp(1) << w) - 1; }
 129       };
 130
 131     template<typename _Tp, _Tp w>
 132       struct Max_w<_Tp, w, true>
 133       {
 134         static _Tp
 135         value()
 136         { return std::numeric_limits<_Tp>::max(); }
 137       };
 138
 139   } // namespace _Private
 140
 141
 142   /**
 143    * Constructs the LCR engine with integral seed @p x0.
 144    */
 145   template<class UIntType, UIntType a, UIntType c, UIntType m>
 146     linear_congruential<UIntType, a, c, m>::
 147     linear_congruential(unsigned long x0)
 148     { this->seed(x0); }
 149
 150   /**
 151    * Constructs the LCR engine with seed generated from @p g.
 152    */
 153   template<class UIntType, UIntType a, UIntType c, UIntType m>
 154     template<class Gen>
 155       linear_congruential<UIntType, a, c, m>::
 156       linear_congruential(Gen& g)
 157       { this->seed(g); }
 158
 159   /**
 160    * Seeds the LCR with integral value @p x0, adjusted so that the
 161    * ring identity is never a member of the convergence set.
 162    */
 163   template<class UIntType, UIntType a, UIntType c, UIntType m>
 164     void
 165     linear_congruential<UIntType, a, c, m>::
 166     seed(unsigned long x0)
 167     {
 168       if ((_Private::mod<UIntType, 1, 0, m>(c) == 0)
 169           && (_Private::mod<UIntType, 1, 0, m>(x0) == 0))
 170         m_x = _Private::mod<UIntType, 1, 0, m>(1);
 171       else
 172         m_x = _Private::mod<UIntType, 1, 0, m>(x0);
 173     }
 174
 175   /**
 176    * Seeds the LCR engine with a value generated by @p g.
 177    */
 178   template<class UIntType, UIntType a, UIntType c, UIntType m>
 179     template<class Gen>
 180       void
 181       linear_congruential<UIntType, a, c, m>::
 182       seed(Gen& g, false_type)
 183       {
 184         UIntType x0 = g();
 185         if ((_Private::mod<UIntType, 1, 0, m>(c) == 0)
 186             && (_Private::mod<UIntType, 1, 0, m>(x0) == 0))
 187           m_x = _Private::mod<UIntType, 1, 0, m>(1);
 188         else
 189           m_x = _Private::mod<UIntType, 1, 0, m>(x0);
 190       }
 191
 192   /**
 193    * Returns a value that is less than or equal to all values potentially
 194    * returned by operator(). The return value of this function does not
 195    * change during the lifetime of the object..
 196    *
 197    * The minumum depends on the @p c parameter: if it is zero, the
 198    * minimum generated must be > 0, otherwise 0 is allowed.
 199    */
 200   template<class UIntType, UIntType a, UIntType c, UIntType m>
 201     typename linear_congruential<UIntType, a, c, m>::result_type
 202     linear_congruential<UIntType, a, c, m>::
 203     min() const
 204     { return (_Private::mod<UIntType, 1, 0, m>(c) == 0) ? 1 : 0; }
 205
 206   /**
 207    * Gets the maximum possible value of the generated range.
 208    *
 209    * For a linear congruential generator, the maximum is always @p m - 1.
 210    */
 211   template<class UIntType, UIntType a, UIntType c, UIntType m>
 212     typename linear_congruential<UIntType, a, c, m>::result_type
 213     linear_congruential<UIntType, a, c, m>::
 214     max() const
 215     { return (m == 0) ? std::numeric_limits<UIntType>::max() : (m - 1); }
 216
 217   /**
 218    * Gets the next generated value in sequence.
 219    */
 220   template<class UIntType, UIntType a, UIntType c, UIntType m>
 221     typename linear_congruential<UIntType, a, c, m>::result_type
 222     linear_congruential<UIntType, a, c, m>::
 223     operator()()
 224     {
 225       m_x = _Private::mod<UIntType, a, c, m>(m_x);
 226       return m_x;
 227     }
 228
 229
 230   template<class _UInt, int w, int n, int m, int r,
 231            _UInt a, int u, int s,
 232            _UInt b, int t, _UInt c, int l>
 233     void
 234     mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
 235     seed(unsigned long value)
 236     {
 237       _M_x[0] = _Private::mod_w<_UInt, w>(value);
 238
 239       for (int i = 1; i < n; ++i)
 240         {
 241           _UInt x = _M_x[i - 1];
 242           x ^= x >> (w - 2);
 243           x *= 1812433253ul;
 244           x += i;
 245           _M_x[i] = _Private::mod_w<_UInt, w>(x);
 246         }
 247       _M_p = n;
 248     }
 249
 250   template<class _UInt, int w, int n, int m, int r,
 251            _UInt a, int u, int s,
 252            _UInt b, int t, _UInt c, int l>
 253     template<class Gen>
 254       void
 255       mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
 256       seed(Gen& gen, false_type)
 257       {
 258         for (int i = 0; i < n; ++i)
 259           _M_x[i] = _Private::mod_w<_UInt, w>(gen());
 260         _M_p = n;
 261       }
 262
 263   template<class _UInt, int w, int n, int m, int r,
 264            _UInt a, int u, int s,
 265            _UInt b, int t, _UInt c, int l>
 266     typename
 267     mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::result_type
 268     mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
 269     max() const
 270     {
 271       using _Private::Max_w;
 272       using std::numeric_limits;
 273       return Max_w<_UInt, w, w == numeric_limits<_UInt>::digits>::value();
 274     }
 275
 276   template<class _UInt, int w, int n, int m, int r,
 277            _UInt a, int u, int s,
 278            _UInt b, int t, _UInt c, int l>
 279     typename
 280     mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::result_type
 281     mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
 282     operator()()
 283     {
 284       // Reload the vector - cost is O(n) amortized over n calls.
 285       if (_M_p >= n)
 286         {
 287           const _UInt upper_mask = (~_UInt()) << r;
 288           const _UInt lower_mask = ~upper_mask;
 289
 290           for (int k = 0; k < (n - m); ++k)
 291             {
 292               _UInt y = (_M_x[k] & upper_mask) | (_M_x[k + 1] & lower_mask);
 293               _M_x[k] = _M_x[k + m] ^ (y >> 1) ^ ((y & 0x01) ? a : 0);
 294             }
 295
 296           for (int k = (n - m); k < (n - 1); ++k)
 297             {
 298               _UInt y = (_M_x[k] & upper_mask) | (_M_x[k + 1] & lower_mask);
 299               _M_x[k] = _M_x[k + (m - n)] ^ (y >> 1) ^ ((y & 0x01) ? a : 0);
 300             }
 301
 302           _M_p = 0;
 303         }
 304
 305       // Calculate o(x(i)).
 306       result_type z = _M_x[_M_p++];
 307       z ^= (z >> u);
 308       z ^= (z << s) & b;
 309       z ^= (z << t) & c;
 310       z ^= (z >> l);
 311
 312       return z;
 313     }
 314
 315
 316   template<typename _IntType, _IntType m, int s, int r>
 317     void
 318     subtract_with_carry<_IntType, m, s, r>::
 319     seed(_IntType __value)
 320     {
 321       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
 322         lcg(__value);
 323
 324       for (int i = 0; i < long_lag; ++i)
 325         _M_x[i] = _Private::mod<_IntType, 1, 0, modulus>(lcg());
 326
 327       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 328       _M_p = 0;
 329     }
 330
 331   //
 332   // This implementation differs from the tr1 spec because the tr1 spec refused
 333   // to make any sense to me:  the exponent of the factor in the spec goes from
 334   // 1 to (n-1), but it would only make sense to me if it went from 0 to (n-1).
 335   //
 336   // This algorithm is still problematic because it can overflow left right and
 337   // center.
 338   //
 339   template<typename _IntType, _IntType __m, int __s, int __r>
 340     template<class Gen>
 341     void
 342     subtract_with_carry<_IntType, __m, __s, __r>::
 343     seed(Gen& gen, false_type)
 344     {
 345       const int n = (std::numeric_limits<_IntType>::digits + 31) / 32;
 346       for (int i = 0; i < long_lag; ++i)
 347         {
 348           _M_x[i] = 0;
 349           unsigned long factor = 1;
 350           for (int j = 0; j < n; ++j)
 351             {
 352               _M_x[i] += gen() * factor;
 353               factor *= 0x80000000;
 354             }
 355           _M_x[i] = _Private::mod<_IntType, 1, 0, modulus>(_M_x[i]);
 356         }
 357       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 358       _M_p = 0;
 359     }
 360
 361   template<typename _IntType, _IntType __m, int __s, int __r>
 362     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
 363     subtract_with_carry<_IntType, __m, __s, __r>::
 364     operator()()
 365     {
 366       // Derive short lag index from current index.
 367       int ps = _M_p - short_lag;
 368       if (ps < 0)
 369         ps += long_lag;
 370
 371       // Calculate new x(i) without overflow or division.
 372       _IntType xi;
 373       if (_M_x[ps] >= _M_x[_M_p] + _M_carry)
 374         {
 375           xi = _M_x[ps] - _M_x[_M_p] - _M_carry;
 376           _M_carry = 0;
 377         }
 378       else
 379         {
 380           xi = modulus - _M_x[_M_p] - _M_carry + _M_x[ps];
 381           _M_carry = 1;
 382         }
 383       _M_x[_M_p++] = xi;
 384
 385       // Adjust current index to loop around in ring buffer.
 386       if (_M_p >= long_lag)
 387         _M_p = 0;
 388
 389       return xi;
 390     }
 391
 392
 393   template<class _E, int __p, int __r>
 394     typename discard_block<_E, __p, __r>::result_type
 395     discard_block<_E, __p, __r>::
 396     operator()()
 397     {
 398       if (_M_n >= used_block)
 399         {
 400           while (_M_n < block_size)
 401             {
 402               _M_b();
 403               ++_M_n;
 404             }
 405           _M_n = 0;
 406         }
 407       ++_M_n;
 408       return _M_b();
 409     }
 410
 411 _GLIBCXX_END_NAMESPACE
 412 }