d6r4.cpp

Visual C++ 常用数值算法集 源代码

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

static int inextp,iff,inext;
static double ma[57];
double  ran3(int& idum)
{
	double mk,b,fac,mj;
	int i,ii,k,mbig,mz;
	long mseed;
    mbig = 1000000000; mseed = 161803398; mz = 0; fac = 0.000000001;
    if ((idum < 0) || (iff == 0))
	{
        iff = 1;
        mj = mseed - fabs(idum);
        mj = mj - mbig * int(mj / mbig);
        ma[56] = mj;
        mk = 1;
        for (i = 1; i<=54; i++)
		{
            ii = (21 * i) - 55 * int((21.0 * i)/ 55.0);
            ma[ii] = mk;
            mk = mj - mk;
            if (mk < mz) mk = mk + mbig;
            mj = ma[ii];
        }
        for (k = 1; k<=4; k++) 
		{
            for (i = 1; i<=55; i++)
			{
                ma[i] = ma[i] - ma[1 + i + 30 - 55 * int((i + 30) / 55)];
                if (ma[i] < mz) ma[i] = ma[i] + mbig;
			}
        }
        inext = 0;
        inextp = 31;
        idum = 1;
    }
    inext = inext + 1;
    if (inext == 56 ) inext = 1;
    inextp = inextp + 1;
    if (inextp == 56)  inextp = 1;
    mj = ma[inext] - ma[inextp];
    if (mj < mz ) mj = mj + mbig;
    ma[inext] = mj;
    b = mj * fac;
	return b;
}

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