///////////////////////////////////////////////////////////////////////
//basicFM2.cpp
#include <iostream>

using std::cout; using std::endl;

unsigned long power(int i, int k){
  unsigned long iout = (unsigned long)i;
  if(k > 1)
    iout = iout*power(iout, k - 1);
  return iout; 
}

int main(int argc, char* argv[]){
  unsigned long a = 630360016;//multiplier
  unsigned long seed=123;//seed
  unsigned long m = power(2,31) - 1;//modulus
  cout << "Your random number is ";
  seed = (a*seed)%m;
  double u = ((double)seed)/((double)m);
  cout << u << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//basicFM2.cpp:with alternative main function
#include <iostream>

using std::cout; using std::endl;

unsigned long power(int i, int k){
  unsigned long iout = (unsigned long)i;
  if(k > 1)
    iout = iout*power(iout, k - 1);
  return iout; 
}

int main(int argc, char* argv[]){
  unsigned long a = 630360016;//multiplier
  unsigned long seed = 0;
  int n = 0;
  cout << "Choose a value for the sample size:" << endl;
  std::cin >> n;
  cout << "Choose a value for the seed:" << endl;
  std::cin >> seed;
  unsigned long m = power(2, 31) - 1;

  cout <<  "Your random numbers are " << endl;
  for(int i = 0; i <  n; i++){
    seed = (a*seed)%m;
    double u = ((double)seed)/((double)m);
    cout << " " << u;
  }

  /*
  double u[n];
  std::cout << "Your random numbers are" << std::endl;
  for(int i = 0; i < n; i++){
    seed = (a*seed)%m;
    u[i] = ((double)seed)/((double)m);
    std::cout << "u[" << i <<"]= " << u[i] << " ";
  }

 double* pU = new double[n];
  std::cout << "Your random numbers are:" << std::endl;
  for(int i = 0; i < n; i++){
    seed = (a*seed)%m;
    pU[i] = ((double)seed)/((double)m);
    std::cout << "u[" << i <<"]= " << pU[i] << " ";
  }
  */ 
  cout << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//FM2.h
#ifndef FM2_H
#define FM2_H

class FM2{
  unsigned long a;//multiplier
  int expFor2;//exponent for 2 in modulus
  unsigned long power(int i, int k) const;
  
public:

  FM2();//default constructor
  ~FM2(){};//destructor
  void rangen(unsigned long seed, int n, double* pu) const;	
};
#endif

//FM2.cpp
#include "FM2.h"

unsigned long FM2::power(int i, int k) const {
  unsigned long iout = (unsigned long)i;
  if(k > 1)
    iout = iout*power(iout, k - 1);    
  return iout; 
}

FM2::FM2(){
  a = 630360016; expFor2 = 31;
}
      
void FM2::rangen(unsigned long seed, int n, double* pU) const {
  unsigned long m = power(2,expFor2)-1;
  for(int i = 0; i <  n; i++){
    seed = (a*seed)%m;
    pU[i] = ((double)seed)/((double)m);
  }
}

//driverFM2.cpp
#include <iostream>
#include <cstdlib>
#include <new>
#include "FM2.h"

using std::cout; using std::cin; using std::endl;

int main(int argc, char* argv[]){
  int n=0;
  unsigned long seed = 0;
  cout << "Choose a value for the sample size:" << endl;
  cin >> n;
  cout << "Choose a value for the seed:" << endl;
  cin >> seed;
  FM2* pFM2 = new(std::nothrow) FM2();
  if(pFM2 == 0){
    cout<< "Memory allocation for pFM2 failed"<<endl;
    exit(1);
  }
  double* pU = new(std::nothrow) double[n];
  if(pU == 0){
    cout<< "Memory allocation for pu failed"<<endl;
    exit(1);
  }
  pFM2->rangen(seed, n, pU);
  cout << "Your random numbers are" << endl;
  for(int i = 0; i < n; i++)
    cout << "u[" << i << "]= "<< pU[i] << " ";
  cout << endl;

  delete pFM2;
  delete[] pU;
  return 0;
}	

driverFM2 : FM2.o driverFM2.o 
	g++ -Wall FM2.o driverFM2.o -o driverFM2 

FM2.o : FM2.cpp FM2.h
	g++ -Wall -c FM2.cpp

driverFM2.o : driverFM2.cpp FM2.h
	g++ -Wall -c driverFM2.cpp
///////////////////////////////////////////////////////////////////////
//gfsr.cpp
#include <iostream>

void intToBin(unsigned long val, bool* pb, int r){
  for(int i = 0; i < r; i++){
    if(val & (1 << i))
      pb[r - i - 1] = 1;
    else
      pb[r - i - 1] = 0;
  }
}

unsigned long power(int i, int k){
  if(k == 0)
    return 1;
  else{
    int iout = i;
    if(k > 1)
      iout = iout*power(iout, k - 1);  
    return iout; 
  }
}

unsigned long binToInt(bool* pb, int r){
  unsigned long val = 0;
  for(int i = 0; i < r; i++)
    val += pb[i]*power(2, r - i - 1);
  return val;
}

void binMult(bool** pb, int p, int q, int r){
  //use the extra row of pb here
  for(int j = 0; j < r; j++){
    pb[p][j] = pb[0][j]^pb[q][j];
  }
  //Update pb
  for(int i = 0; i < p; i++)
    for(int j = 0; j < r; j++)
      pb[i][j] = pb[i + 1][j];
}

int main(){
  int r;
  std::cout << "Input the value for r:" << std::endl;
  std::cin >> r;
  int p;
  std::cout << "Input the value for p:" << std::endl;
  std::cin >> p;
  int q;
  std::cout << "Input the value for q:" << std::endl;
  std::cin >> q;
  unsigned long* x0 = new unsigned long[p];
  std::cout << "Input a " << p << " integer seed:" << std::endl;
  for(int i = 0; i < p; i++)
    std::cin >> x0[i];
  
  unsigned long modulus = power(2, p) - 1;

  //matrix of Booleans to hold the bits
  bool** pb = new bool*[p + 1];
  //rows will hold the integer bit-wise representations
  //extra row added for temporary storage
  for(int i = 0; i < (p + 1); i++)
    pb[i] = new bool[r];
  
  //initialize the bit pattern using x0
  for(int i = 0; i < p; i++)
    intToBin(x0[i], pb[i], r);
  
  //now generate the "data"
  for(int i = 0; i < (modulus - p); i++){ 
    binMult(pb, p, q, r); 
    std::cout << " " << binToInt(pb[p - 1], r) << " ";
  }
  std::cout << std::endl;
    
  return 0;
}
///////////////////////////////////////////////////////////////////////
//WH.h
#ifndef WH_H
#define WH_H

class WH{
  unsigned long a1, a2, a3;
  unsigned long m1, m2, m3;
  
 public:

  WH();
  ~WH(){};
  void rangen(unsigned long seed1, unsigned long seed2, 
	      unsigned long seed3, int n, double* pu) const;	
};
#endif

// WH.cpp
#include "WH.h"

WH::WH(){
  a1 = 171; m1 = 30269; a2 = 172; m2 = 30307; a3 = 170; m3 = 30323;
}

void WH::rangen(unsigned long seed1, unsigned long seed2, 
		unsigned long seed3, int n, double* pu) const {
  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;
  }
}
///////////////////////////////////////////////////////////////////////
//rngStreamEx.cpp
#include <iostream>
#include "RngStream.h"

using std::cout; using std::endl;

int main(){
  unsigned long seed[6] = {1, 2, 3, 4, 5, 6};
  RngStream::SetPackageSeed (seed);
  RngStream rngObj;
  cout << "Your random number is ";
  cout << rngObj.RandU01() << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//driverExp.cpp
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <new>
#include "WH.h"

using std::cout; using std::cin; using std::endl;

double Finv(double u, double lambda){
  return (-log(1 - u)/lambda);
}

int main(int argc, char* argv[]){
  double lambda; int n;
    
  cout << "Input the value for lambda:" << endl;
  cin >> lambda;
    
  unsigned long seed1, seed2, seed3;
  cout << "Choose a value for the sample size:" << endl;
  cin >> n;
  cout << "Choose a value for seed1:" << endl;
  cin >> seed1;
  cout << "Choose a value for seed2:" << endl;
  cin >> seed2;
  cout << "Choose a value for seed3:" << endl;
  cin >> seed3;

  WH* pWH = new(std::nothrow) WH();
  if(pWH == 0){
    cout << "Memory allocation for WH failed" <<endl;
    exit(1);
  }
  double* pU = new(std::nothrow) double[n];
  if(pU == 0){
    cout <<"Memory allocation for double pointer failed" <<endl;
    exit(1);
  }
  pWH->rangen(seed1, seed2, seed3, n, pU);
  cout << "Your random numbers are " << endl;
  for(int i=0; i< n; i++)
    cout << Finv(pU[i], lambda) << " ";
  cout << endl;

  delete pWH;
  delete[] pU;
  return 0;
}
///////////////////////////////////////////////////////////////////////
// WHNCauchy.h
#ifndef WHNCAUCHY_H
#define WHNCAUCHY_H

class WHNCauchy{
  unsigned long a1, a2, a3;
  unsigned long m1, m2, m3;
  double g(double x) const;
  double Ginv(double x) const;
  double f(double x) const;

 public:

  WHNCauchy();
  ~WHNCauchy(){};
  void rangen(unsigned long seed1, unsigned long seed2, 
	      unsigned long seed3, int n, double* pNorm) const;
};
#endif

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

using std::cout; using std::endl;

double pi = 2.*acos(0.);
double sqrtwopi = sqrt(2.*pi);
double c = sqrt(2*pi/exp(1.));

WHNCauchy::WHNCauchy(){
  a1 = 171; m1 = 30269; a2 = 172; m2 = 30307; a3 = 170; m3 = 30323;
}

void WHNCauchy::rangen(unsigned long seed1, 
		       unsigned long seed2, unsigned long seed3, 
		       int n, double* pNorm) const {
  double temp, y;
  double u[2];
  int nFill = 0, iAccept, nReject = 0;
 
  while(nFill < n){
    iAccept = 0;
    while(iAccept == 0){
      //generate two uniforms
      for(int i = 0; i < 2; 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.;
        u[i] = temp;
      }
      y = Ginv(u[1]);
      if(c*g(y)*u[0] <= f(y)){//accept the result
        pNorm[nFill] = y;
        nFill += 1;//augment the number of samples
        iAccept = 1;//end the interior while loop
      }
      else nReject += 1;
    }
  }
  cout << "The Number of rejections was " << nReject << endl; 
}	
	
double WHNCauchy::g(double x) const {
  return 1/((1 + x*x)*pi);
}

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

double WHNCauchy::f(double x) const {
  return exp(-.5*x*x)/sqrtwopi;
}
///////////////////////////////////////////////////////////////////////
//WHBM.h
#ifndef WHN_H
#define WHN_H

class WHBM{
  unsigned long a1, a2, a3;
  unsigned long m1, m2, m3;
  
 public:

  WHBM();
  ~WHBM(){};
  void rangen(unsigned long seed1, unsigned long seed2, 
        unsigned long seed3, int n, double* pu) const;	
};
#endif

//WHBM.cpp
#include "WHBM.h"
#include <cmath>

WHBM::WHBM(){
  a1 = 171; m1 = 30269; a2 = 172; m2 = 30307; a3 = 170; m3 = 30323;
}

void WHBM::rangen(unsigned long seed1, unsigned long seed2,
        unsigned long seed3, int n, double* pNorm) const {
  double temp = 0, spare = 0, mult = 0, u[2] = {0, 0}, Rsqr = 0;
  
  int nFill = 0, iAccept = 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
        for(int i = 0; i < 2; 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.;
          u[i] = 2*temp - 1; //translate into [-1, 1]
        }
        Rsqr = u[0]*u[0] + u[1]*u[1];
        if(Rsqr < 1 && Rsqr != 0){//accept the result
          mult = sqrt(-2.*log(Rsqr)/Rsqr);
          pNorm[nFill] = u[0]*mult;
          spare = u[1]*mult;//save this one for next time
          nFill += 1;
          iAccept = 1;//now we have one left over
        } } } } }	
