1 /* boost random/lagged_fibonacci.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: lagged_fibonacci.hpp,v 1.28 2005/05/21 15:57:00 dgregor Exp $

    11  *

    12  * Revision history

    13  *  2001-02-18  moved to individual header files

    14  */

    15 

    16 #ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP

    17 #define BOOST_RANDOM_LAGGED_FIBONACCI_HPP

    18 

    19 #include <cmath>

    20 #include <iostream>

    21 #include <algorithm>     // std::max

    22 #include <iterator>

    23 #include <cmath>         // std::pow

    24 #include <boost/config.hpp>

    25 #include <boost/limits.hpp>

    26 #include <boost/cstdint.hpp>

    27 #include <boost/detail/workaround.hpp>

    28 #include <boost/random/linear_congruential.hpp>

    29 #include <boost/random/uniform_01.hpp>

    30 #include <boost/random/detail/pass_through_engine.hpp>

    31 

    32 namespace boost {

    33 namespace random {

    34 

    35 #if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300

    36 #  define BOOST_RANDOM_EXTRACT_LF

    37 #endif

    38 

    39 #if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)

    40 #  define BOOST_RANDOM_EXTRACT_LF

    41 #endif

    42 

    43 #  ifdef BOOST_RANDOM_EXTRACT_LF

    44 namespace detail

    45 {

    46   template<class IStream, class F, class RealType>

    47   IStream&

    48   extract_lagged_fibonacci_01(

    49       IStream& is

    50       , F const& f

    51       , unsigned int& i

    52       , RealType* x

    53       , RealType modulus)

    54   {

    55         is >> i >> std::ws;

    56         for(unsigned int i = 0; i < f.long_lag; ++i)

    57         {

    58             RealType value;

    59             is >> value >> std::ws;

    60             x[i] = value / modulus;

    61         }

    62         return is;

    63   }

    64 

    65   template<class IStream, class F, class UIntType>

    66   IStream&

    67   extract_lagged_fibonacci(

    68       IStream& is

    69       , F const& f

    70       , unsigned int& i

    71       , UIntType* x)

    72   {

    73       is >> i >> std::ws;

    74       for(unsigned int i = 0; i < f.long_lag; ++i)

    75           is >> x[i] >> std::ws;

    76       return is;

    77   }

    78 }

    79 #  endif

    80 

    81 template<class UIntType, int w, unsigned int p, unsigned int q,

    82          UIntType val = 0>

    83 class lagged_fibonacci

    84 {

    85 public:

    86   typedef UIntType result_type;

    87   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);

    88   BOOST_STATIC_CONSTANT(int, word_size = w);

    89   BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);

    90   BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);

    91 

    92   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }

    93   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }

    94 

    95   lagged_fibonacci() { init_wordmask(); seed(); }

    96   explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }

    97   template<class It> lagged_fibonacci(It& first, It last)

    98   { init_wordmask(); seed(first, last); }

    99   // compiler-generated copy ctor and assignment operator are fine

   100 

   101 private:

   102   void init_wordmask()

   103   {

   104     wordmask = 0;

   105     for(int i = 0; i < w; ++i)

   106       wordmask |= (1u << i);

   107   }

   108 

   109 public:

   110   void seed(uint32_t value = 331u)

   111   {

   112     minstd_rand0 gen(value);

   113     for(unsigned int j = 0; j < long_lag; ++j)

   114       x[j] = gen() & wordmask;

   115     i = long_lag;

   116   }

   117 

   118   template<class It>

   119   void seed(It& first, It last)

   120   {

   121     // word size could be smaller than the seed values

   122     unsigned int j;

   123     for(j = 0; j < long_lag && first != last; ++j, ++first)

   124       x[j] = *first & wordmask;

   125     i = long_lag;

   126     if(first == last && j < long_lag)

   127       throw std::invalid_argument("lagged_fibonacci::seed");

   128   }

   129 

   130   result_type operator()()

   131   {

   132     if(i >= long_lag)

   133       fill();

   134     return x[i++];

   135   }

   136 

   137   static bool validation(result_type x)

   138   {

   139     return x == val;

   140   }

   141   

   142 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE

   143 

   144 #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS

   145   template<class CharT, class Traits>

   146   friend std::basic_ostream<CharT,Traits>&

   147   operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)

   148   {

   149     os << f.i << " ";

   150     for(unsigned int i = 0; i < f.long_lag; ++i)

   151       os << f.x[i] << " ";

   152     return os;

   153   }

   154 

   155   template<class CharT, class Traits>

   156   friend std::basic_istream<CharT, Traits>&

   157   operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)

   158   {

   159 # ifdef BOOST_RANDOM_EXTRACT_LF

   160       return detail::extract_lagged_fibonacci(is, f, f.i, f.x);

   161 # else

   162       is >> f.i >> std::ws;

   163       for(unsigned int i = 0; i < f.long_lag; ++i)

   164           is >> f.x[i] >> std::ws;

   165       return is;

   166 # endif 

   167   }

   168 #endif

   169 

   170   friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)

   171   { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }

   172   friend bool operator!=(const lagged_fibonacci& x,

   173                          const lagged_fibonacci& y)

   174   { return !(x == y); }

   175 #else

   176   // Use a member function; Streamable concept not supported.

   177   bool operator==(const lagged_fibonacci& rhs) const

   178   { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }

   179   bool operator!=(const lagged_fibonacci& rhs) const

   180   { return !(*this == rhs); }

   181 #endif

   182 

   183 private:

   184   void fill();

   185   UIntType wordmask;

   186   unsigned int i;

   187   UIntType x[long_lag];

   188 };

   189 

   190 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION

   191 //  A definition is required even for integral static constants

   192 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>

   193 const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;

   194 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>

   195 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;

   196 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>

   197 const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;

   198 #endif

   199 

   200 template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>

   201 void lagged_fibonacci<UIntType, w, p, q, val>::fill()

   202 {

   203   // two loops to avoid costly modulo operations

   204   {  // extra scope for MSVC brokenness w.r.t. for scope

   205   for(unsigned int j = 0; j < short_lag; ++j)

   206     x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;

   207   }

   208   for(unsigned int j = short_lag; j < long_lag; ++j)

   209     x[j] = (x[j] + x[j-short_lag]) & wordmask;

   210   i = 0;

   211 }

   212 

   213 

   214 

   215 // lagged Fibonacci generator for the range [0..1)

   216 // contributed by Matthias Troyer

   217 // for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958

   218 

   219 template<class T, unsigned int p, unsigned int q>

   220 struct fibonacci_validation

   221 {

   222   BOOST_STATIC_CONSTANT(bool, is_specialized = false);

   223   static T value() { return 0; }

   224   static T tolerance() { return 0; }

   225 };

   226 

   227 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION

   228 //  A definition is required even for integral static constants

   229 template<class T, unsigned int p, unsigned int q>

   230 const bool fibonacci_validation<T, p, q>::is_specialized;

   231 #endif

   232 

   233 #define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \

   234 template<> \

   235 struct fibonacci_validation<T, P, Q>  \

   236 {                                     \

   237   BOOST_STATIC_CONSTANT(bool, is_specialized = true);     \

   238   static T value() { return V; }      \

   239   static T tolerance()                \

   240 { return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \

   241 };

   242 // (The extra static_cast<T> in the std::max call above is actually

   243 // unnecessary except for HP aCC 1.30, which claims that

   244 // numeric_limits<double>::epsilon() doesn't actually return a double.)

   245 

   246 BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)

   247 BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)

   248 BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)

   249 BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)

   250 BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)

   251 BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)

   252 BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)

   253 BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)

   254 BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)

   255 

   256 #undef BOOST_RANDOM_FIBONACCI_VAL

   257 

   258 template<class RealType, int w, unsigned int p, unsigned int q>

   259 class lagged_fibonacci_01

   260 {

   261 public:

   262   typedef RealType result_type;

   263   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);

   264   BOOST_STATIC_CONSTANT(int, word_size = w);

   265   BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);

   266   BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);

   267 

   268   lagged_fibonacci_01() { init_modulus(); seed(); }

   269   explicit lagged_fibonacci_01(uint32_t value) { init_modulus(); seed(value); }

   270   template<class Generator>

   271   explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); }

   272   template<class It> lagged_fibonacci_01(It& first, It last)

   273   { init_modulus(); seed(first, last); }

   274   // compiler-generated copy ctor and assignment operator are fine

   275 

   276 private:

   277   void init_modulus()

   278   {

   279 #ifndef BOOST_NO_STDC_NAMESPACE

   280     // allow for Koenig lookup

   281     using std::pow;

   282 #endif

   283     _modulus = pow(RealType(2), word_size);

   284   }

   285 

   286 public:

   287   void seed(uint32_t value = 331u)

   288   {

   289     minstd_rand0 intgen(value);

   290     seed(intgen);

   291   }

   292 

   293   // For GCC, moving this function out-of-line prevents inlining, which may

   294   // reduce overall object code size.  However, MSVC does not grok

   295   // out-of-line template member functions.

   296   template<class Generator>

   297   void seed(Generator & gen)

   298   {

   299     // use pass-by-reference, but wrap argument in pass_through_engine

   300     typedef detail::pass_through_engine<Generator&> ref_gen;

   301     uniform_01<ref_gen, RealType> gen01 =

   302       uniform_01<ref_gen, RealType>(ref_gen(gen));

   303     // I could have used std::generate_n, but it takes "gen" by value

   304     for(unsigned int j = 0; j < long_lag; ++j)

   305       x[j] = gen01();

   306     i = long_lag;

   307   }

   308 

   309   template<class It>

   310   void seed(It& first, It last)

   311   {

   312 #ifndef BOOST_NO_STDC_NAMESPACE

   313     // allow for Koenig lookup

   314     using std::fmod;

   315     using std::pow;

   316 #endif

   317     unsigned long mask = ~((~0u) << (w%32));   // now lowest w bits set

   318     RealType two32 = pow(RealType(2), 32);

   319     unsigned int j;

   320     for(j = 0; j < long_lag && first != last; ++j, ++first) {

   321       x[j] = RealType(0);

   322       for(int i = 0; i < w/32 && first != last; ++i, ++first)

   323         x[j] += *first / pow(two32,i+1);

   324       if(first != last && mask != 0)

   325         x[j] += fmod((*first & mask) / _modulus, RealType(1));

   326     }

   327     i = long_lag;

   328     if(first == last && j < long_lag)

   329       throw std::invalid_argument("lagged_fibonacci_01::seed");

   330   }

   331 

   332   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }

   333   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }

   334 

   335   result_type operator()()

   336   {

   337     if(i >= long_lag)

   338       fill();

   339     return x[i++];

   340   }

   341 

   342   static bool validation(result_type x)

   343   {

   344     result_type v = fibonacci_validation<result_type, p, q>::value();

   345     result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();

   346     // std::abs is a source of trouble: sometimes, it's not overloaded

   347     // for double, plus the usual namespace std noncompliance -> avoid it

   348     // using std::abs;

   349     // return abs(x - v) < 5 * epsilon

   350     return x > v - epsilon && x < v + epsilon;

   351   }

   352   

   353 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE

   354 

   355 #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS

   356   template<class CharT, class Traits>

   357   friend std::basic_ostream<CharT,Traits>&

   358   operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)

   359   {

   360 #ifndef BOOST_NO_STDC_NAMESPACE

   361     // allow for Koenig lookup

   362     using std::pow;

   363 #endif

   364     os << f.i << " ";

   365     std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); 

   366     for(unsigned int i = 0; i < f.long_lag; ++i)

   367       os << f.x[i] * f._modulus << " ";

   368     os.flags(oldflags);

   369     return os;

   370   }

   371 

   372   template<class CharT, class Traits>

   373   friend std::basic_istream<CharT, Traits>&

   374   operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)

   375     {

   376 # ifdef BOOST_RANDOM_EXTRACT_LF

   377         return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);

   378 # else

   379         is >> f.i >> std::ws;

   380         for(unsigned int i = 0; i < f.long_lag; ++i) {

   381             typename lagged_fibonacci_01::result_type value;

   382             is >> value >> std::ws;

   383             f.x[i] = value / f._modulus;

   384         }

   385         return is;

   386 # endif 

   387     }

   388 #endif

   389 

   390   friend bool operator==(const lagged_fibonacci_01& x,

   391                          const lagged_fibonacci_01& y)

   392   { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }

   393   friend bool operator!=(const lagged_fibonacci_01& x,

   394                          const lagged_fibonacci_01& y)

   395   { return !(x == y); }

   396 #else

   397   // Use a member function; Streamable concept not supported.

   398   bool operator==(const lagged_fibonacci_01& rhs) const

   399   { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }

   400   bool operator!=(const lagged_fibonacci_01& rhs) const

   401   { return !(*this == rhs); }

   402 #endif

   403 

   404 private:

   405   void fill();

   406   unsigned int i;

   407   RealType x[long_lag];

   408   RealType _modulus;

   409 };

   410 

   411 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION

   412 //  A definition is required even for integral static constants

   413 template<class RealType, int w, unsigned int p, unsigned int q>

   414 const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;

   415 template<class RealType, int w, unsigned int p, unsigned int q>

   416 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;

   417 template<class RealType, int w, unsigned int p, unsigned int q>

   418 const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;

   419 template<class RealType, int w, unsigned int p, unsigned int q>

   420 const int lagged_fibonacci_01<RealType,w,p,q>::word_size;

   421 

   422 #endif

   423 

   424 template<class RealType, int w, unsigned int p, unsigned int q>

   425 void lagged_fibonacci_01<RealType, w, p, q>::fill()

   426 {

   427   // two loops to avoid costly modulo operations

   428   {  // extra scope for MSVC brokenness w.r.t. for scope

   429   for(unsigned int j = 0; j < short_lag; ++j) {

   430     RealType t = x[j] + x[j+(long_lag-short_lag)];

   431     if(t >= RealType(1))

   432       t -= RealType(1);

   433     x[j] = t;

   434   }

   435   }

   436   for(unsigned int j = short_lag; j < long_lag; ++j) {

   437     RealType t = x[j] + x[j-short_lag];

   438     if(t >= RealType(1))

   439       t -= RealType(1);

   440     x[j] = t;

   441   }

   442   i = 0;

   443 }

   444 

   445 } // namespace random

   446 

   447 typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;

   448 typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;

   449 typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;

   450 typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;

   451 typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;

   452 typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;

   453 typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;

   454 typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;

   455 typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;

   456 

   457 

   458 // It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>

   459 // to help the compiler generate efficient code.  For GCC, this seems useless,

   460 // because GCC optimizes (x-0)/(1-0) to (x-0).  This is good enough for now.

   461 

   462 } // namespace boost

   463 

   464 #endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP