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>
  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,
 120                       __w == std::numeric_limits<_Tp>::digits>::__calc(__x); }
 121
 122     // Selector to return the maximum value possible that will fit in
 123     // @p __w bits of @p _Tp.
 124     template<typename _Tp, _Tp __w, bool>
 125       struct _Max_w
 126       {
 127         static _Tp
 128         __value()
 129         { return (_Tp(1) << __w) - 1; }
 130       };
 131
 132     template<typename _Tp, _Tp __w>
 133       struct _Max_w<_Tp, __w, true>
 134       {
 135         static _Tp
 136         __value()
 137         { return std::numeric_limits<_Tp>::max(); }
 138       };
 139
 140   } // namespace _Private
 141
 142
 143   /**
 144    * Constructs the LCR engine with integral seed @p __x0.
 145    */
 146   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 147     linear_congruential<_UIntType, __a, __c, __m>::
 148     linear_congruential(unsigned long __x0)
 149     { this->seed(__x0); }
 150
 151   /**
 152    * Constructs the LCR engine with seed generated from @p __g.
 153    */
 154   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 155     template<class _Gen>
 156       linear_congruential<_UIntType, __a, __c, __m>::
 157       linear_congruential(_Gen& __g)
 158       { this->seed(__g); }
 159
 160   /**
 161    * Seeds the LCR with integral value @p __x0, adjusted so that the
 162    * ring identity is never a member of the convergence set.
 163    */
 164   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 165     void
 166     linear_congruential<_UIntType, __a, __c, __m>::
 167     seed(unsigned long __x0)
 168     {
 169       if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
 170           && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 171         _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
 172       else
 173         _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
 174     }
 175
 176   /**
 177    * Seeds the LCR engine with a value generated by @p __g.
 178    */
 179   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 180     template<class _Gen>
 181       void
 182       linear_congruential<_UIntType, __a, __c, __m>::
 183       seed(_Gen& __g, false_type)
 184       {
 185         _UIntType __x0 = __g();
 186         if ((_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0)
 187             && (_Private::__mod<_UIntType, 1, 0, __m>(__x0) == 0))
 188           _M_x = _Private::__mod<_UIntType, 1, 0, __m>(1);
 189         else
 190           _M_x = _Private::__mod<_UIntType, 1, 0, __m>(__x0);
 191       }
 192
 193   /**
 194    * Returns a value that is less than or equal to all values potentially
 195    * returned by operator(). The return value of this function does not
 196    * change during the lifetime of the object..
 197    *
 198    * The minumum depends on the @p __c parameter: if it is zero, the
 199    * minimum generated must be > 0, otherwise 0 is allowed.
 200    */
 201   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 202     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 203     linear_congruential<_UIntType, __a, __c, __m>::
 204     min() const
 205     { return (_Private::__mod<_UIntType, 1, 0, __m>(__c) == 0) ? 1 : 0; }
 206
 207   /**
 208    * Gets the maximum possible value of the generated range.
 209    *
 210    * For a linear congruential generator, the maximum is always @p __m - 1.
 211    */
 212   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 213     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 214     linear_congruential<_UIntType, __a, __c, __m>::
 215     max() const
 216     { return (__m == 0) ? std::numeric_limits<_UIntType>::max() : (__m - 1); }
 217
 218   /**
 219    * Gets the next generated value in sequence.
 220    */
 221   template<class _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
 222     typename linear_congruential<_UIntType, __a, __c, __m>::result_type
 223     linear_congruential<_UIntType, __a, __c, __m>::
 224     operator()()
 225     {
 226       _M_x = _Private::__mod<_UIntType, __a, __c, __m>(_M_x);
 227       return _M_x;
 228     }
 229
 230
 231   template<class _UIntType, int __w, int __n, int __m, int __r,
 232            _UIntType __a, int __u, int __s,
 233            _UIntType __b, int __t, _UIntType __c, int __l>
 234     void
 235     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 236                      __b, __t, __c, __l>::
 237     seed(unsigned long __value)
 238     {
 239       _M_x[0] = _Private::__mod_w<_UIntType, __w>(__value);
 240
 241       for (int __i = 1; __i < state_size; ++__i)
 242         {
 243           _UIntType __x = _M_x[__i - 1];
 244           __x ^= __x >> (__w - 2);
 245           __x *= 1812433253ul;
 246           __x += __i;
 247           _M_x[__i] = _Private::__mod_w<_UIntType, __w>(__x);
 248         }
 249       _M_p = state_size;
 250     }
 251
 252   template<class _UIntType, int __w, int __n, int __m, int __r,
 253            _UIntType __a, int __u, int __s,
 254            _UIntType __b, int __t, _UIntType __c, int __l>
 255     template<class _Gen>
 256       void
 257       mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 258                        __b, __t, __c, __l>::
 259       seed(_Gen& __gen, false_type)
 260       {
 261         for (int __i = 0; __i < state_size; ++__i)
 262           _M_x[__i] = _Private::__mod_w<_UIntType, __w>(__gen());
 263         _M_p = state_size;
 264       }
 265
 266   template<class _UIntType, int __w, int __n, int __m, int __r,
 267            _UIntType __a, int __u, int __s,
 268            _UIntType __b, int __t, _UIntType __c, int __l>
 269     typename
 270     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 271                      __b, __t, __c, __l>::result_type
 272     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 273                      __b, __t, __c, __l>::
 274     max() const
 275     {
 276       using _Private::_Max_w;
 277       using std::numeric_limits;
 278       return _Max_w<_UIntType, __w,
 279                     __w == numeric_limits<_UIntType>::digits>::__value();
 280     }
 281
 282   template<class _UIntType, int __w, int __n, int __m, int __r,
 283            _UIntType __a, int __u, int __s,
 284            _UIntType __b, int __t, _UIntType __c, int __l>
 285     typename
 286     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 287                      __b, __t, __c, __l>::result_type
 288     mersenne_twister<_UIntType, __w, __n, __m, __r, __a, __u, __s,
 289                      __b, __t, __c, __l>::
 290     operator()()
 291     {
 292       // Reload the vector - cost is O(n) amortized over n calls.
 293       if (_M_p >= state_size)
 294         {
 295           const _UIntType __upper_mask = (~_UIntType()) << __r;
 296           const _UIntType __lower_mask = ~__upper_mask;
 297
 298           for (int __k = 0; __k < (__n - __m); ++__k)
 299             {
 300               _UIntType __y = ((_M_x[__k] & __upper_mask)
 301                                |(_M_x[__k + 1] & __lower_mask));
 302               _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
 303                            ^ ((__y & 0x01) ? __a : 0));
 304             }
 305
 306           for (int __k = (__n - __m); __k < (__n - 1); ++__k)
 307             {
 308               _UIntType __y = ((_M_x[__k] & __upper_mask)
 309                                | (_M_x[__k + 1] & __lower_mask));
 310               _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
 311                            ^ ((__y & 0x01) ? __a : 0));
 312             }
 313
 314           _M_p = 0;
 315         }
 316
 317       // Calculate o(x(i)).
 318       result_type __z = _M_x[_M_p++];
 319       __z ^= (__z >> __u);
 320       __z ^= (__z << __s) & __b;
 321       __z ^= (__z << __t) & __c;
 322       __z ^= (__z >> __l);
 323
 324       return __z;
 325     }
 326
 327
 328   template<typename _IntType, _IntType __m, int __s, int __r>
 329     void
 330     subtract_with_carry<_IntType, __m, __s, __r>::
 331     seed(unsigned long __value)
 332     {
 333       std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
 334         __lcg(__value);
 335
 336       for (int __i = 0; __i < long_lag; ++__i)
 337         _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(__lcg());
 338
 339       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 340       _M_p = 0;
 341     }
 342
 343   //
 344   // This implementation differs from the tr1 spec because the tr1 spec refused
 345   // to make any sense to me:  the exponent of the factor in the spec goes from
 346   // 1 to (n-1), but it would only make sense to me if it went from 0 to (n-1).
 347   //
 348   // This algorithm is still problematic because it can overflow left right and
 349   // center.
 350   //
 351   template<typename _IntType, _IntType __m, int __s, int __r>
 352     template<class _Gen>
 353     void
 354     subtract_with_carry<_IntType, __m, __s, __r>::
 355     seed(_Gen& __gen, false_type)
 356     {
 357       const int __n = (std::numeric_limits<_IntType>::digits + 31) / 32;
 358       for (int __i = 0; __i < long_lag; ++__i)
 359         {
 360           _M_x[__i] = 0;
 361           unsigned long __factor = 1;
 362           for (int __j = 0; __j < __n; ++__j)
 363             {
 364               _M_x[__i] += __gen() * __factor;
 365               __factor *= 0x80000000;
 366             }
 367           _M_x[__i] = _Private::__mod<_IntType, 1, 0, modulus>(_M_x[__i]);
 368         }
 369       _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
 370       _M_p = 0;
 371     }
 372
 373   template<typename _IntType, _IntType __m, int __s, int __r>
 374     typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
 375     subtract_with_carry<_IntType, __m, __s, __r>::
 376     operator()()
 377     {
 378       // Derive short lag index from current index.
 379       int __ps = _M_p - short_lag;
 380       if (__ps < 0)
 381         __ps += long_lag;
 382
 383       // Calculate new x(i) without overflow or division.
 384       _IntType __xi;
 385       if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
 386         {
 387           __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
 388           _M_carry = 0;
 389         }
 390       else
 391         {
 392           __xi = modulus - _M_x[_M_p] - _M_carry + _M_x[__ps];
 393           _M_carry = 1;
 394         }
 395       _M_x[_M_p++] = __xi;
 396
 397       // Adjust current index to loop around in ring buffer.
 398       if (_M_p >= long_lag)
 399         _M_p = 0;
 400
 401       return __xi;
 402     }
 403
 404
 405   template<class _UniformRandomNumberGenerator, int __p, int __r>
 406     typename discard_block<_UniformRandomNumberGenerator,
 407                            __p, __r>::result_type
 408     discard_block<_UniformRandomNumberGenerator, __p, __r>::
 409     operator()()
 410     {
 411       if (_M_n >= used_block)
 412         {
 413           while (_M_n < block_size)
 414             {
 415               _M_b();
 416               ++_M_n;
 417             }
 418           _M_n = 0;
 419         }
 420       ++_M_n;
 421       return _M_b();
 422     }
 423
 424 _GLIBCXX_END_NAMESPACE
 425 }