d13r5.cpp

使用VC++编写的大量数学算法的源代码

#include "iostream.h"
#include "math.h"
#include "stdlib.h"

void avevar(double data[], int n, double& ave, double& var)
{
    ave = 0.0;
    var = 0.0;
	int j;
	double s;
    for (j = 1; j<=n; j++)
	{
        ave = ave + data[j];
    }
    ave = ave / n;
    for (j = 1; j<=n; j++)
	{
        s = data[j] - ave;
        var = var + s * s;
    }
    var = var / (n - 1);
}

double ran1(long& idum)
{
    int j,iff=-1;
	static long ix1,ix2,ix3;
	static double r[98];
    long m1 = 259200; long m2 = 134456; long m3 = 243000;
	long ia1 = 7141; long ia2 = 8121; long ia3 = 4561;
	long ic1 = 54773; long ic2 = 28411; long ic3 = 51349;
	double rm1 = 0.0000038580247; double rm2 = 0.0000074373773;    
    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 in ran1"<<endl;
		exit(1);
	}
	double temp=r[j];
    r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
	return temp;
}

double gasdev(long& idum)
{
    static int iset;
	static double gset;
	double v1,v2,r,fac;
    if (iset == 0)
	{
		do
		{
			v1 = 2.0 * ran1(idum) - 1.0;
			v2 = 2.0 * ran1(idum) - 1.0;
			r = v1 * v1 + v2 * v2;
		}while (r >= 1.0 || r == 0);
		fac = sqrt(-2.0 * log(r) / r);
		gset = v1 * fac;
		iset = 1;
		return v2 * fac;
	}
    else
	{
		iset = 0;
		return gset;
    }
}

void main()
{
     //program d13r5
     //driver for routine avevar
     int npts = 1000;
     double ave,var,shift,eps = 0.1;
     double data[1001];
     //generate gaussian distributed data
     long idum = -5;
     cout<<endl;
     cout<<"     Shift    Average   Variance"<<endl;
	 cout.setf(ios::fixed|ios::right);
	 cout.precision(2);
     for (int i = 1; i<=11; i++)
	 {
         shift = (i - 1) * eps;
         for (int j = 1; j<=npts; j++)
		 {
             data[j] = shift + i * gasdev(idum);
         }
         avevar(data, npts, ave, var);
		 cout.width(10);
         cout<<shift;
		 cout.width(10);
         cout<<ave;
		 cout.width(10);
         cout<<var<<endl;
	 }
}