boost/random/faure.hpp

0001 /* boost random/faure.hpp header file
0002  *
0003  * Copyright Justinas Vygintas Daugmaudis 2010-2018
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  */
0008
0009 #ifndef BOOST_RANDOM_FAURE_HPP
0010 #define BOOST_RANDOM_FAURE_HPP
0011
0012 #include <boost/random/detail/qrng_base.hpp>
0013
0014 #include <cmath>
0015 #include <vector>
0016 #include <algorithm>
0017
0018 #include <boost/assert.hpp>
0019
0020 namespace boost {
0021 namespace random {
0022
0023 /** @cond */
0024 namespace detail {
0025
0026 namespace qrng_tables {
0027
0028 // There is no particular reason why 187 first primes were chosen
0029 // to be put into this table. The only reason was, perhaps, that
0030 // the number of dimensions for Faure generator would be around
0031 // the same order of magnitude as the number of dimensions supported
0032 // by the Sobol qrng.
0033 struct primes
0034 {
0035   typedef unsigned short value_type;
0036
0037   BOOST_STATIC_CONSTANT(int, number_of_primes = 187);
0038
0039   // A function that returns lower bound prime for a given n
0040   static value_type lower_bound(std::size_t n)
0041   {
0042     static const value_type prim_a[number_of_primes] = {
0043       2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
0044       59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
0045       127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
0046       191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
0047       257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
0048       331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
0049       401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
0050       467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557,
0051       563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619,
0052       631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
0053       709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
0054       797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
0055       877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953,
0056       967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031,
0057       1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093,
0058       1097, 1103, 1109, 1117 };
0059
0060     qrng_detail::dimension_assert("Faure", n, prim_a[number_of_primes - 1]);
0061
0062     return *std::lower_bound(prim_a, prim_a + number_of_primes, n);
0063   }
0064 };
0065
0066 } // namespace qrng_tables
0067 } // namespace detail
0068
0069 namespace qrng_detail {
0070 namespace fr {
0071
0072 // Returns the integer part of the logarithm base Base of arg.
0073 // In erroneous situations, e.g., integer_log(base, 0) the function
0074 // returns 0 and does not report the error. This is the intended
0075 // behavior.
0076 template <typename T>
0077 inline T integer_log(T base, T arg)
0078 {
0079   T ilog = T();
0080   while (base <= arg)
0081   {
0082     arg /= base; ++ilog;
0083   }
0084   return ilog;
0085 }
0086
0087 // Perform exponentiation by squaring (potential for code reuse in multiprecision::powm)
0088 template <typename T>
0089 inline T integer_pow(T base, T e)
0090 {
0091   T result = static_cast<T>(1);
0092   while (e)
0093   {
0094     if (e & static_cast<T>(1))
0095       result *= base;
0096     e >>= 1;
0097     base *= base;
0098   }
0099   return result;
0100 }
0101
0102 } // namespace fr
0103
0104 // Computes a table of binomial coefficients modulo qs.
0105 template<typename RealType, typename SeqSizeT, typename PrimeTable>
0106 struct binomial_coefficients
0107 {
0108   typedef RealType value_type;
0109   typedef SeqSizeT size_type;
0110
0111   // Binomial values modulo qs_base will never be bigger than qs_base.
0112   // We can choose an appropriate integer type to hold modulo values and
0113   // shave off memory footprint.
0114   typedef typename PrimeTable::value_type packed_uint_t;
0115
0116   // default copy c-tor is fine
0117
0118   explicit binomial_coefficients(std::size_t dimension)
0119   {
0120     resize(dimension);
0121   }
0122
0123   void resize(std::size_t dimension)
0124   {
0125     qs_base = PrimeTable::lower_bound(dimension);
0126
0127     // Throw away previously computed coefficients.
0128     // This will trigger recomputation on next update
0129     coeff.clear();
0130   }
0131
0132   template <typename Iterator>
0133   void update(size_type seq, Iterator first, Iterator last)
0134   {
0135     if (first != last)
0136     {
0137       const size_type ilog = fr::integer_log(static_cast<size_type>(qs_base), seq);
0138       const size_type hisum = ilog + 1;
0139       if (coeff.size() != size_hint(hisum)) {
0140         ytemp.resize(static_cast<std::size_t>(hisum)); // cast safe because log is small
0141         compute_coefficients(hisum);
0142         qs_pow = fr::integer_pow(static_cast<size_type>(qs_base), ilog);
0143       }
0144
0145       *first = compute_recip(seq, ytemp.rbegin());
0146
0147       // Find other components using the Faure method.
0148       ++first;
0149       for ( ; first != last; ++first)
0150       {
0151         *first = RealType();
0152         RealType r = static_cast<RealType>(1);
0153
0154         for (size_type i = 0; i != hisum; ++i)
0155         {
0156           RealType ztemp = ytemp[static_cast<std::size_t>(i)] * upper_element(i, i, hisum);
0157           for (size_type j = i + 1; j != hisum; ++j)
0158             ztemp += ytemp[static_cast<std::size_t>(j)] * upper_element(i, j, hisum);
0159
0160           // Sum ( J <= I <= HISUM ) ( old ytemp(i) * binom(i,j) ) mod QS.
0161           ytemp[static_cast<std::size_t>(i)] = std::fmod(ztemp, static_cast<RealType>(qs_base));
0162           r *= static_cast<RealType>(qs_base);
0163           *first += ytemp[static_cast<std::size_t>(i)] / r;
0164         }
0165       }
0166     }
0167   }
0168
0169 private:
0170   inline static size_type size_hint(size_type n)
0171   {
0172     return n * (n + 1) / 2;
0173   }
0174
0175   packed_uint_t& upper_element(size_type i, size_type j, size_type dim)
0176   {
0177     BOOST_ASSERT( i < dim );
0178     BOOST_ASSERT( j < dim );
0179     BOOST_ASSERT( i <= j );
0180     return coeff[static_cast<std::size_t>((i * (2 * dim - i + 1)) / 2 + j - i)];
0181   }
0182
0183   template<typename Iterator>
0184   RealType compute_recip(size_type seq, Iterator out) const
0185   {
0186     // Here we do
0187     //   Sum ( 0 <= J <= HISUM ) YTEMP(J) * QS**J
0188     //   Sum ( 0 <= J <= HISUM ) YTEMP(J) / QS**(J+1)
0189     // in one go
0190     RealType r = RealType();
0191     size_type m, k = qs_pow;
0192     for( ; k != 0; ++out, seq = m, k /= qs_base )
0193     {
0194       m  = seq % k;
0195       RealType v  = static_cast<RealType>((seq - m) / k); // RealType <- size type
0196       r += v;
0197       r /= static_cast<RealType>(qs_base);
0198       *out = v; // saves double dereference
0199     }
0200     return r;
0201   }
0202
0203   void compute_coefficients(const size_type n)
0204   {
0205     // Resize and initialize to zero
0206     coeff.resize(static_cast<std::size_t>(size_hint(n)));
0207     std::fill(coeff.begin(), coeff.end(), packed_uint_t());
0208
0209     // The first row and the diagonal is assigned to 1
0210     upper_element(0, 0, n) = 1;
0211     for (size_type i = 1; i < n; ++i)
0212     {
0213       upper_element(0, i, n) = 1;
0214       upper_element(i, i, n) = 1;
0215     }
0216
0217     // Computes binomial coefficients MOD qs_base
0218     for (size_type i = 1; i < n; ++i)
0219     {
0220       for (size_type j = i + 1; j < n; ++j)
0221       {
0222         upper_element(i, j, n) = ( upper_element(i, j-1, n) +
0223                                    upper_element(i-1, j-1, n) ) % qs_base;
0224       }
0225     }
0226   }
0227
0228 private:
0229   packed_uint_t qs_base;
0230
0231   // here we cache precomputed data; note that binomial coefficients have
0232   // to be recomputed iff the integer part of the logarithm of seq changes,
0233   // which happens relatively rarely.
0234   std::vector<packed_uint_t> coeff; // packed upper (!) triangular matrix
0235   std::vector<RealType> ytemp;
0236   size_type qs_pow;
0237 };
0238
0239 } // namespace qrng_detail
0240
0241 typedef detail::qrng_tables::primes default_faure_prime_table;
0242
0243 /** @endcond */
0244
0245 //!Instantiations of class template faure_engine model a \quasi_random_number_generator.
0246 //!The faure_engine uses the algorithm described in
0247 //! \blockquote
0248 //!Henri Faure,
0249 //!Discrepance de suites associees a un systeme de numeration (en dimension s),
0250 //!Acta Arithmetica,
0251 //!Volume 41, 1982, pages 337-351.
0252 //! \endblockquote
0253 //
0254 //! \blockquote
0255 //!Bennett Fox,
0256 //!Algorithm 647:
0257 //!Implementation and Relative Efficiency of Quasirandom
0258 //!Sequence Generators,
0259 //!ACM Transactions on Mathematical Software,
0260 //!Volume 12, Number 4, December 1986, pages 362-376.
0261 //! \endblockquote
0262 //!
0263 //!In the following documentation @c X denotes the concrete class of the template
0264 //!faure_engine returning objects of type @c RealType, u and v are the values of @c X.
0265 //!
0266 //!Some member functions may throw exceptions of type @c std::bad_alloc.
0267 template<typename RealType, typename SeqSizeT, typename PrimeTable = default_faure_prime_table>
0268 class faure_engine
0269   : public qrng_detail::qrng_base<
0270       faure_engine<RealType, SeqSizeT, PrimeTable>
0271     , qrng_detail::binomial_coefficients<RealType, SeqSizeT, PrimeTable>
0272     , SeqSizeT
0273     >
0274 {
0275   typedef faure_engine<RealType, SeqSizeT, PrimeTable> self_t;
0276
0277   typedef qrng_detail::binomial_coefficients<RealType, SeqSizeT, PrimeTable> lattice_t;
0278   typedef qrng_detail::qrng_base<self_t, lattice_t, SeqSizeT> base_t;
0279
0280   friend class qrng_detail::qrng_base<self_t, lattice_t, SeqSizeT>;
0281
0282 public:
0283   typedef RealType result_type;
0284
0285   /** @copydoc boost::random::niederreiter_base2_engine::min() */
0286   static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
0287   { return static_cast<result_type>(0); }
0288
0289   /** @copydoc boost::random::niederreiter_base2_engine::max() */
0290   static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
0291   { return static_cast<result_type>(1); }
0292
0293   //!Effects: Constructs the `s`-dimensional default Faure quasi-random number generator.
0294   //!
0295   //!Throws: bad_alloc, invalid_argument.
0296   explicit faure_engine(std::size_t s)
0297     : base_t(s) // initialize the binomial table here
0298   {}
0299
0300   /** @copydetails boost::random::niederreiter_base2_engine::seed(UIntType)
0301    * Throws: bad_alloc.
0302    */
0303   void seed(SeqSizeT init = 0)
0304   {
0305     compute_seq(init);
0306     base_t::reset_seq(init);
0307   }
0308
0309 #ifdef BOOST_RANDOM_DOXYGEN
0310   //=========================Doxygen needs this!==============================
0311
0312   /** @copydoc boost::random::niederreiter_base2_engine::dimension() */
0313   std::size_t dimension() const { return base_t::dimension(); }
0314
0315   /** @copydoc boost::random::niederreiter_base2_engine::operator()() */
0316   result_type operator()()
0317   {
0318     return base_t::operator()();
0319   }
0320
0321   /** @copydoc boost::random::niederreiter_base2_engine::discard(boost::uintmax_t)
0322    * Throws: bad_alloc.
0323    */
0324   void discard(boost::uintmax_t z)
0325   {
0326     base_t::discard(z);
0327   }
0328
0329   /** Returns true if the two generators will produce identical sequences of outputs. */
0330   BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(faure_engine, x, y)
0331   { return static_cast<const base_t&>(x) == y; }
0332
0333   /** Returns true if the two generators will produce different sequences of outputs. */
0334   BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(faure_engine)
0335
0336   /** Writes the textual representation of the generator to a @c std::ostream. */
0337   BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, faure_engine, s)
0338   { return os << static_cast<const base_t&>(s); }
0339
0340   /** Reads the textual representation of the generator from a @c std::istream. */
0341   BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, faure_engine, s)
0342   { return is >> static_cast<base_t&>(s); }
0343
0344 #endif // BOOST_RANDOM_DOXYGEN
0345
0346 private:
0347 /** @cond hide_private_members */
0348   void compute_seq(SeqSizeT seq)
0349   {
0350     qrng_detail::check_seed_sign(seq);
0351     this->lattice.update(seq, this->state_begin(), this->state_end());
0352   }
0353 /** @endcond */
0354 };
0355
0356 /**
0357  * @attention This specialization of \faure_engine supports up to 1117 dimensions.
0358  *
0359  * However, it is possible to provide your own prime table to \faure_engine should the default one be insufficient.
0360  */
0361 typedef faure_engine<double, boost::uint_least64_t, default_faure_prime_table> faure;
0362
0363 } // namespace random
0364
0365 } // namespace boost
0366
0367 #endif // BOOST_RANDOM_FAURE_HPP