#include <math.h>
#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h> // use builtin rand() instead of using "_random"
#include <vector>
#include <time.h>

typedef double runtype;
typedef std::vector<runtype> Array;
typedef std::vector<Array> Matrix;
typedef std::vector<Matrix> array3d;


ofstream outfile1("loaded_die.txt");
ofstream outfile2("birthday.txt");
ofstream outfile3("needle.txt");
ofstream outfile4("two_dice.txt");
ofstream outfile5("shielding.txt");

void loaded_die(ofstream&);
runtype probably(int, runtype);
void needle(ofstream&);
void two_dice(ofstream&);
void shielding(ofstream&);

void main()
{
	srand( (unsigned)time( NULL ) );
	loaded_die(outfile1);

// Chapter 11.3, Birthday Probability
	runtype npts = 1000;
	outfile2<<"Birthday Problem Simulation with n persons"<<endl;
	outfile2<<setw(3)<<" t "<<setw(8)<<"simulation"<<endl;
	outfile2<<setw(3)<<"5" <<setw(8)<< probably(5, npts)<<endl;
	outfile2<<setw(3)<<"10"<<setw(8)<< probably(10, npts)<<endl;
	outfile2<<setw(3)<<"15"<<setw(8)<< probably(15, npts)<<endl;
	outfile2<<setw(3)<<"20"<<setw(8)<< probably(20, npts)<<endl;
	outfile2<<setw(3)<<"22"<<setw(8)<< probably(22, npts)<<endl;
	outfile2<<setw(3)<<"23"<<setw(8)<< probably(23, npts)<<endl;
	outfile2<<setw(3)<<"25"<<setw(8)<< probably(25, npts)<<endl;
	outfile2<<setw(3)<<"30"<<setw(8)<< probably(30, npts)<<endl;
	outfile2<<setw(3)<<"35"<<setw(8)<< probably(35, npts)<<endl;
	outfile2<<setw(3)<<"40"<<setw(8)<< probably(40, npts)<<endl;
	outfile2<<setw(3)<<"45"<<setw(8)<< probably(45, npts)<<endl;
	outfile2<<setw(3)<<"50"<<setw(8)<< probably(50, npts)<<endl;
	outfile2<<setw(3)<<"55"<<setw(8)<< probably(55, npts)<<endl;

	needle(outfile3);

	two_dice(outfile4);

	shielding(outfile5);
}




void loaded_die(ofstream& out)
{
// Chapter 11.3, Loaded Die
	const int n=5000;
	const int iprt=1000;

	Array y(7, 0);
	y[1] = 0.2;
	y[2] = 0.34;
	y[3] = 0.56;
	y[4] = 0.72;
	y[5] = 0.89;
	y[6] = 1.0 ;
	Array m(7, 0);
	m[1] = 0.0;
	m[2] = 0.0;
	m[3] = 0.0;
	m[4] = 0.0;
	m[5] = 0.0;
	m[6] = 0.0;
	Array r(n+1, 0);

	int i, j;

	for (j=1; j<= n; j++)
	{
		r[j] = (runtype) rand()/(runtype)RAND_MAX;
	}

	for (i = 1; i <= n; i++)
		for (j = 1; j <= 6; j++)
			if (r[i] <= y[j])
			{
				m[j] = m[j] + 1;
				j=7;
			}
	for (i = 1; i <= 6; i++)
		out<<"m["<< i << "] = " << m[i]<<endl;
}


bool birthday(int n)
{
	std::vector<double> r(n+1, 0);
	std::vector<bool> days(365+1, 0);
    int i, number;

	for (i=1; i<= n; i++)
		r[i] = (double) rand()/(double)RAND_MAX;
 
	for (i=1;i<= 365; i++)
		days[i] = false;
	for(i=1;i<=n;i++)
	{
		number = (int)(365*r[i] + 1);
		if (days[number])
			return true;
		days[number] = true;
	}
	return false;
}


runtype probably(int n, runtype npts)
{
	int i;
	runtype sum;

	sum = 0;
	for (i = 1; i<= npts; i++)
		if(birthday(n))
			sum++;
	return sum/npts;
}



void needle(ofstream& out)
{
// Chapter 11.3, Needle
	const int n=5000;
	const int iprt=1000;
	Matrix r(n+1, Array(2+1, 0) );
	runtype pi, prob, t, v, w;
	int m, i;

	m = 0;
	pi = 3.141592654;
	t = 2/pi;
	for (i = 1; i <= n; i++)
	{
		r[i][1] = (runtype)rand()/(runtype)RAND_MAX;
		r[i][2] = (runtype)rand()/(runtype)RAND_MAX;
	}
	for (i = 1; i <= n; i++)
	{
		w = r[i][1];
		v = (pi/2)*r[i][2];
		if (w <= sin(v))
			m++;
		if ( (i%iprt) == 0)
		{
			prob = (runtype)m/(runtype)i;
			out<<"i = "<<setw(5)<<i<<"  prob = "<<setw(8)<<prob<<"   2/pi ="<<setw(8)<<t<<endl;
		}      
	} 
}


void two_dice(ofstream& out)
{
// Chapter 11.3, Two Dice
	const int n=5000;
	const int iprt=1000;
	array3d r(n+1, Matrix(24+1, Array(2+1, 0)));
	int i, j, m;

  
	int i1, i2;
	runtype prob;

	for (i = 1; i <= n; i++)
		for (j = 1; j <= 24; j++)
		{
			r[i][j][1] = (runtype)rand()/RAND_MAX;
			r[i][j][2] = (runtype)rand()/RAND_MAX;
		}
	m = 0;
	for (i = 1; i <= n; i++)
	{
		j=1;
		while (j<=24)
		{
			i1 = (int)(6*r[i][j][1] + 1);
			i2 = (int)(6*r[i][j][2] + 1); 
			if ( (i1 + i2) == 12)
			{
				m++;
				j = 25;
			}
			j++;
		}
		if ( (i%1000) == 0)
		{
			prob = (runtype)m/(runtype)i;
			out<<"i = "<<setw(5)<<i<<"  prob ="<<setw(8)<<prob<<endl;
		}
	}
}


void shielding(ofstream& out)
{
// Chapter 11.3, Shielding
	const int n=5000;
	const int iprt=1000;

	Matrix r(n+1, Array(7+1,0));
	runtype x, per, pi;
	int i, j, m;

	m = 0;
	for (i = 1; i <= n; i++)
		for (j = 1; j <= 7; j++)
			r[i][j] = (runtype)rand() /(runtype)RAND_MAX;
	for (i=1; i<=n; i++)
	{
		x = 1;
		j=1;
		while (j<=7)
		{
			x += cos(pi*r[i][j]);
			if (x <= 0)
			j = n +20;//instead of break
			if (x >= 5)
			{
				m++;
				j= n+20;//instead of break
			}
			j++;
		}
		if ( (i%iprt) == 0)
		{
			per = 100*(runtype)m/(runtype)i;
			out<<"i = "<<setw(5)<<i<<"  per ="<<setw(8)<<per<<endl;
		}
	}
}