#include <math.h>
#include <fstream.h>
#include <iomanip.h>
#include <vector>
#include <stdlib.h>
#include <time.h>

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

inline runtype sqr(runtype x) { return x*x; }
inline runtype f(runtype x, runtype y) { return sin(sqrt(log(x + y + 1))); }

void double_integral(ofstream&);
void volume_region(ofstream&);
void cone(ofstream&);

ofstream outfile1("double_integral.txt");
ofstream outfile2("volume_region.txt");
ofstream outfile3("cone.txt");

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

	double_integral(outfile1);
	volume_region(outfile2);
	cone(outfile3);
}


void double_integral(ofstream& out)
{
// Chapter 11.2, Double Integral
	const int n = 5000, iprt = 1000;
	Matrix r(2+1, Array(n+1, 0));
	int i, j;
	runtype sum, vol, x, y;
 
	for (i = 1; i <= 2; i++)
	{
		for (j=1; j<= n; j++)
		{
			r[i][j] = (runtype) rand()/(runtype)RAND_MAX;
		}
	}
	j = 0;
	sum = 0;
	for (i = 1; i <= n; i++)
	{
		x = r[1][i];
		y = r[2][i];
		if ( sqr(x - 0.5) + sqr(y - 0.5) <= 0.25 )
		{
			j++;
			sum += f(x, y);
			if ( (j%iprt) == 0)
			{
				vol = atan(1.0)*sum/(runtype)j;
				out<<"j = "<<setw(5)<< j <<"    vol =" << vol<< endl;
			}
		}
	}
	vol = atan(1.0)*sum/(runtype)j;
	out << "j = "<<setw(5)<< j << "    vol =" <<vol <<endl;
}


void volume_region(ofstream& out)
{
// Chapter 11.2, Volume of a Region
	const int n=5000;
	const int iprt=1000;
	Matrix r(n+1, Array(3+1, 0));
	runtype vol, x, y, z;
	int i, j, m;

	for (i = 1; i <= n; i++)
		for (j = 1; j <= 3; j++)   
			r[i][j] = (double)rand()/(double)RAND_MAX;    
	m = 0;
	vol = 0;
	for (i = 1; i <= n; i++)
	{
		x = r[i][1];
		y = r[i][2];
		z = r[i][3];
		if ( ( x*x + sin(y) <= z) && (x - z + exp(y) <= 1))
			m++;
		if ( (i%iprt) == 0)
		{
			vol = (double)m/(double)i;
			out<<"i = "<<setw(5)<<i<<"  vol ="<<setw(8)<<vol<<endl;
		}
	}
}


void cone(ofstream& out)
{
// Chapter 11.2, Cone
	const int n = 5000, iprt = 1000;
	Matrix r(n+1, Array(3+1, 0) );
	int i, j, m;
	runtype vol, x, y, z;

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

	for (i = 1; i <= n; i++)
	{
		x = 2*r[i][1] - 1;
		y = 2*r[i][2] - 1;
		z = 2*r[i][3];
		if ( (x*x + y*y <= z*z) && (x*x + y*y <= z*(2 - z)))
			m++;
		if ( (i%iprt) == 0)
		{
			vol = 8*(runtype)m/(runtype)i;
			out<<"i = "<<setw(5)<<i<<"   vol = "<<setw(8)<<vol<<endl;
		}        	   
	}
}