# include<math.h>
# include<iomanip.h>
# include<iostream.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 d6r4
    //driver for routine ran3
    //calculates pi statistically using volume of unit n-sphere
    double  yprob[4],a,b,pi,x1,x2,x3,x4;
	int iy[4],idum,j,i,m,k;
    idum = -1;
	pi = 3.1415926;
	a=0.0;
	b=0.0;
    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;
    for (j = 1; j<=15; j++)
	{
        for (k =int( pow(2,(j - 1))); int(k<= pow(2,j)); k++)
		{
            x1=ran3(idum); 
            x2=ran3(idum);
            x3=ran3(idum);
            x4=ran3(idum);
            if (fnc(x1, x2, a, b) < 1.0)  iy[1] = iy[1] + 1;
            if (fnc(x1, x2, x3, a) < 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<< setw(10)<<setprecision(0)<<setiosflags(ios::fixed)<<pow(2 , j);
        for (m = 1; m<=3; m++)
            cout<< setprecision(6)<<setiosflags(ios::fixed)<<setw(13)<<yprob[m];
		cout<<endl;
     }
    cout<< "    actual"<<setprecision(6)<<setiosflags(ios::fixed)<<setw(13)<<pi;
    cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(13)<<4.0 * pi / 3.0;
    cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(13)<<0.5*pow(pi,2)<<endl;
}