epoc32/include/stdapis/boost/random/lagged_fibonacci.hpp
author William Roberts <williamr@symbian.org>
Tue, 16 Mar 2010 16:12:26 +0000
branchSymbian2
changeset 2 2fe1408b6811
permissions -rw-r--r--
Final list of Symbian^2 public API header files
     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