///////////////////////////////////////////////////////////////////////
//RmathEx1.cpp
#define MATHLIB_STANDALONE 
#include <iostream>
#include "Rmath.h"

using std::cout; using std::cin; using std::endl;

int main()
{
  double x1, x2;
  cout << "Enter a argument for the normal cdf:" << endl;
  cin >> x1;
  cout << "Enter a argument for the chi-squared cdf:" << endl;
  cin >> x2;

  cout << "Prob(Z <= " << x1 << ") = " << 
    pnorm(x1, 0, 1, 1, 0)  << endl;
  cout << "Prob(Chi^2 <= " << x2 << ")= " << 
    pchisq(x2, 1, 1, 0) << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//RmathEx2.cpp
#define MATHLIB_STANDALONE 
#include <iostream>
#include "Rmath.h"

using std::cin; using std::cout; using std::endl;

int main()
{
  unsigned int seed1, seed2, n, df;
  cout << "Enter a value for the sample size:" << endl;
  cin >> n;
  cout << "Enter a value for the first seed:" << endl;
  cin >> seed1;
  cout << "Enter a value for the second seed:" << endl;
  cin >> seed2;
  cout << "Enter a value for degrees-of-freedom:" << endl;
  cin >> df;

  set_seed(seed1,seed2);
  cout << "Here is a sample of standard normals" << endl;
  for(int i = 0; i < n; i++)
    cout << rnorm(0, 1) << " ";
  cout << endl;
  cout << "Here is a sample of chi-squares with " << df
       << " degrees-of-freedom" << endl;
  for(int i = 0; i < n; i++)
    cout << rchisq(df) << " ";
  cout << endl;
  
  unsigned int seed3, seed4;
  get_seed(&seed3, &seed4);
  cout << "The first seed is now " << seed3 << endl;
  cout << "The second seed is now " << seed4 << endl;

  set_seed(seed1,seed2);
  cout << "Here is your sample of standard normals again " << endl;
  for(int i = 0; i < n; i++)
    cout << rnorm(0, 1) << " ";
  cout << endl;

  return 0;
}