d6r3.cpp

这个是常用的数值算法用VC编写的。相信对大家有用哦。

#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 ran2(int& idum)
{
	static int iy,iff,ir[97];
	double rm,ran;
    int j,m,ic,ia;
    m = 714025; ia = 1366; ic = 150889; 
    if (idum < 0 || iff == 0) 
	{
        iff = 1;
        idum = (ic - idum) % m;
        for( j = 1;j<=97;j++)
		{
            idum = (ia*idum + ic) % m;
            ir[j] = idum;
        }

        idum = (ia * idum + ic) % m;
        iy = idum;
    }
    j = 1 + int((97 * iy) / m);
   if(j > 97 || j < 1)
   {
	   cout<<"abnormal exit"; 
	   exit(1);
   }
    rm = 0.0000014005112;
    iy = ir[j];
    ran = iy * rm;
    idum = (ia * idum + ic) % m;
    ir[j] = idum;
	return ran;
}

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