/* boost random/lagged_fibonacci.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id: lagged_fibonacci.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $
 *
 * Revision history
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
#define BOOST_RANDOM_LAGGED_FIBONACCI_HPP

#include <boost/config/no_tr1/cmath.hpp>
#include <iostream>
#include <algorithm>     // std::max
#include <iterator>
#include <boost/config/no_tr1/cmath.hpp>         // std::pow
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/cstdint.hpp>
#include <boost/detail/workaround.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/pass_through_engine.hpp>

namespace boost {
namespace random {

#if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
#  define BOOST_RANDOM_EXTRACT_LF
#endif

#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
#  define BOOST_RANDOM_EXTRACT_LF
#endif

#  ifdef BOOST_RANDOM_EXTRACT_LF
namespace detail
{
  template<class IStream, class F, class RealType>
  IStream&
  extract_lagged_fibonacci_01(
      IStream& is
      , F const& f
      , unsigned int& i
      , RealType* x
      , RealType modulus)
  {
        is >> i >> std::ws;
        for(unsigned int i = 0; i < f.long_lag; ++i)
        {
            RealType value;
            is >> value >> std::ws;
            x[i] = value / modulus;
        }
        return is;
  }

  template<class IStream, class F, class UIntType>
  IStream&
  extract_lagged_fibonacci(
      IStream& is
      , F const& f
      , unsigned int& i
      , UIntType* x)
  {
      is >> i >> std::ws;
      for(unsigned int i = 0; i < f.long_lag; ++i)
          is >> x[i] >> std::ws;
      return is;
  }
}
#  endif

/** 
 * Instantiations of class template \lagged_fibonacci model a
 * \pseudo_random_number_generator. It uses a lagged Fibonacci
 * algorithm with two lags @c p and @c q:
 * x(i) = x(i-p) + x(i-q) (mod 2<sup>w</sup>) with p > q.
 */
template<class UIntType, int w, unsigned int p, unsigned int q,
         UIntType val = 0>
class lagged_fibonacci
{
public:
  typedef UIntType result_type;
  BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  BOOST_STATIC_CONSTANT(int, word_size = w);
  BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
  BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);

  /**
   * Returns: the smallest value that the generator can produce
   */
  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
  /**
   * Returns: the largest value that the generator can produce
   */
  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }

  /**
   * Creates a new @c lagged_fibonacci generator and calls @c seed()
   */
  lagged_fibonacci() { init_wordmask(); seed(); }
  /**
   * Creates a new @c lagged_fibonacci generator and calls @c seed(value)
   */
  explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }
  /**
   * Creates a new @c lagged_fibonacci generator and calls @c seed(first, last)
   */
  template<class It> lagged_fibonacci(It& first, It last)
  { init_wordmask(); seed(first, last); }
  // compiler-generated copy ctor and assignment operator are fine

private:
  /// \cond hide_private_members
  void init_wordmask()
  {
    wordmask = 0;
    for(int j = 0; j < w; ++j)
      wordmask |= (1u << j);
  }
  /// \endcond

public:
  /**
   * Sets the state of the generator to values produced by
   * a \minstd_rand generator.
   */
  void seed(uint32_t value = 331u)
  {
    minstd_rand0 gen(value);
    for(unsigned int j = 0; j < long_lag; ++j)
      x[j] = gen() & wordmask;
    i = long_lag;
  }

  /**
   * Sets the state of the generator to values from the iterator
   * range [first, last).  If there are not enough elements in the
   * range [first, last) throws @c std::invalid_argument.
   */
  template<class It>
  void seed(It& first, It last)
  {
    // word size could be smaller than the seed values
    unsigned int j;
    for(j = 0; j < long_lag && first != last; ++j, ++first)
      x[j] = *first & wordmask;
    i = long_lag;
    if(first == last && j < long_lag)
      throw std::invalid_argument("lagged_fibonacci::seed");
  }

  /**
   * Returns: the next value of the generator
   */
  result_type operator()()
  {
    if(i >= long_lag)
      fill();
    return x[i++];
  }

  static bool validation(result_type x)
  {
    return x == val;
  }
  
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  template<class CharT, class Traits>
  friend std::basic_ostream<CharT,Traits>&
  operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)
  {
    os << f.i << " ";
    for(unsigned int i = 0; i < f.long_lag; ++i)
      os << f.x[i] << " ";
    return os;
  }

  template<class CharT, class Traits>
  friend std::basic_istream<CharT, Traits>&
  operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
  {
# ifdef BOOST_RANDOM_EXTRACT_LF
      return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
# else
      is >> f.i >> std::ws;
      for(unsigned int i = 0; i < f.long_lag; ++i)
          is >> f.x[i] >> std::ws;
      return is;
# endif 
  }
