d13r8.cpp

C语言的编程合集！是我自己编的！很有用的！大家可以下载！有意见麻烦大家提出来以便本人及时修改！

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

double gammln(double xx)
{
	double cof[7];
	cof[1] = 76.18009173;
	cof[2] = -86.50532033;
	cof[3] = 24.01409822;
	cof[4] = -1.231739516;
	cof[5] = 0.00120858003;
	cof[6] = -0.00000536382;
	double ser,x,tmp,stp = 2.50662827465;
	double half = 0.5;
	double one = 1.0;
	double fpf = 5.5;
	int j;
	x = xx - one;
	tmp = x + fpf;
	tmp = (x + half) * log(tmp) - tmp;
	ser = one;
	for (j = 1; j<=6; j++)
	{
		x = x + one;
		ser = ser + cof[j] / x;
	}
	return tmp + log(stp * ser);
}
void gcf(double& gammcf, double a, double x, double gln)
{
	int an,n,itmax = 100;
	double anf,g,ana,eps = 0.0000003;
	gln = gammln(a);
	double gold = 0.0;
	double a1,a0 = 1.0;
	a1 = x;
	double b0 = 0.0;
	double b1 = 1.0;
	double fac = 1.0;
	for (n = 1; n<=itmax; n++)
	{
		an = n;
		ana = an - a;
		a0 = (a1 + a0 * ana) * fac;
		b0 = (b1 + b0 * ana) * fac;
		anf = an * fac;
		a1 = x * a0 + anf * a1;
		b1 = x * b0 + anf * b1;
		if (a1 != 0.0)
		{
			fac = 1.0 / a1;
			g = b1 * fac;
			if (fabs((g - gold) / g) < eps)
			{
				gammcf = exp(-x + a * log(x) - gln) * g;
				return;
			}
			gold = g;
		}
	}
	cout<< "a too large, itmax too small"<<endl;
	gammcf = exp(-x + a * log(x) - gln) * g;
}
void gser(double& gamser, double a, double x, double gln)
{
	int n,itmax = 100;
	double ap,sum,del,eps = 0.0000003;
	gln = gammln(a);
	if (x <= 0.0)
	{
		if (x < 0.0)
		{
			cout<<"pause in gser"<<endl;
			return;
		}
		gamser = 0.0;
		return;
	}
	ap = a;
	sum = 1.0 / a;
	del = sum;
	for (n = 1; n<=itmax; n++)
	{
		ap = ap + 1.0;
		del = del * x / ap;
		sum = sum + del;
		if (fabs(del) < fabs(sum) * eps)
		{
			gamser = sum * exp(-x + a * log(x) - gln);
			return;
		}
	}
	cout<< "a too large, itmax too small";
	gamser = sum * exp(-x + a * log(x) - gln);
}
double gammq(double a, double x)
{
	double gamser,gammcf,gln=0;
	if (x < 0.0 || a <= 0.0)
	{
		cout<< "pause"<<endl;
		exit(1);
	}
	if (x < a + 1.0)
	{
		gser(gamser, a, x, gln);
		return 1.0 - gamser;
	}
	else
	{
		gcf(gammcf, a, x, gln);
		return gammcf;
	}
}
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 expdev(long& idum)
{
    double dum = 0.0;
    while (dum == 0.0)
	{
        dum = ran1(idum);
    }
    return -log(dum);
}

void chsone(double bins[], double ebins[], int nbins, int knstrn, double& df, double& chsq, double& prob)
{
    df = double(nbins - 1 - knstrn);
    chsq = 0.0;
    for (int j = 1; j<=nbins; j++)
	{
        if (ebins[j] <= 0.0)
		{
            cout<<" bad expected number"<<endl;
            return;
        }
        chsq = chsq + (bins[j] - ebins[j]) * (bins[j] - ebins[j]) / ebins[j];
    }
    prob = gammq(0.5 * df, 0.5 * chsq);
}

void main()
{
    //program d13r8
    //driver for routine chsone
    int i,j,npts = 2000;
    int nbins = 10;
    double bins[11], ebins[11];
    long idum = -15;
    for (j = 1; j<=nbins; j++)
	{
        bins[j] = 0.0;
    }
	double x;
	int ibin;
    for (i = 1; i<=npts; i++)
	{
        x = expdev(idum);
        ibin = int(x * nbins / 3.0) + 1;
        if (ibin <= nbins)
		{
			bins[ibin] = bins[ibin] + 1.0;
		}
    }
    for (i = 1; i<=nbins; i++)
	{
        ebins[i] = 3.0 * npts / nbins * exp(-3.0 * (i - 0.5) / nbins);
    }
	double df,chsq,prob;
    chsone(bins, ebins, nbins, -1, df, chsq, prob);
    cout<<endl;
    cout<<"     Expected    Observed"<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(2);
    for (i = 1; i<=nbins; i++)
	{
		cout.width(12);
        cout<<ebins[i];
		cout.width(12);
        cout<<bins[i]<<endl;
    }
    cout<<endl;
	cout.unsetf(ios::fixed);
	cout.setf(ios::scientific);
	cout.precision(5);
    cout<<"Chi-squared: "<<chsq<<endl;
    cout<<"Probability: "<<prob<<endl;
}