d6r2.cpp

C语言的编程合集！是我自己编的！很有用的！大家可以下载！有意见麻烦大家提出来以便本人及时修改！

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

static double r[98];
static int iff,ix1,ix2,ix3;
double ran1(int& idum)
{
    double rm1,rm2,t;
	int m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3,j;
	 
    m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 0.0000038580247,
    m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 0.0000074373773,
    m3 = 243000, ia3= 4561, ic3 = 51349;
	
    if ((idum < 0) || (iff == 0))
	{
        iff = 1;
        ix1 =  (ic1 - idum) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1)% m1;
            ix2 = (ia2 * ix2 + ic2)% m2;
            r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
        }
        idum = 1;
	}
      ix1 = (ia1 * ix1 + ic1) % m1;
      ix2 = (ia2 * ix2 + ic2) % m2;
      ix3 = (ia3 * ix3 + ic3) % m3;
      j = 1 + int((97 * ix3) / m3);
      if ((j > 97) || (j < 1)) 
	  {
		  cout<< "abnormal exit"<<endl;
	      return 1;
	  }
      t = r[j];
      r[j] = (float(ix1) + float(ix2) * rm2) * rm1;
	  return t;
}

void main()
{
    //program d6r2
    //driver for routine ran1
    //calculates pi statistically using volume of unit n-sphere
    double  yprob[4],a,b,pi,x1,x2,k,x3,x4;
	int iy[4],idum,j,i,m;
    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 = pow(2,(j - 1)); k<= pow(2,j); k++)
		{
            x1 = ran1(idum); 
            x2 = ran1(idum);
            x3 = ran1(idum);
            x4 = ran1(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<<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;
}