comparison DEPENDENCIES/generic/include/boost/random/mersenne_twister.hpp @ 101:c530137014c0

Update Boost headers (1.58.0)

  * accompanying file LICENSE_1_0.txt or copy at
  * http://www.boost.org/LICENSE_1_0.txt)
  *
  * See http://www.boost.org for most recent version including documentation.
  *
- * $Id: mersenne_twister.hpp 74867 2011-10-09 23:13:31Z steven_watanabe $
+ * $Id$
  *
  * Revision history
+ *  2013-10-14  fixed some warnings with Wshadow (mgaunard)
  *  2001-02-18  moved to individual header files
  */
 
 #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP
 #define BOOST_RANDOM_MERSENNE_TWISTER_HPP

 #include <boost/random/detail/config.hpp>
 #include <boost/random/detail/ptr_helper.hpp>
 #include <boost/random/detail/seed.hpp>
 #include <boost/random/detail/seed_impl.hpp>
 #include <boost/random/detail/generator_seed_seq.hpp>
+#include <boost/random/detail/polynomial.hpp>
+
+#include <boost/random/detail/disable_warnings.hpp>
 
 namespace boost {
 namespace random {
 
 /**

  *
  *  @blockquote
  *  "Mersenne Twister: A 623-dimensionally equidistributed uniform
  *  pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura,
  *  ACM Transactions on Modeling and Computer Simulation: Special Issue on
- *  Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
+ *  Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
  *  @endblockquote
  *
  * @xmlnote
  * The boost variant has been implemented from scratch and does not
  * derive from or use mt19937.c provided on the above WWW site. However, it
  * was verified that both produce identical output.
  * @endxmlnote
  *
  * The seeding from an integer was changed in April 2005 to address a
  * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html">weakness</a>.
- * 
+ *
  * The quality of the generator crucially depends on the choice of the
  * parameters.  User code should employ one of the sensibly parameterized
  * generators such as \mt19937 instead.
  *
  * The generator requires considerable amounts of memory for the storage of

     BOOST_STATIC_CONSTANT(std::size_t, tempering_t = t);
     BOOST_STATIC_CONSTANT(UIntType, tempering_c = c);
     BOOST_STATIC_CONSTANT(std::size_t, tempering_l = l);
     BOOST_STATIC_CONSTANT(UIntType, initialization_multiplier = f);
     BOOST_STATIC_CONSTANT(UIntType, default_seed = 5489u);
-  
+
     // backwards compatibility
     BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);
     BOOST_STATIC_CONSTANT(std::size_t, output_u = u);
     BOOST_STATIC_CONSTANT(std::size_t, output_s = s);
     BOOST_STATIC_CONSTANT(UIntType, output_b = b);
     BOOST_STATIC_CONSTANT(std::size_t, output_t = t);
     BOOST_STATIC_CONSTANT(UIntType, output_c = c);
     BOOST_STATIC_CONSTANT(std::size_t, output_l = l);
-    
+
     // old Boost.Random concept requirements
     BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
 
 
     /**

      * (i + f * (x(i-1) xor (x(i-1) rshift w-2))) mod 2<sup>w</sup>
      * for i = 1 .. n-1. x(n) is the first value to be returned by operator().
      */
     BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(mersenne_twister_engine, UIntType, value)
     {
-        // New seeding algorithm from 
+        // New seeding algorithm from
         // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
         // In the previous versions, MSBs of the seed affected only MSBs of the
         // state x[].
         const UIntType mask = (max)();
         x[0] = value & mask;
         for (i = 1; i < n; i++) {
             // See Knuth "The Art of Computer Programming"
             // Vol. 2, 3rd ed., page 106
             x[i] = (f * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
         }
-    }
-    
+
+        normalize_state();
+    }
+
     /**
      * Seeds a mersenne_twister_engine using values produced by seq.generate().
      */
     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mersenne_twister_engine, SeeqSeq, seq)
     {
         detail::seed_array_int<w>(seq, x);
         i = n;
 
-        // fix up the state if it's all zeroes.
-        if((x[0] & (~static_cast<UIntType>(0) << r)) == 0) {
-            for(std::size_t j = 1; j < n; ++j) {
-                if(x[j] != 0) return;
-            }
-            x[0] = static_cast<UIntType>(1) << (w-1);
-        }
+        normalize_state();
     }
 
     /** Sets the state of the generator using values from an iterator range. */
     template<class It>
     void seed(It& first, It last)
     {
         detail::fill_array_int<w>(first, last, x);
         i = n;
 
-        // fix up the state if it's all zeroes.
-        if((x[0] & (~static_cast<UIntType>(0) << r)) == 0) {
-            for(std::size_t j = 1; j < n; ++j) {
-                if(x[j] != 0) return;
-            }
-            x[0] = static_cast<UIntType>(1) << (w-1);
-        }
-    }
-  
+        normalize_state();
+    }
+
     /** Returns the smallest value that the generator can produce. */
     static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
     { return 0; }
     /** Returns the largest value that the generator can produce. */
     static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
     { return boost::low_bits_mask_t<w>::sig_bits; }
