gcf.cpp

代码

void gcf(double& gammcf, double& a, double& x, double& gln)
{
	int itmax,n;
	double eps,a0,a1,b0,b1,fac,an,ana,anf,gold,g;
    itmax = 100;
    eps = 0.0000003;
    gln = gammln(a);
    gold = 0.0;
    a0 = 1.0;
    a1 = x;
    b0 = 0.0;
    b1 = 1.0;
    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;
				break;
			}
            gold = g;
        }
    }
}