d13r11.cpp

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

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

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 sort(int n, double ra[])
{
	int i,j,l,ir;
	double rra;
	l = n / 2 + 1;
	ir = n;
	do
	{
		if (l > 1)
		{
			l = l - 1;
			rra = ra[l];
		}
		else
		{
			rra = ra[ir];
			ra[ir] = ra[1];
			ir = ir - 1;
			if (ir == 1)
			{
				ra[1] = rra;
				return;
			}
		}
		i = l;
		j = l + l;
		while (j <= ir)
		{		
			if (j < ir)
			{
				if (ra[j] < ra[j + 1])
				{
					j = j + 1;
				}
			}
			if (rra < ra[j])
			{
				ra[i] = ra[j];
				i = j;
				j = j + j;
			}
			else
			{
				j = ir + 1;
			}
		}
		ra[i] = rra;
	}while(1);
}

double probks(double alam)
{
	double eps1 = 0.001;
    double eps2 = 0.00000001;
    double a2 = -2.0 * alam * alam;
    double fac = 2.0;
    double aaa,temp = 0.0;
    double term,termbf = 0.0;
    for (int j = 1; j<=100; j++)
	{
        term = fac * exp(a2 * j * j);
        temp = temp + term;
        aaa = fabs(term);
        if (aaa < eps1 * termbf || aaa < eps2 * temp)
		{
			return temp;
        }
        fac = -fac;
        termbf = fabs(term);
    }
    return 1.0;
}

void kstwo(double data1[], int n1, double data2[], int n2, double& d, double& prob)
{
    sort(n1, data1);
    sort(n2, data2);
    int en1 = n1;
    int en2 = n2;
    int j1 = 1;
    int j2 = 1;
    double am,fn1,fo1 = 0.0;
    double dt,fn2,fo2 = 0.0;
    d = 0.0;
    while(j1 <= n1 && j2 <= n2)
	{
        if (data1[j1] < data2[j2])
		{
            fn1 = (double)j1 / (double)en1;
            if (fabs(fn1 - fo2) > fabs(fo1 - fo2))
			{
                am = fabs(fn1 - fo2);
			}
            else
			{
                am = fabs(fo1 - fo2);
            }
            dt = am;
            if (dt > d)
			{
				d = dt;
			}
            fo1 = fn1;
            j1 = j1 + 1;
		}
        else
		{
            fn2 = (double)j2 / (double)en2;
            if (fabs(fn2 - fo1) > fabs(fo2 - fo1))
			{
                am = fabs(fn2 - fo1);
			}
            else
			{
                am = fabs(fo2 - fo1);
            }
            dt = am;
            if (dt > d)
			{
				d = dt;
			}
            fo2 = fn2;
            j2 = j2 + 1;
        }
    }
    prob = probks(sqrt((double)en1 * (double)en2 / (en1 + en2)) * d);
}

void main()
{
    //program d13r11
    //driver for routine kstwo
    int j,n1 = 2000;
    int n2 = 1000;
    double d,prob,var,factr,eps = 0.1;
    double data1[2001], data2[1001];
    long idum = -1357;
    for (j = 1; j<=n1; j++)
	{
        data1[j] = gasdev(idum);
    }
    cout<<endl;
    cout<<"Variance ratio  K-S statistic      Probability"<<endl;
	cout.setf(ios::fixed|ios::left);
	cout.precision(4);
    idum = -2468;
    for (int i = 1; i<=11; i++)
	{
        var = 1.0 + (i - 1) * eps;
        factr = sqrt(var);
        for (j = 1; j<=n2; j++)
		{
            data2[j] = factr * gasdev(idum);
        }
        kstwo(data1, n1, data2, n2, d, prob);
		cout<<"   ";
		cout.width(16);
        cout<<var;
		cout.width(16);
        cout<<d;
		cout.unsetf(ios::fixed);
		cout.setf(ios::scientific);
        cout<<prob<<endl;
		cout.unsetf(ios::scientific);
		cout.setf(ios::fixed);
    }
}