-    
+
     /** Produces the next value of the generator. */
     result_type operator()();
 
     /** Fills a range with random values */
     template<class Iter>

      * }
      * @endcode
      */
     void discard(boost::uintmax_t z)
     {
-        for(boost::uintmax_t j = 0; j < z; ++j) {
-            (*this)();
+#ifndef BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD
+#define BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD 10000000
+#endif
+        if(z > BOOST_RANDOM_MERSENNE_TWISTER_DISCARD_THRESHOLD) {
+            discard_many(z);
+        } else {
+            for(boost::uintmax_t j = 0; j < z; ++j) {
+                (*this)();
+            }
         }
     }
 
 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
     /** Writes a mersenne_twister_engine to a @c std::ostream */

                const mersenne_twister_engine& mt)
     {
         mt.print(os);
         return os;
     }
-    
+
     /** Reads a mersenne_twister_engine from a @c std::istream */
     template<class CharT, class Traits>
     friend std::basic_istream<CharT,Traits>&
     operator>>(std::basic_istream<CharT,Traits>& is,
                mersenne_twister_engine& mt)

 
     /**
      * Returns true if the two generators are in the same state,
      * and will thus produce identical sequences.
      */
-    friend bool operator==(const mersenne_twister_engine& x,
-                           const mersenne_twister_engine& y)
-    {
-        if(x.i < y.i) return x.equal_imp(y);
-        else return y.equal_imp(x);
-    }
-    
+    friend bool operator==(const mersenne_twister_engine& x_,
+                           const mersenne_twister_engine& y_)
+    {
+        if(x_.i < y_.i) return x_.equal_imp(y_);
+        else return y_.equal_imp(x_);
+    }
+
     /**
      * Returns true if the two generators are in different states.
      */
-    friend bool operator!=(const mersenne_twister_engine& x,
-                           const mersenne_twister_engine& y)
-    { return !(x == y); }
+    friend bool operator!=(const mersenne_twister_engine& x_,
+                           const mersenne_twister_engine& y_)
+    { return !(x_ == y_); }
 
 private:
     /// \cond show_private
 
     void twist();

             y0 = y1;
         }
     }
 
     /**
+     * Converts an arbitrary array into a valid generator state.
+     * First we normalize x[0], so that it contains the same
+     * value we would get by running the generator forwards
+     * and then in reverse.  (The low order r bits are redundant).
+     * Then, if the state consists of all zeros, we set the
+     * high order bit of x[0] to 1.  This function only needs to
+     * be called by seed, since the state transform preserves
+     * this relationship.
+     */
+    void normalize_state()
+    {
+        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
+        const UIntType lower_mask = ~upper_mask;
+        UIntType y0 = x[m-1] ^ x[n-1];
+        if(y0 & (static_cast<UIntType>(1) << (w-1))) {
+            y0 = ((y0 ^ a) << 1) | 1;
+        } else {
+            y0 = y0 << 1;
+        }
+        x[0] = (x[0] & upper_mask) | (y0 & lower_mask);
+
+        // fix up the state if it's all zeroes.
+        for(std::size_t j = 0; j < n; ++j) {
+            if(x[j] != 0) return;
+        }
+        x[0] = static_cast<UIntType>(1) << (w-1);
+    }
+
+    /**
      * Given a pointer to the last element of the rewind array,
      * and the current size of the rewind array, finds an element
      * relative to the next available slot in the rewind array.
      */
     UIntType

         } else {
             return *(last - (n - 1 - index));
         }
     }
 