#endif

  friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
  { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
  friend bool operator!=(const lagged_fibonacci& x,
                         const lagged_fibonacci& y)
  { return !(x == y); }
#else
  // Use a member function; Streamable concept not supported.
  bool operator==(const lagged_fibonacci& rhs) const
  { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
  bool operator!=(const lagged_fibonacci& rhs) const
  { return !(*this == rhs); }
#endif

private:
  /// \cond hide_private_members
  void fill();
  /// \endcond

  UIntType wordmask;
  unsigned int i;
  UIntType x[long_lag];
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
#endif

/// \cond hide_private_members

template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
void lagged_fibonacci<UIntType, w, p, q, val>::fill()
{
  // two loops to avoid costly modulo operations
  {  // extra scope for MSVC brokenness w.r.t. for scope
  for(unsigned int j = 0; j < short_lag; ++j)
    x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;
  }
  for(unsigned int j = short_lag; j < long_lag; ++j)
    x[j] = (x[j] + x[j-short_lag]) & wordmask;
  i = 0;
}


// lagged Fibonacci generator for the range [0..1)
// contributed by Matthias Troyer
// for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958

template<class T, unsigned int p, unsigned int q>
struct fibonacci_validation
{
  BOOST_STATIC_CONSTANT(bool, is_specialized = false);
  static T value() { return 0; }
  static T tolerance() { return 0; }
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class T, unsigned int p, unsigned int q>
const bool fibonacci_validation<T, p, q>::is_specialized;
#endif

#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
template<> \
struct fibonacci_validation<T, P, Q>  \
{                                     \
  BOOST_STATIC_CONSTANT(bool, is_specialized = true);     \
  static T value() { return V; }      \
  static T tolerance()                \
{ return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
};
// (The extra static_cast<T> in the std::max call above is actually
// unnecessary except for HP aCC 1.30, which claims that
// numeric_limits<double>::epsilon() doesn't actually return a double.)

BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)

#undef BOOST_RANDOM_FIBONACCI_VAL

/// \endcond

/**
 * Instantiations of class template @c lagged_fibonacci_01 model a
 * \pseudo_random_number_generator. It uses a lagged Fibonacci
 * algorithm with two lags @c p and @c q, evaluated in floating-point
 * arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See
 *
 *  @blockquote
 *  "Uniform random number generators for supercomputers", Richard Brent,
 *  Proc. of Fifth Australian Supercomputer Conference, Melbourne,
 *  Dec. 1992, pp. 704-706.
 *  @endblockquote
 *
 * @xmlnote
 * The quality of the generator crucially depends on the choice
 * of the parameters. User code should employ one of the sensibly
 * parameterized generators such as \lagged_fibonacci607 instead.
 * @endxmlnote
 *
 * The generator requires considerable amounts of memory for the storage
 * of its state array. For example, \lagged_fibonacci607 requires about
 * 4856 bytes and \lagged_fibonacci44497 requires about 350 KBytes.
 */
template<class RealType, int w, unsigned int p, unsigned int q>
class lagged_fibonacci_01
{
public:
  typedef RealType result_type;
  BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  BOOST_STATIC_CONSTANT(int, word_size = w);
  BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
  BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);

  /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(). */
  lagged_fibonacci_01() { init_modulus(); seed(); }
  /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(value). */
  BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01, uint32_t, value)
  { init_modulus(); seed(value); }
  /** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(gen). */
  BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(lagged_fibonacci_01, Generator, gen)
  { init_modulus(); seed(gen); }
  template<class It> lagged_fibonacci_01(It& first, It last)
  { init_modulus(); seed(first, last); }
  // compiler-generated copy ctor and assignment operator are fine

private:
  /// \cond hide_private_members
  void init_modulus()
  {
#ifndef BOOST_NO_STDC_NAMESPACE
    // allow for Koenig lookup
    using std::pow;
#endif
    _modulus = pow(RealType(2), word_size);
  }
  /// \endcond

public:
  /** Calls seed(331u). */
  void seed() { seed(331u); }
  /**
   * Constructs a \minstd_rand0 generator with the constructor parameter
   * value and calls seed with it. Distinct seeds in the range
   * [1, 2147483647) will produce generators with different states. Other
   * seeds will be equivalent to some seed within this range. See
   * \linear_congruential for details.
   */
  BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_01, uint32_t, value)
  {
    minstd_rand0 intgen(value);
    seed(intgen);
  }

  /**
   * Sets the state of this @c lagged_fibonacci_01 to the values returned
   * by p invocations of \uniform_01<code>\<RealType\>()(gen)</code>.
   *
   * Complexity: Exactly p invocations of gen.
   */
  BOOST_RANDOM_DETAIL_GENERATOR_SEED(lagged_fibonacci, Generator, gen)
  {
    // use pass-by-reference, but wrap argument in pass_through_engine
    typedef detail::pass_through_engine<Generator&> ref_gen;
    uniform_01<ref_gen, RealType> gen01 =
      uniform_01<ref_gen, RealType>(ref_gen(gen));
    // I could have used std::generate_n, but it takes "gen" by value
    for(unsigned int j = 0; j < long_lag; ++j)
      x[j] = gen01();
    i = long_lag;
  }

  template<class It>
  void seed(It& first, It last)
  {
#ifndef BOOST_NO_STDC_NAMESPACE
    // allow for Koenig lookup
    using std::fmod;
    using std::pow;
#endif
    unsigned long mask = ~((~0u) << (w%32));   // now lowest w bits set
    RealType two32 = pow(RealType(2), 32);
    unsigned int j;
    for(j = 0; j < long_lag && first != last; ++j) {
      x[j] = RealType(0);
      for(int k = 0; k < w/32 && first != last; ++k, ++first)
        x[j] += *first / pow(two32,k+1);
      if(first != last && mask != 0) {
        x[j] += fmod((*first & mask) / _modulus, RealType(1));
        ++first;
      }
    }
    i = long_lag;
    if(first == last && j < long_lag)
      throw std::invalid_argument("lagged_fibonacci_01::seed");
  }

  result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
  result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }

  result_type operator()()
  {
    if(i >= long_lag)
      fill();
    return x[i++];
  }

  static bool validation(result_type x)
  {
    result_type v = fibonacci_validation<result_type, p, q>::value();
    result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
    // std::abs is a source of trouble: sometimes, it's not overloaded
    // for double, plus the usual namespace std noncompliance -> avoid it
    // using std::abs;
    // return abs(x - v) < 5 * epsilon
    return x > v - epsilon && x < v + epsilon;
  }
  
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE

