//ranGen.h
#ifndef RANGEN_H
#define RANGEN_H

enum unifType {FM2, WH};
enum normType {BM, C};

class ranGen{
  
  unifType T;
  normType NT;
  unsigned long seed1, seed2, seed3;
  
  double gCauchy(double x) const;
  double GinvCauchy(double u) const;
  double fNorm(double x) const;

 public:

 ranGen(unifType t = WH, normType nt = BM) : T(t), NT(nt) {}
  ~ranGen(){};

  void ranUnif(int n, double* pu);
  void ranNorm(int n, double* pNorm, double mu = 0, 
	       double sig = 1);

  void setSeed(unsigned long s1);
};
#endif

//ranGen.cpp
#include "ranGen.h"
#include <iostream>
#include <cmath>

//anonymous namespace
namespace {
  double pi = 2.*acos(0.);
  double sqrtwopi = sqrt(2.*pi);
  unsigned long a1 = 171, m1 = 30269, a2 = 172, m2 = 30307,
    a3 = 170, m3 = 30323, m4 = 2147483647, a4 = 630360016;
}

void ranGen::ranUnif(int n, double* pu) {
  if(T == WH){    
    double temp = 0;
    for(int i=0; i< n; i++){
      seed1=(a1*seed1)%m1;
      seed2=(a2*seed2)%m2;
      seed3=(a3*seed3)%m3;
      temp = (((double)seed1)/((double)m1) + ((double)seed2)/((double)m2)
              + ((double)seed3)/((double)m3));
      if(temp > 2.) temp -= 2.;
      if(temp > 1.) temp -= 1.;
      pu[i] = temp;
    }
  }
  else{
    for(int i=0; i< n; i++){
      seed1=(a4*seed1)%m4;
      pu[i]=((double)seed1)/((double)m4);
    }
  }
}

void ranGen::ranNorm(int n, double* pNorm, double mu, 
		     double sig) {
  double* pu = new double(2);
  int nFill = 0, iAccept = 0, nReject = 0;
  
  if(NT == BM){
    double spare = 0, mult = 0, Rsqr = 0;
    
    while(nFill < n){
      if(iAccept == 1){//if we have one left, use it now
        pNorm[nFill] = spare;
        iAccept = 0;//next time we need to generate another pair 
        
        nFill += 1;
      }
      else{
        while(iAccept == 0){
          
          //generate two uniforms on [-1, 1]
          ranUnif(2,pu);
          pu[0] = 2*pu[0]- 1;
          pu[1] = 2*pu[1] - 1;
          
          Rsqr = pu[0]*pu[0] + pu[1]*pu[1];
          if(Rsqr < 1 && Rsqr != 0){//accept the result
            mult = sqrt(-2.*log(Rsqr)/Rsqr);
            pNorm[nFill] = mu + sig*pu[0]*mult;
            spare = mu + sig*pu[1]*mult;//save this one for next time
            nFill += 1;
            iAccept = 1;//now we have one left over
          }
          else nReject += 2;
        }
      }
    }
  }	
  
  else{
    double y;
    double c = sqrt(2*pi/exp(1.));
    
    while(nFill < n){
      iAccept = 0;
      while(iAccept == 0){
        ranUnif(2, pu);
        y = GinvCauchy(pu[1]);
        if(c*gCauchy(y)*pu[0] <= fNorm(y)){//accept the result
          pNorm[nFill] = mu + sig*y;
          nFill += 1;//augment the number of samples
          iAccept = 1;//end the interior while loop
        }
        else nReject += 1;
      }
    }
  }	
}

void ranGen::setSeed(unsigned long s1){
  if(T == WH){
    seed1 = s1;
    seed2=(a4*seed1)%m4;
    seed3=(a4*seed2)%m4;
  }
  else 
    seed1 = s1;
}
	
double ranGen::gCauchy(double x) const {
  return 1./((1. + x*x)*pi);
}

double ranGen::GinvCauchy(double u) const {
  return tan(pi*(u - .5));
}

double ranGen::fNorm(double x) const {
  return exp(-.5*x*x)/sqrtwopi;
}