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

double ran0(int idum)
{
    int iff,iseed,j;
	double v[98],y,t1,t2;
	iff =0;
    if (idum < 0 || iff ==0 )
	{
        iseed = int(fabs(idum));
        iff = 1;
        idum = 1;
        for (j = 1;j<=97;j++)
            v[j] = double(rand())/double(RAND_MAX);
        t1=rand(); t2=RAND_MAX;
        y =t1/t2;
    }
    j = 1 + int(97 * y);
    if (j > 97 || j < 1 )
	{
		return y;
	}
    y = v[j];
	return y;
    v[j] = double(rand())/double(RAND_MAX);
}

void main()
{
    //program d6r1
    //driver for routine ran0
    //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 = ran0(idum);
            x2 = ran0(idum);
            x3 = ran0(idum);
            x4 = ran0(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(10)<< pow(2,j);
        for (m = 1; m<=3; m++)
            cout<<setprecision(5)<<setiosflags(ios::fixed)<<setw(14)<<yprob[m];
		cout<<endl;
    }
	cout<<setprecision(5)<<setiosflags(ios::fixed)<<setw(8)<<"    actual       "<<pi;
    cout<<setprecision(5)<<setiosflags(ios::fixed)<<setw(14)<<4.0 * pi / 3.0;
    cout<<setprecision(5)<<setiosflags(ios::fixed)<<setw(14)<<0.5 * pow(pi,2)<<endl;
}