#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>

double fnc(double &x1, double &x2, double &x3, double &x4)
{
	double temp;
    temp = sqrt(x1*x1 + x2*x2 + x3*x3 + x4*x4);
	return temp;
}

void main()
{
    //program d6r3
    //driver for routine ran2
    //calculates pi statistically using volume of unit n-sphere
    int  idum,i,j,m;
	double iy[4],yprob[4],x1,x2,x3,x4,a,b,k;
    const double pi = 3.1415926; 
	a=0.0;
	b=0.0;
    idum = -1;
    for (i = 1; i<=3; i++)
        iy[i] = 0;
    cout<<"Volume of unit n-sphere, n=2,3,4"<<endl;
    cout<<"# points     pi     (4/3)*pi    (1/2)*pi^2"<<endl;
    //randize
    for (j = 1; j<=15; j++)
	{
        for( k = pow(2 ,(j - 1)); k<=pow(2,j); k++)
		{
            x1 = ran2(idum);
            x2 = ran2(idum);
            x3 = ran2(idum);
            x4 = ran2(idum);
            if (fnc(x1, x2, a,b) < 1.0)
				iy[1] = iy[1] + 1;
            if (fnc(x1, x2, x3, b) < 1.0) 
				iy[2] = iy[2] + 1;
            if (fnc(x1, x2, x3, x4) < 1.0) 
				iy[3] = iy[3] + 1;
        }
        for (i = 1; i<=3; i++)
            yprob[i] = 1.0 * pow(2,(i + 1)) * iy[i] / (pow(2,j));
        cout<<setprecision(0)<<setiosflags(ios::fixed)<<setw(8)<<pow(2,j);
        for (m = 1; m<=3; m++)
		{
            cout<< setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<yprob[m];
		}
        cout<<endl;
     }
    cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(8)<<"  actual      "<<pi;
    cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<4.0 * pi / 3.0;
    cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<0.5 * pow(pi,2)<<endl;
}