/*
 * Author: Sebastian Weber 
 *
 * uniform random number generator
 * based on a four-tap shift-register algorithm
 *
 * implements Boost UniformRandomNumberGenerator concept
 *
 * some code pieces are taken from the gsl-library
 *
 */

#ifndef FOUR_TAP2_HPP
#define FOUR_TAP2_HPP

#include <boost/static_assert.hpp>

namespace boost {
  namespace random {

#include <stdexcept>
#include <iostream>
#include <algorithm>


template<class UIntType>
class four_tap {
public:
  typedef four_tap<UIntType> self;
  typedef UIntType result_type;

  // enfore unsigned integers
  BOOST_STATIC_ASSERT(std::numeric_limits<UIntType>::is_integer && !std::numeric_limits<UIntType>::is_signed);

  four_tap() { seed(); }
  four_tap(result_type s) { seed(s);  }

  // copy constructor
  four_tap(const self& ft) {
    // copy over random data
    std::copy(&ft.m_random_array[0], &ft.m_random_array[m_size], &m_random_array[0]);
    // reeint the pointers correctly
    m_random_which1 = &m_random_array[0] - (ft.m_random_array - ft.m_random_which1);
    m_random_which2 = &m_random_array[0] - (ft.m_random_array - ft.m_random_which2);
    m_random_which3 = &m_random_array[0] - (ft.m_random_array - ft.m_random_which3);
    m_random_which4 = &m_random_array[0] - (ft.m_random_array - ft.m_random_which4);
    m_random_which5 = &m_random_array[0] - (ft.m_random_array - ft.m_random_which5);
    m_random_end    = &m_random_array[m_size];
  }

  // copy assignment
  self& operator=(const self& ft) { 
    return self(ft);
  }

  void seed() { seed(result_type(4357)); }

  void seed(result_type s);

  result_type operator()();

   static const bool has_fixed_range = false;

  // actually we have a fixed range, but that depends upon the UIntType
  //static const bool has_fixed_range = true;
  //static const result_type min_value = std::numeric_limits<UIntType>::min();;
  //static const result_type max_value = std::numeric_limits<UIntType>::max();

  result_type min() const { return std::numeric_limits<UIntType>::min(); }
  result_type max() const { return std::numeric_limits<UIntType>::max(); }

private:
  static const result_type m_multiplier = result_type(16807);
  static const std::size_t m_size = 16384;

  void initPtrs() {
    m_random_which1 = &m_random_array[0];
    m_random_which2 = &m_random_array[471];
    m_random_which3 = &m_random_array[1586];
    m_random_which4 = &m_random_array[6988];
    m_random_which5 = &m_random_array[9689];
    m_random_end    = &m_random_array[m_size];
  }

  // the array of random numbers
  result_type m_random_array[m_size];

  // next random numbers
  result_type *m_random_which1, *m_random_which2, *m_random_which3, *m_random_which4, *m_random_which5;

  // end of array
  result_type *m_random_end;
};

template<class UIntType>
void four_tap<UIntType>::seed(UIntType s) {
  const std::size_t digits = static_cast<std::size_t>(std::numeric_limits<UIntType>::digits);

  /*
   * init-code taken from gsl, adapted to be compatible with 32bit and
   * 64bit integers
   */

  /* Masks for turning on the diagonal bit and turning off the
     leftmost bits */
  UIntType msb = UIntType(1); msb <<= digits-1;
  UIntType mask = std::numeric_limits<UIntType>::max();


  /* We use the congruence s_{n+1} = (69069*s_n) mod 2^32 to
     initialize the state. This works because ANSI-C unsigned long
     integer arithmetic is automatically modulo 2^32 (or a higher
     power of two), so we can safely ignore overflow.
  */

#define LCG(n) ((UIntType(69069) * n) & std::numeric_limits<UIntType>::max())
  //#define LCG(n) ((UIntType(69069) * n) & 0xffffffffUL)

  /* Brian Gough suggests this to avoid low-order bit correlations */
  for (std::size_t i = 0; i < m_size; ++i) {
    UIntType t = 0 ;
    UIntType bit = msb ;
    for (std::size_t j = 0; j < std::size_t(digits); ++j) {
      s = LCG(s);
      if (s & msb) 
	t |= bit ;
      bit >>= 1 ;
    }
    m_random_array[i] = t;
  }

  /* Perform the "orthogonalization" of the matrix */
  /* Based on the orthogonalization used in r250, as suggested initially
   * by Kirkpatrick and Stoll, and pointed out to me by Brian Gough
   */

  /* BJG: note that this orthogonalisation doesn't have any effect
     here because the the initial 6695 elements do not participate in
     the calculation.  For practical purposes this orthogonalisation
     is somewhat irrelevant, because the probability of the original
     sequence being degenerate should be exponentially small. */

  for (std::size_t i=0; i<std::size_t(digits); ++i) {
    std::size_t k=7+i*3;
    m_random_array[k] &= mask;     /* Turn off bits left of the diagonal */
    m_random_array[k] |= msb;      /* Turn on the diagonal bit           */
    mask >>= 1;
    msb >>= 1;
  }

  // set pointers correct
  initPtrs();	
}


template<class UIntType>
UIntType four_tap<UIntType>::operator()() {

  /* create random number */
  if( ++m_random_which1 == m_random_end) {
    m_random_which1= m_random_array;
  }
  if( ++m_random_which2 == m_random_end) {
    m_random_which2= m_random_array;
  }
  if( ++m_random_which3 == m_random_end) {
    m_random_which3= m_random_array;
  }
  if( ++m_random_which4 == m_random_end) {
    m_random_which4= m_random_array;
  }
  if( ++m_random_which5 == m_random_end) {
    m_random_which5= m_random_array;
  }

  return( *m_random_which1=
  	  *m_random_which2 ^ *m_random_which3 ^ *m_random_which4 ^ *m_random_which5 );
}

  }; /* random namespace */

  typedef random::four_tap<uint32_t> four_tap32_t;
  typedef random::four_tap<uint64_t> four_tap64_t;
  typedef four_tap32_t four_tap_t;

}; /* boost namespace */

#endif /* FOUR_TAP2_HPP */