#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  template<class CharT, class Traits>
  friend std::basic_ostream<CharT,Traits>&
  operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)
  {
#ifndef BOOST_NO_STDC_NAMESPACE
    // allow for Koenig lookup
    using std::pow;
#endif
    os << f.i << " ";
    std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); 
    for(unsigned int i = 0; i < f.long_lag; ++i)
      os << f.x[i] * f._modulus << " ";
    os.flags(oldflags);
    return os;
  }

  template<class CharT, class Traits>
  friend std::basic_istream<CharT, Traits>&
  operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)
    {
# ifdef BOOST_RANDOM_EXTRACT_LF
        return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
# else
        is >> f.i >> std::ws;
        for(unsigned int i = 0; i < f.long_lag; ++i) {
            typename lagged_fibonacci_01::result_type value;
            is >> value >> std::ws;
            f.x[i] = value / f._modulus;
        }
        return is;
# endif 
    }
#endif

  friend bool operator==(const lagged_fibonacci_01& x,
                         const lagged_fibonacci_01& y)
  { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
  friend bool operator!=(const lagged_fibonacci_01& x,
                         const lagged_fibonacci_01& y)
  { return !(x == y); }
#else
  // Use a member function; Streamable concept not supported.
  bool operator==(const lagged_fibonacci_01& rhs) const
  { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
  bool operator!=(const lagged_fibonacci_01& rhs) const
  { return !(*this == rhs); }
#endif

private:
  /// \cond hide_private_members
  void fill();
  /// \endcond
  unsigned int i;
  RealType x[long_lag];
  RealType _modulus;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class RealType, int w, unsigned int p, unsigned int q>
const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
template<class RealType, int w, unsigned int p, unsigned int q>
const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
template<class RealType, int w, unsigned int p, unsigned int q>
const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
template<class RealType, int w, unsigned int p, unsigned int q>
const int lagged_fibonacci_01<RealType,w,p,q>::word_size;

#endif

/// \cond hide_private_members
template<class RealType, int w, unsigned int p, unsigned int q>
void lagged_fibonacci_01<RealType, w, p, q>::fill()
{
  // two loops to avoid costly modulo operations
  {  // extra scope for MSVC brokenness w.r.t. for scope
  for(unsigned int j = 0; j < short_lag; ++j) {
    RealType t = x[j] + x[j+(long_lag-short_lag)];
    if(t >= RealType(1))
      t -= RealType(1);
    x[j] = t;
  }
  }
  for(unsigned int j = short_lag; j < long_lag; ++j) {
    RealType t = x[j] + x[j-short_lag];
    if(t >= RealType(1))
      t -= RealType(1);
    x[j] = t;
  }
  i = 0;
}
/// \endcond

} // namespace random

#ifdef BOOST_RANDOM_DOXYGEN
namespace detail {
/**
 * The specializations lagged_fibonacci607 ... lagged_fibonacci44497
 * use well tested lags.
 *
 * See
 *
 *  @blockquote
 *  "On the Periods of Generalized Fibonacci Recurrences", Richard P. Brent
 *  Computer Sciences Laboratory Australian National University, December 1992
 *  @endblockquote
 *
 * The lags used here can be found in
 *
 *  @blockquote
 *  "Uniform random number generators for supercomputers", Richard Brent,
 *  Proc. of Fifth Australian Supercomputer Conference, Melbourne,
 *  Dec. 1992, pp. 704-706.
 *  @endblockquote
 */
struct lagged_fibonacci_doc {};
}
#endif

/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
/**
 * @copydoc boost::detail::lagged_fibonacci_doc
 */
typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;


// It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>
// to help the compiler generate efficient code.  For GCC, this seems useless,
// because GCC optimizes (x-0)/(1-0) to (x-0).  This is good enough for now.

} // namespace boost

#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP