#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;

ofstream outfile1("coarse_check.txt");
ofstream outfile2("ellipse.txt");
ofstream outfile3("ellipse_err.txt");
ofstream outfile4("circle_err.txt");

void coarse_check(ofstream&);
void ellipse(ofstream&);
void ellipse_err(ofstream&);
void circle_err(ofstream&);

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

	coarse_check(outfile1);

	ellipse(outfile2);

	ellipse_err(outfile3);

	circle_err(outfile4);

}



void coarse_check(ofstream& out)
{
// Chapter 11.1, Coarse Check
	const int n = 10000;
	int i, j, m;
	runtype per;
	Array r(n+1, 0);

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

	for (i = 1; i <= n; i++)
	{
		if (r[i] <= 0.5)
			m++;
		if ( (i%1000) == 0)
		{
			per = 100*(runtype)m/(runtype)n;
			out << setw(8) <<"t = " << i << setw(15) << "per = "<< per <<endl; 
		}
	}
}

void ellipse(ofstream& out)
{
// Chapter 11.1, Ellipse
	const int n = 1000;
	const int npts = 2000;
	Array x(n+1, 0);
	Array  y(n+1, 0);
	Matrix   r(2+1, Array(npts+1, 0));
	runtype u, v, w;
	int i, j;

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

	j = 1;
	i = 1;
	while(i <= npts)
	{
		u = 4*r[1][i] - 2;
		v = 2*r[2][i] - 1;
		w = u*u + 4*v*v;
		if (w <= 4)
		{
			if (j == n)
				i = npts+20;
			x[j] = u;
			y[j] = v;
			out<< "x["<<setw(5)<< j <<"] = "<<setw(15)<<x[j]<<"   y["<<setw(5)<< j <<"] = "<<setw(15)<<y[j]<<endl;
			j ++;
		}
		i++;
	}
}


void ellipse_err(ofstream& out)
{
// Chapter 11.1, Ellipse Erroneous
	const int n = 1000;
	Array x(n+1, 0);
	Array y(n+1, 0);
	Matrix  r(2+1, Array(n+1, 0));
	int i, j;

	for (i = 1; i <= 2; i++)
	{
		for (j=1; j<= n; j++)
		{
			r[i][j] = (runtype) rand()/(runtype)RAND_MAX;
		}
	}
	for (i = 1; i <= n; i++)
	{
		x[i] = 4*r[1][i] - 2;
		y[i] = ((2*r[2][i] - 1)/2)*sqrt(4 - x[i]*x[i]);
		out<< "x["<<setw(4)<< i <<"] = "<<setw(15)<<x[i]<<"   y["<<setw(4)<< i <<"] = "<<setw(15)<<y[i]<<endl;
	}
}


void circle_err(ofstream& out)
{
// Chapter 11.1, Circle Erroneous
	const int n = 1000;
	Array x(n+1, 0);
	Array y(n+1, 0);
	Matrix  r(2+1, Array(n+1, 0));
	runtype pi;
	int i, j;

	pi = 3.141592654;

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

	for (i = 1; i <= n; i++)
	{
		x[i] = r[1][i]*cos(2*pi*r[2][i]);
		y[i] = r[1][i]*sin(2*pi*r[2][i]);
		out<< "x["<<setw(4)<< i <<"] = "<<setw(15)<<x[i]<<"   y["<< setw(4)<<i <<"] = "<<setw(15)<<y[i]<<endl;
	}  
}