#ifndef BOOST_RANDOM_NUMBERS_HPP
#define BOOST_RANDOM_NUMBERS_HPP
#include <boost/random.hpp>
#include <ctime> 
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
//typedef boost::mt19937 rng;
typedef boost::lagged_fibonacci607 rng;
rng gener;
void initrng(rng & gen) {
  // initialize generator
  static bool first=true;
  if(first) {
    std::cout<<"  initrng: initializing generator\n";
    gen.seed(static_cast<boost::uint32_t> (std::time(0)));
    first = false;
  }
}
double unirand(void){
  // U[0,1)
  static bool first=true;
  if(first) {
    initrng(gener);
    first =false;
  }
  typedef boost::uniform_real<> dist_type;
  dist_type dist(0,1);
  boost::variate_generator<rng&, dist_type> random(gener, dist);
  return random();
}
double gaussrand(){
  // Gaussian (Normal) distribution
  static bool first=true;
  if(first) {
    initrng(gener);
    first =false;
  }
  typedef boost::normal_distribution<> dist_type;
  dist_type  dist(0.0,1.0);
  boost::variate_generator<rng&, dist_type> random(gener, dist);
  return random();
}
#endif