EIC code displayed by LXR master/​include/​boost/​random/​niederreiter_base2.hpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 09:51:12

 /* boost random/nierderreiter_base2.hpp header file
  *
  * Copyright Justinas Vygintas Daugmaudis 2010-2018
  * Distributed under the Boost Software License, Version 1.0. (See
  * accompanying file LICENSE_1_0.txt or copy at
  * http://www.boost.org/LICENSE_1_0.txt)
  */
 
 #ifndef BOOST_RANDOM_NIEDERREITER_BASE2_HPP
 #define BOOST_RANDOM_NIEDERREITER_BASE2_HPP
 
 #include <boost/random/detail/niederreiter_base2_table.hpp>
 #include <boost/random/detail/gray_coded_qrng.hpp>
 
 #include <boost/dynamic_bitset.hpp>
 
 namespace boost {
 namespace random {
 
 /** @cond */
 namespace qrng_detail {
 namespace nb2 {
 
 // Return the base 2 logarithm for a given bitset v
 template <typename DynamicBitset>
 inline typename DynamicBitset::size_type
 bitset_log2(const DynamicBitset& v)
 {
   if (v.none())
     boost::throw_exception( std::invalid_argument("bitset_log2") );
 
   typename DynamicBitset::size_type hibit = v.size() - 1;
   while (!v.test(hibit))
     --hibit;
   return hibit;
 }
 
 
 // Multiply polynomials over Z_2.
 template <typename PolynomialT, typename DynamicBitset>
 inline void modulo2_multiply(PolynomialT P, DynamicBitset& v, DynamicBitset& pt)
 {
   pt.reset(); // pt == 0
   for (; P; P >>= 1, v <<= 1)
     if (P & 1) pt ^= v;
   pt.swap(v);
 }
 
 
 // Calculate the values of the constants V(J,R) as
 // described in BFN section 3.3.
 //
 // pb = polynomial defined in section 2.3 of BFN.
 template <typename DynamicBitset>
 inline void calculate_v(const DynamicBitset& pb,
   typename DynamicBitset::size_type kj,
   typename DynamicBitset::size_type pb_degree,
   DynamicBitset& v)
 {
   typedef typename DynamicBitset::size_type size_type;
 
   // Now choose values of V in accordance with
   // the conditions in section 3.3.
   size_type r = 0;
   for ( ; r != kj; ++r)
     v.reset(r);
 
   // Quoting from BFN: "Our program currently sets each K_q
   // equal to eq. This has the effect of setting all unrestricted
   // values of v to 1."
   for ( ; r < pb_degree; ++r)
     v.set(r);
 
   // Calculate the remaining V's using the recursion of section 2.3,
   // remembering that the B's have the opposite sign.
   for ( ; r != v.size(); ++r)
   {
     bool term = false;
     for (typename DynamicBitset::size_type k = 0; k < pb_degree; ++k)
     {
       term ^= pb.test(k) & v[r + k - pb_degree];
     }
     v[r] = term;
   }
 }
 
 } // namespace nb2
 
 template<typename UIntType, unsigned w, typename Nb2Table>
 struct niederreiter_base2_lattice
 {
   typedef UIntType value_type;
 
   BOOST_STATIC_ASSERT(w > 0u);
   BOOST_STATIC_CONSTANT(unsigned, bit_count = w);
 
 private:
   typedef std::vector<value_type> container_type;
 
 public:
   explicit niederreiter_base2_lattice(std::size_t dimension)
   {
     resize(dimension);
   }
 
   void resize(std::size_t dimension)
   {
     typedef boost::dynamic_bitset<> bitset_type;
 
     dimension_assert("Niederreiter base 2", dimension, Nb2Table::max_dimension);
 
     // Initialize the bit array
     container_type cj(bit_count * dimension);
 
     // Reserve temporary space for lattice computation
     bitset_type v, pb, tmp;
 
     // Compute Niedderreiter base 2 lattice
     for (std::size_t dim = 0; dim != dimension; ++dim)
     {
       const typename Nb2Table::value_type poly = Nb2Table::polynomial(dim);
       if (poly > (std::numeric_limits<value_type>::max)()) {
         boost::throw_exception( std::range_error("niederreiter_base2: polynomial value outside the given value type range") );
       }
 
       const unsigned degree = qrng_detail::msb(poly); // integer log2(poly)
       const unsigned space_required = degree * ((bit_count / degree) + 1); // ~ degree + bit_count
 
       v.resize(degree + bit_count - 1);
 
       // For each dimension, we need to calculate powers of an
       // appropriate irreducible polynomial, see Niederreiter
       // page 65, just below equation (19).
       // Copy the appropriate irreducible polynomial into PX,
       // and its degree into E. Set polynomial B = PX ** 0 = 1.
       // M is the degree of B. Subsequently B will hold higher
       // powers of PX.
       pb.resize(space_required); tmp.resize(space_required);
 
       typename bitset_type::size_type kj, pb_degree = 0;
       pb.reset(); // pb == 0
       pb.set(pb_degree); // set the proper bit for the pb_degree
 
       value_type j = high_bit_mask_t<bit_count - 1>::high_bit;
       do
       {
         // Now choose a value of Kj as defined in section 3.3.
         // We must have 0 <= Kj < E*J = M.
         // The limit condition on Kj does not seem to be very relevant
         // in this program.
         kj = pb_degree;
 
         // Now multiply B by PX so B becomes PX**J.
         // In section 2.3, the values of Bi are defined with a minus sign :
         // don't forget this if you use them later!
         nb2::modulo2_multiply(poly, pb, tmp);
         pb_degree += degree;
         if (pb_degree >= pb.size()) {
           // Note that it is quite possible for kj to become bigger than
           // the new computed value of pb_degree.
           pb_degree = nb2::bitset_log2(pb);
         }
 
         // If U = 0, we need to set B to the next power of PX
         // and recalculate V.
         nb2::calculate_v(pb, kj, pb_degree, v);
 
         // Niederreiter (page 56, after equation (7), defines two
         // variables Q and U.  We do not need Q explicitly, but we
         // do need U.
 
         // Advance Niederreiter's state variables.
         for (unsigned u = 0; j && u != degree; ++u, j >>= 1)
         {
           // Now C is obtained from V. Niederreiter
           // obtains A from V (page 65, near the bottom), and then gets
           // C from A (page 56, equation (7)).  However this can be done
           // in one step.  Here CI(J,R) corresponds to
           // Niederreiter's C(I,J,R), whose values we pack into array
           // CJ so that CJ(I,R) holds all the values of C(I,J,R) for J from 1 to NBITS.
           for (unsigned r = 0; r != bit_count; ++r) {
             value_type& num = cj[dimension * r + dim];
             // set the jth bit in num
             num = (num & ~j) | (-v[r + u] & j);
           }
         }
       } while (j != 0);
     }
 
     bits.swap(cj);
   }
 
   typename container_type::const_iterator iter_at(std::size_t n) const
   {
     BOOST_ASSERT(!(n > bits.size()));
     return bits.begin() + n;
   }
 
 private:
   container_type bits;
 };
 
 } // namespace qrng_detail
 
 typedef detail::qrng_tables::niederreiter_base2 default_niederreiter_base2_table;
 
 /** @endcond */
 
 //!Instantiations of class template niederreiter_base2_engine model a \quasi_random_number_generator.
 //!The niederreiter_base2_engine uses the algorithm described in
 //! \blockquote
 //!Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992).
 //! \endblockquote
 //!
 //!\attention niederreiter_base2_engine skips trivial zeroes at the start of the sequence. For example,
 //!the beginning of the 2-dimensional Niederreiter base 2 sequence in @c uniform_01 distribution will look
 //!like this:
 //!\code{.cpp}
 //!0.5, 0.5,
 //!0.75, 0.25,
 //!0.25, 0.75,
 //!0.375, 0.375,
 //!0.875, 0.875,
 //!...
 //!\endcode
 //!
 //!In the following documentation @c X denotes the concrete class of the template
 //!niederreiter_base2_engine returning objects of type @c UIntType, u and v are the values of @c X.
 //!
 //!Some member functions may throw exceptions of type std::range_error. This
 //!happens when the quasi-random domain is exhausted and the generator cannot produce
 //!any more values. The length of the low discrepancy sequence is given by
 //! \f$L=Dimension \times (2^{w} - 1)\f$.
 template<typename UIntType, unsigned w, typename Nb2Table = default_niederreiter_base2_table>
 class niederreiter_base2_engine
   : public qrng_detail::gray_coded_qrng<
       qrng_detail::niederreiter_base2_lattice<UIntType, w, Nb2Table>
     >
 {
   typedef qrng_detail::niederreiter_base2_lattice<UIntType, w, Nb2Table> lattice_t;
   typedef qrng_detail::gray_coded_qrng<lattice_t> base_t;
 
 public:
   //!Effects: Constructs the default `s`-dimensional Niederreiter base 2 quasi-random number generator.
   //!
   //!Throws: bad_alloc, invalid_argument, range_error.
   explicit niederreiter_base2_engine(std::size_t s)
     : base_t(s) // initialize lattice here
   {}
 
 #ifdef BOOST_RANDOM_DOXYGEN
   //=========================Doxygen needs this!==============================
   typedef UIntType result_type;
 
   //!Returns: Tight lower bound on the set of values returned by operator().
   //!
   //!Throws: nothing.
   static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
   { return (base_t::min)(); }
 
   //!Returns: Tight upper bound on the set of values returned by operator().
   //!
   //!Throws: nothing.
   static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
   { return (base_t::max)(); }
 
   //!Returns: The dimension of of the quasi-random domain.
   //!
   //!Throws: nothing.
   std::size_t dimension() const { return base_t::dimension(); }
 
   //!Effects: Resets the quasi-random number generator state to
   //!the one given by the default construction. Equivalent to u.seed(0).
   //!
   //!\brief Throws: nothing.
   void seed()
   {
     base_t::seed();
   }
 
   //!Effects: Effectively sets the quasi-random number generator state to the `init`-th
   //!vector in the `s`-dimensional quasi-random domain, where `s` == X::dimension().
   //!\code
   //!X u, v;
   //!for(int i = 0; i < N; ++i)
   //!    for( std::size_t j = 0; j < u.dimension(); ++j )
   //!        u();
   //!v.seed(N);
   //!assert(u() == v());
   //!\endcode
   //!
   //!\brief Throws: range_error.
   void seed(UIntType init)
   {
     base_t::seed(init);
   }
 
   //!Returns: Returns a successive element of an `s`-dimensional
   //!(s = X::dimension()) vector at each invocation. When all elements are
   //!exhausted, X::operator() begins anew with the starting element of a
   //!subsequent `s`-dimensional vector.
   //!
   //!Throws: range_error.
   result_type operator()()
   {
     return base_t::operator()();
   }
 
   //!Effects: Advances *this state as if `z` consecutive
   //!X::operator() invocations were executed.
   //!\code
   //!X u = v;
   //!for(int i = 0; i < N; ++i)
   //!    u();
   //!v.discard(N);
   //!assert(u() == v());
   //!\endcode
   //!
   //!Throws: range_error.
   void discard(boost::uintmax_t z)
   {
     base_t::discard(z);
   }
 
   //!Returns true if the two generators will produce identical sequences of outputs.
   BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(niederreiter_base2_engine, x, y)
   { return static_cast<const base_t&>(x) == y; }
 
   //!Returns true if the two generators will produce different sequences of outputs.
   BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(niederreiter_base2_engine)
 
   //!Writes the textual representation of the generator to a @c std::ostream.
   BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, niederreiter_base2_engine, s)
   { return os << static_cast<const base_t&>(s); }
 
   //!Reads the textual representation of the generator from a @c std::istream.
   BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, niederreiter_base2_engine, s)
   { return is >> static_cast<base_t&>(s); }
 
 #endif // BOOST_RANDOM_DOXYGEN
 };
 
 
 /**
  * @attention This specialization of \niederreiter_base2_engine supports up to 4720 dimensions.
  *
  * Binary irreducible polynomials (primes in the ring `GF(2)[X]`, evaluated at `X=2`) were generated
  * while condition `max(prime)` < 2<sup>16</sup> was satisfied.
  *
  * There are exactly 4720 such primes, which yields a Niederreiter base 2 table for 4720 dimensions.
  *
  * However, it is possible to provide your own table to \niederreiter_base2_engine should the default one be insufficient.
  */
 typedef niederreiter_base2_engine<boost::uint_least64_t, 64u, default_niederreiter_base2_table> niederreiter_base2;
 
 } // namespace random
 
 } // namespace boost
 
 #endif // BOOST_RANDOM_NIEDERREITER_BASE2_HPP

This page was automatically generated by the 2.3.7 LXR engine. The LXR team