ceph/src/boost/boost/random/faure.hpp

   1 /* boost random/faure.hpp header file
   2  *
   3  * Copyright Justinas Vygintas Daugmaudis 2010-2018
   4  * Distributed under the Boost Software License, Version 1.0. (See
   5  * accompanying file LICENSE_1_0.txt or copy at
   6  * http://www.boost.org/LICENSE_1_0.txt)
   7  */
   8
   9 #ifndef BOOST_RANDOM_FAURE_HPP
  10 #define BOOST_RANDOM_FAURE_HPP
  11
  12 #include <boost/random/detail/qrng_base.hpp>
  13
  14 #include <cmath>
  15 #include <vector>
  16 #include <algorithm>
  17
  18 #include <boost/assert.hpp>
  19
  20 namespace boost {
  21 namespace random {
  22
  23 /** @cond */
  24 namespace detail {
  25
  26 namespace qrng_tables {
  27
  28 // There is no particular reason why 187 first primes were chosen
  29 // to be put into this table. The only reason was, perhaps, that
  30 // the number of dimensions for Faure generator would be around
  31 // the same order of magnitude as the number of dimensions supported
  32 // by the Sobol qrng.
  33 struct primes
  34 {
  35   typedef unsigned short value_type;
  36
  37   BOOST_STATIC_CONSTANT(int, number_of_primes = 187);
  38
  39   // A function that returns lower bound prime for a given n
  40   static value_type lower_bound(std::size_t n)
  41   {
  42     static const value_type prim_a[number_of_primes] = {
  43       2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
  44       59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
  45       127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181,
  46       191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
  47       257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
  48       331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
  49       401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
  50       467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557,
  51       563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619,
  52       631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
  53       709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
  54       797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
  55       877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953,
  56       967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031,
  57       1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093,
  58       1097, 1103, 1109, 1117 };
  59
  60     qrng_detail::dimension_assert("Faure", n, prim_a[number_of_primes - 1]);
  61
  62     return *std::lower_bound(prim_a, prim_a + number_of_primes, n);
  63   }
  64 };
  65
  66 } // namespace qrng_tables
  67 } // namespace detail
  68
  69 namespace qrng_detail {
  70 namespace fr {
  71
  72 // Returns the integer part of the logarithm base Base of arg.
  73 // In erroneous situations, e.g., integer_log(base, 0) the function
  74 // returns 0 and does not report the error. This is the intended
  75 // behavior.
  76 template <typename T>
  77 inline T integer_log(T base, T arg)
  78 {
  79   T ilog = T();
  80   while (base <= arg)
  81   {
  82     arg /= base; ++ilog;
  83   }
  84   return ilog;
  85 }
  86
  87 // Perform exponentiation by squaring (potential for code reuse in multiprecision::powm)
  88 template <typename T>
  89 inline T integer_pow(T base, T e)
  90 {
  91   T result = static_cast<T>(1);
  92   while (e)
  93   {
  94     if (e & static_cast<T>(1))
  95       result *= base;
  96     e >>= 1;
  97     base *= base;
  98   }
  99   return result;
 100 }
 101
 102 } // namespace fr
 103
 104 // Computes a table of binomial coefficients modulo qs.
 105 template<typename RealType, typename SeqSizeT, typename PrimeTable>
 106 struct binomial_coefficients
 107 {
 108   typedef RealType value_type;
 109   typedef SeqSizeT size_type;
 110
 111   // Binomial values modulo qs_base will never be bigger than qs_base.
 112   // We can choose an appropriate integer type to hold modulo values and
 113   // shave off memory footprint.
 114   typedef typename PrimeTable::value_type packed_uint_t;
 115
 116   // default copy c-tor is fine
 117
 118   explicit binomial_coefficients(std::size_t dimension)
 119   {
 120     resize(dimension);
 121   }
 122
 123   void resize(std::size_t dimension)
 124   {
 125     qs_base = PrimeTable::lower_bound(dimension);
 126
 127     // Throw away previously computed coefficients.
 128     // This will trigger recomputation on next update
 129     coeff.clear();
 130   }
 131
 132   template <typename Iterator>
 133   void update(size_type seq, Iterator first, Iterator last)
 134   {
 135     if (first != last)
 136     {
 137       const size_type ilog = fr::integer_log(static_cast<size_type>(qs_base), seq);
 138       const size_type hisum = ilog + 1;
 139       if (coeff.size() != size_hint(hisum)) {
 140         ytemp.resize(static_cast<std::size_t>(hisum)); // cast safe because log is small
 141         compute_coefficients(hisum);
 142         qs_pow = fr::integer_pow(static_cast<size_type>(qs_base), ilog);
 143       }
 144
 145       *first = compute_recip(seq, ytemp.rbegin());
 146
 147       // Find other components using the Faure method.
 148       ++first;
 149       for ( ; first != last; ++first)
 150       {
 151         *first = RealType();
 152         RealType r = static_cast<RealType>(1);
 153
 154         for (size_type i = 0; i != hisum; ++i)
 155         {
 156           RealType ztemp = ytemp[static_cast<std::size_t>(i)] * upper_element(i, i, hisum);
 157           for (size_type j = i + 1; j != hisum; ++j)
 158             ztemp += ytemp[static_cast<std::size_t>(j)] * upper_element(i, j, hisum);
 159
 160           // Sum ( J <= I <= HISUM ) ( old ytemp(i) * binom(i,j) ) mod QS.
 161           ytemp[static_cast<std::size_t>(i)] = std::fmod(ztemp, static_cast<RealType>(qs_base));
 162           r *= static_cast<RealType>(qs_base);
 163           *first += ytemp[static_cast<std::size_t>(i)] / r;
 164         }
 165       }
 166     }
 167   }
 168
 169 private:
 170   inline static size_type size_hint(size_type n)
 171   {
 172     return n * (n + 1) / 2;
 173   }
 174
 175   packed_uint_t& upper_element(size_type i, size_type j, size_type dim)
 176   {
 177     BOOST_ASSERT( i < dim );
 178     BOOST_ASSERT( j < dim );
 179     BOOST_ASSERT( i <= j );
 180     return coeff[static_cast<std::size_t>((i * (2 * dim - i + 1)) / 2 + j - i)];
 181   }
 182
 183   template<typename Iterator>
 184   RealType compute_recip(size_type seq, Iterator out) const
 185   {
 186     // Here we do
 187     //   Sum ( 0 <= J <= HISUM ) YTEMP(J) * QS**J
 188     //   Sum ( 0 <= J <= HISUM ) YTEMP(J) / QS**(J+1)
 189     // in one go
 190     RealType r = RealType();
 191     size_type m, k = qs_pow;
 192     for( ; k != 0; ++out, seq = m, k /= qs_base )
 193     {
 194       m  = seq % k;
 195       RealType v  = static_cast<RealType>((seq - m) / k); // RealType <- size type
 196       r += v;
 197       r /= static_cast<RealType>(qs_base);
 198       *out = v; // saves double dereference
 199     }
 200     return r;
 201   }
 202
 203   void compute_coefficients(const size_type n)
 204   {
 205     // Resize and initialize to zero
 206     coeff.resize(static_cast<std::size_t>(size_hint(n)));
 207     std::fill(coeff.begin(), coeff.end(), packed_uint_t());
 208
 209     // The first row and the diagonal is assigned to 1
 210     upper_element(0, 0, n) = 1;
 211     for (size_type i = 1; i < n; ++i)
 212     {
 213       upper_element(0, i, n) = 1;
 214       upper_element(i, i, n) = 1;
 215     }
 216
 217     // Computes binomial coefficients MOD qs_base
 218     for (size_type i = 1; i < n; ++i)
 219     {
 220       for (size_type j = i + 1; j < n; ++j)
 221       {
 222         upper_element(i, j, n) = ( upper_element(i, j-1, n) +
 223                                    upper_element(i-1, j-1, n) ) % qs_base;
 224       }
 225     }
 226   }
 227
 228 private:
 229   packed_uint_t qs_base;
 230
 231   // here we cache precomputed data; note that binomial coefficients have
 232   // to be recomputed iff the integer part of the logarithm of seq changes,
 233   // which happens relatively rarely.
 234   std::vector<packed_uint_t> coeff; // packed upper (!) triangular matrix
 235   std::vector<RealType> ytemp;
 236   size_type qs_pow;
 237 };
 238
 239 } // namespace qrng_detail
 240
 241 typedef detail::qrng_tables::primes default_faure_prime_table;
 242
 243 /** @endcond */
 244
 245 //!Instantiations of class template faure_engine model a \quasi_random_number_generator.
 246 //!The faure_engine uses the algorithm described in
 247 //! \blockquote
 248 //!Henri Faure,
 249 //!Discrepance de suites associees a un systeme de numeration (en dimension s),
 250 //!Acta Arithmetica,
 251 //!Volume 41, 1982, pages 337-351.
 252 //! \endblockquote
 253 //
 254 //! \blockquote
 255 //!Bennett Fox,
 256 //!Algorithm 647:
 257 //!Implementation and Relative Efficiency of Quasirandom
 258 //!Sequence Generators,
 259 //!ACM Transactions on Mathematical Software,
 260 //!Volume 12, Number 4, December 1986, pages 362-376.
 261 //! \endblockquote
 262 //!
 263 //!In the following documentation @c X denotes the concrete class of the template
 264 //!faure_engine returning objects of type @c RealType, u and v are the values of @c X.
 265 //!
 266 //!Some member functions may throw exceptions of type @c std::bad_alloc.
 267 template<typename RealType, typename SeqSizeT, typename PrimeTable = default_faure_prime_table>
 268 class faure_engine
 269   : public qrng_detail::qrng_base<
 270       faure_engine<RealType, SeqSizeT, PrimeTable>
 271     , qrng_detail::binomial_coefficients<RealType, SeqSizeT, PrimeTable>
 272     , SeqSizeT
 273     >
 274 {
 275   typedef faure_engine<RealType, SeqSizeT, PrimeTable> self_t;
 276
 277   typedef qrng_detail::binomial_coefficients<RealType, SeqSizeT, PrimeTable> lattice_t;
 278   typedef qrng_detail::qrng_base<self_t, lattice_t, SeqSizeT> base_t;
 279
 280   friend class qrng_detail::qrng_base<self_t, lattice_t, SeqSizeT>;
 281
 282 public:
 283   typedef RealType result_type;
 284
 285   /** @copydoc boost::random::niederreiter_base2_engine::min() */
 286   static BOOST_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION ()
 287   { return static_cast<result_type>(0); }
 288
 289   /** @copydoc boost::random::niederreiter_base2_engine::max() */
 290   static BOOST_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION ()
 291   { return static_cast<result_type>(1); }
 292
 293   //!Effects: Constructs the `s`-dimensional default Faure quasi-random number generator.
 294   //!
 295   //!Throws: bad_alloc, invalid_argument.
 296   explicit faure_engine(std::size_t s)
 297     : base_t(s) // initialize the binomial table here
 298   {}
 299
 300   /** @copydetails boost::random::niederreiter_base2_engine::seed(UIntType)
 301    * Throws: bad_alloc.
 302    */
 303   void seed(SeqSizeT init = 0)
 304   {
 305     compute_seq(init);
 306     base_t::reset_seq(init);
 307   }
 308
 309 #ifdef BOOST_RANDOM_DOXYGEN
 310   //=========================Doxygen needs this!==============================
 311
 312   /** @copydoc boost::random::niederreiter_base2_engine::dimension() */
 313   std::size_t dimension() const { return base_t::dimension(); }
 314
 315   /** @copydoc boost::random::niederreiter_base2_engine::operator()() */
 316   result_type operator()()
 317   {
 318     return base_t::operator()();
 319   }
 320
 321   /** @copydoc boost::random::niederreiter_base2_engine::discard(boost::uintmax_t)
 322    * Throws: bad_alloc.
 323    */
 324   void discard(boost::uintmax_t z)
 325   {
 326     base_t::discard(z);
 327   }
 328
 329   /** Returns true if the two generators will produce identical sequences of outputs. */
 330   BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(faure_engine, x, y)
 331   { return static_cast<const base_t&>(x) == y; }
 332
 333   /** Returns true if the two generators will produce different sequences of outputs. */
 334   BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(faure_engine)
 335
 336   /** Writes the textual representation of the generator to a @c std::ostream. */
 337   BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, faure_engine, s)
 338   { return os << static_cast<const base_t&>(s); }
 339
 340   /** Reads the textual representation of the generator from a @c std::istream. */
 341   BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, faure_engine, s)
 342   { return is >> static_cast<base_t&>(s); }
 343
 344 #endif // BOOST_RANDOM_DOXYGEN
 345
 346 private:
 347 /** @cond hide_private_members */
 348   void compute_seq(SeqSizeT seq)
 349   {
 350     qrng_detail::check_seed_sign(seq);
 351     this->lattice.update(seq, this->state_begin(), this->state_end());
 352   }
 353 /** @endcond */
 354 };
 355
 356 /**
 357  * @attention This specialization of \faure_engine supports up to 1117 dimensions.
 358  *
 359  * However, it is possible to provide your own prime table to \faure_engine should the default one be insufficient.
 360  */
 361 typedef faure_engine<double, boost::uint_least64_t, default_faure_prime_table> faure;
 362
 363 } // namespace random
 364
 365 } // namespace boost
 366
 367 #endif // BOOST_RANDOM_FAURE_HPP