+    /**
+     * Optimized algorithm for large jumps.
+     *
+     * Hiroshi Haramoto, Makoto Matsumoto, and Pierre L'Ecuyer. 2008.
+     * A Fast Jump Ahead Algorithm for Linear Recurrences in a Polynomial
+     * Space. In Proceedings of the 5th international conference on
+     * Sequences and Their Applications (SETA '08).
+     * DOI=10.1007/978-3-540-85912-3_26
+     */
+    void discard_many(boost::uintmax_t z)
+    {
+        // Compute the minimal polynomial, phi(t)
+        // This depends only on the transition function,
+        // which is constant.  The characteristic
+        // polynomial is the same as the minimal
+        // polynomial for a maximum period generator
+        // (which should be all specializations of
+        // mersenne_twister.)  Even if it weren't,
+        // the characteristic polynomial is guaranteed
+        // to be a multiple of the minimal polynomial,
+        // which is good enough.
+        detail::polynomial phi = get_characteristic_polynomial();
+
+        // calculate g(t) = t^z % phi(t)
+        detail::polynomial g = mod_pow_x(z, phi);
+
+        // h(s_0, t) = \sum_{i=0}^{2k-1}o(s_i)t^{2k-i-1}
+        detail::polynomial h;
+        const std::size_t num_bits = w*n - r;
+        for(std::size_t j = 0; j < num_bits * 2; ++j) {
+            // Yes, we're advancing the generator state
+            // here, but it doesn't matter because
+            // we're going to overwrite it completely
+            // in reconstruct_state.
+            if(i >= n) twist();
+            h[2*num_bits - j - 1] = x[i++] & UIntType(1);
+        }
+        // g(t)h(s_0, t)
+        detail::polynomial gh = g * h;
+        detail::polynomial result;
+        for(std::size_t j = 0; j <= num_bits; ++j) {
+            result[j] = gh[2*num_bits - j - 1];
+        }
+        reconstruct_state(result);
+    }
+    static detail::polynomial get_characteristic_polynomial()
+    {
+        const std::size_t num_bits = w*n - r;
+        detail::polynomial helper;
+        helper[num_bits - 1] = 1;
+        mersenne_twister_engine tmp;
+        tmp.reconstruct_state(helper);
+        // Skip the first num_bits elements, since we
+        // already know what they are.
+        for(std::size_t j = 0; j < num_bits; ++j) {
+            if(tmp.i >= n) tmp.twist();
+            if(j == num_bits - 1)
+                assert((tmp.x[tmp.i] & 1) == 1);
+            else
+                assert((tmp.x[tmp.i] & 1) == 0);
+            ++tmp.i;
+        }
+        detail::polynomial phi;
+        phi[num_bits] = 1;
+        detail::polynomial next_bits = tmp.as_polynomial(num_bits);
+        for(std::size_t j = 0; j < num_bits; ++j) {
+            int val = next_bits[j] ^ phi[num_bits-j-1];
+            phi[num_bits-j-1] = val;
+            if(val) {
+                for(std::size_t k = j + 1; k < num_bits; ++k) {
+                    phi[num_bits-k-1] ^= next_bits[k-j-1];
+                }
+            }
+        }
+        return phi;
+    }
+    detail::polynomial as_polynomial(std::size_t size) {
+        detail::polynomial result;
+        for(std::size_t j = 0; j < size; ++j) {
+            if(i >= n) twist();
+            result[j] = x[i++] & UIntType(1);
+        }
+        return result;
+    }
+    void reconstruct_state(const detail::polynomial& p)
+    {
+        const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
+        const UIntType lower_mask = ~upper_mask;
+        const std::size_t num_bits = w*n - r;
+        for(std::size_t j = num_bits - n + 1; j <= num_bits; ++j)
+            x[j % n] = p[j];
+
+        UIntType y0 = 0;
+        for(std::size_t j = num_bits + 1; j >= n - 1; --j) {
+            UIntType y1 = x[j % n] ^ x[(j + m) % n];
+            if(p[j - n + 1])
+                y1 = (y1 ^ a) << UIntType(1) | UIntType(1);
+            else
+                y1 = y1 << UIntType(1);
+            x[(j + 1) % n] = (y0 & upper_mask) | (y1 & lower_mask);
+            y0 = y1;
+        }
+        i = 0;
+    }
+
     /// \endcond
 
     // state representation: next output is o(x(i))
     //   x[0]  ... x[k] x[k+1] ... x[n-1]   represents
     //  x(i-k) ... x(i) x(i+1) ... x(i-k+n-1)
 
-    UIntType x[n]; 
+    UIntType x[n];
     std::size_t i;
 };
 
 /// \cond show_private
 

  *  @blockquote
  *  "Mersenne Twister: A 623-dimensionally equidistributed
  *  uniform pseudo-random number generator", Makoto Matsumoto
  *  and Takuji Nishimura, ACM Transactions on Modeling and
  *  Computer Simulation: Special Issue on Uniform Random Number
- *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
+ *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
  *  @endblockquote
  */
 typedef mersenne_twister_engine<uint32_t,32,351,175,19,0xccab8ee7,
     11,0xffffffff,7,0x31b6ab00,15,0xffe50000,17,1812433253> mt11213b;
 

  *  @blockquote
  *  "Mersenne Twister: A 623-dimensionally equidistributed
  *  uniform pseudo-random number generator", Makoto Matsumoto
  *  and Takuji Nishimura, ACM Transactions on Modeling and
  *  Computer Simulation: Special Issue on Uniform Random Number
- *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
+ *  Generation, Vol. 8, No. 1, January 1998, pp. 3-30.
  *  @endblockquote
  */
 typedef mersenne_twister_engine<uint32_t,32,624,397,31,0x9908b0df,
     11,0xffffffff,7,0x9d2c5680,15,0xefc60000,18,1812433253> mt19937;
 

 
 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt11213b)
 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937_64)
 
+#include <boost/random/detail/enable_warnings.hpp>
+
 #endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP