boost/random/inversive_congruential.hpp

   1 /* boost random/inversive_congruential.hpp header file
   2  *
   3  * Copyright Jens Maurer 2000-2001
   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  * See http://www.boost.org for most recent version including documentation.
   9  *
  10  * $Id: inversive_congruential.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $
  11  *
  12  * Revision history
  13  *  2001-02-18  moved to individual header files
  14  */
  15
  16 #ifndef BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP
  17 #define BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP
  18
  19 #include <iostream>
  20 #include <cassert>
  21 #include <stdexcept>
  22 #include <boost/config.hpp>
  23 #include <boost/static_assert.hpp>
  24 #include <boost/random/detail/config.hpp>
  25 #include <boost/random/detail/const_mod.hpp>
  26
  27 namespace boost {
  28 namespace random {
  29
  30 // Eichenauer and Lehn 1986
  31 /**
  32  * Instantiations of class template @c inversive_congruential model a
  33  * \pseudo_random_number_generator. It uses the inversive congruential
  34  * algorithm (ICG) described in
  35  *
  36  *  @blockquote
  37  *  "Inversive pseudorandom number generators: concepts, results and links",
  38  *  Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation
  39  *  Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman
  40  *  (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps
  41  *  @endblockquote
  42  *
  43  * The output sequence is defined by x(n+1) = (a*inv(x(n)) - b) (mod p),
  44  * where x(0), a, b, and the prime number p are parameters of the generator.
  45  * The expression inv(k) denotes the multiplicative inverse of k in the
  46  * field of integer numbers modulo p, with inv(0) := 0.
  47  *
  48  * The template parameter IntType shall denote a signed integral type large
  49  * enough to hold p; a, b, and p are the parameters of the generators. The
  50  * template parameter val is the validation value checked by validation.
  51  *
  52  * @xmlnote
  53  * The implementation currently uses the Euclidian Algorithm to compute
  54  * the multiplicative inverse. Therefore, the inversive generators are about
  55  * 10-20 times slower than the others (see section"performance"). However,
  56  * the paper talks of only 3x slowdown, so the Euclidian Algorithm is probably
  57  * not optimal for calculating the multiplicative inverse.
  58  * @endxmlnote
  59  */
  60 template<class IntType, IntType a, IntType b, IntType p, IntType val>
  61 class inversive_congruential
  62 {
  63 public:
  64   typedef IntType result_type;
  65 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
  66   static const bool has_fixed_range = true;
  67   static const result_type min_value = (b == 0 ? 1 : 0);
  68   static const result_type max_value = p-1;
  69 #else
  70   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  71 #endif
  72   BOOST_STATIC_CONSTANT(result_type, multiplier = a);
  73   BOOST_STATIC_CONSTANT(result_type, increment = b);
  74   BOOST_STATIC_CONSTANT(result_type, modulus = p);
  75
  76   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return b == 0 ? 1 : 0; }
  77   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return p-1; }
  78
  79   /**
  80    * Constructs an inversive_congruential generator with
  81    * @c y0 as the initial state.
  82    */
  83   explicit inversive_congruential(IntType y0 = 1) : value(y0)
  84   {
  85     BOOST_STATIC_ASSERT(b >= 0);
  86     BOOST_STATIC_ASSERT(p > 1);
  87     BOOST_STATIC_ASSERT(a >= 1);
  88     if(b == 0)
  89       assert(y0 > 0);
  90   }
  91   template<class It> inversive_congruential(It& first, It last)
  92   { seed(first, last); }
  93
  94   /** Changes the current state to y0. */
  95   void seed(IntType y0 = 1) { value = y0; if(b == 0) assert(y0 > 0); }
  96   template<class It> void seed(It& first, It last)
  97   {
  98     if(first == last)
  99       throw std::invalid_argument("inversive_congruential::seed");
 100     value = *first++;
 101   }
 102   IntType operator()()
 103   {
 104     typedef const_mod<IntType, p> do_mod;
 105     value = do_mod::mult_add(a, do_mod::invert(value), b);
 106     return value;
 107   }
 108
 109   static bool validation(result_type x) { return val == x; }
 110
 111 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
 112
 113 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
 114   template<class CharT, class Traits>
 115   friend std::basic_ostream<CharT,Traits>&
 116   operator<<(std::basic_ostream<CharT,Traits>& os, inversive_congruential x)
 117   { os << x.value; return os; }
 118
 119   template<class CharT, class Traits>
 120   friend std::basic_istream<CharT,Traits>&
 121   operator>>(std::basic_istream<CharT,Traits>& is, inversive_congruential& x)
 122   { is >> x.value; return is; }
 123 #endif
 124
 125   friend bool operator==(inversive_congruential x, inversive_congruential y)
 126   { return x.value == y.value; }
 127   friend bool operator!=(inversive_congruential x, inversive_congruential y)
 128   { return !(x == y); }
 129 #else
 130   // Use a member function; Streamable concept not supported.
 131   bool operator==(inversive_congruential rhs) const
 132   { return value == rhs.value; }
 133   bool operator!=(inversive_congruential rhs) const
 134   { return !(*this == rhs); }
 135 #endif
 136 private:
 137   IntType value;
 138 };
 139
 140 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
 141 //  A definition is required even for integral static constants
 142 template<class IntType, IntType a, IntType b, IntType p, IntType val>
 143 const bool inversive_congruential<IntType, a, b, p, val>::has_fixed_range;
 144 template<class IntType, IntType a, IntType b, IntType p, IntType val>
 145 const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::min_value;
 146 template<class IntType, IntType a, IntType b, IntType p, IntType val>
 147 const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::max_value;
 148 template<class IntType, IntType a, IntType b, IntType p, IntType val>
 149 const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::multiplier;
 150 template<class IntType, IntType a, IntType b, IntType p, IntType val>
 151 const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::increment;
 152 template<class IntType, IntType a, IntType b, IntType p, IntType val>
 153 const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::modulus;
 154 #endif
 155
 156 } // namespace random
 157
 158 /**
 159  * The specialization hellekalek1995 was suggested in
 160  *
 161  *  @blockquote
 162  *  "Inversive pseudorandom number generators: concepts, results and links",
 163  *  Peter Hellekalek, In: "Proceedings of the 1995 Winter Simulation
 164  *  Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and D. Goldsman
 165  *  (editors), 1995, pp. 255-262. ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps
 166  *  @endblockquote
 167  */
 168 typedef random::inversive_congruential<int32_t, 9102, 2147483647-36884165,
 169   2147483647, 0> hellekalek1995;
 170
 171 } // namespace boost
 172
 173 #endif // BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP