d13r9f.java

JAVA数值计算算法集随书源码

public class d13r9F
{
	double chsone_df;
	double chsone_chsq;
	double chsone_prob;
	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)
			{
                System.out.println("bad Math.expected number");
                System.exit(1);
            }
            chsq = chsq + (bins[j] - ebins[j]) * (bins[j] - ebins[j]) / ebins[j];
        }
        prob = gammq(0.5 * df, 0.5 * chsq);
        chsone_df = df;
        chsone_chsq = chsq;
        chsone_prob = prob;
	}

	double gammln(double xx)
	{
		int j;
		double temp;
		double cof[] = new double[7];
		double stp, half, one, fpf, x, tmp, ser;
		cof[1] = 76.18009173;
		cof[2] = -86.50532033;
		cof[3] = 24.01409822;
        cof[4] = -1.231739516;
        cof[5] = 0.00120858003;
        cof[6] = -0.00000536382;
        stp = 2.50662827465;
        half = 0.5;
        one = 1.0;
        fpf = 5.5;
        x = xx - one;
        tmp = x + fpf;
        tmp = (x + half) * Math.log(tmp) - tmp;
        ser = one;
        for (j = 1; j <= 6; j++)
		{
            x = x + one;
            ser = ser + cof[j] / x;
        }
        temp = tmp + Math.log(stp * ser);
		return temp;
	}

	double gammq_a, gammq_x;
	double gammq(double a, double x)
	{
		double temp, gamser, gammcf, gln;
		gamser = 0; gln = 0; gammcf = 0;
        if (x < 0.0 || a <= 0.0)
		{
            System.out.println("pause in gammq");
            System.exit(1);
        }
        if (x < a + 1.0)
		{
            gser(gamser, a, x, gln);
			gamser = gser_gamser; a  = gser_a; x = gser_x; gln = gser_gln;
            temp= 1.0 - gamser;
        }
        else
		{
            gcf(gammcf, a, x, gln);
			gammcf = gcf_gammcf; a = gcf_a; x = gcf_x; gln = gcf_gln;
            temp = gammcf;
		}
		gammq_a = a;
		gammq_x = x;
		return temp;   
	}

	double gcf_gammcf, gcf_a, gcf_x, gcf_gln;
	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 (Math.abs((g - gold) / g) < eps )
                {
                    gammcf = Math.exp(-x + a * Math.log(x) - gln) * g;
					break;
				}
                gold = g;
            }
        }
		gcf_gammcf = gammcf; gcf_a = a; gcf_x = x; gcf_gln = gln;
	}

	double gser_gamser, gser_a, gser_x, gser_gln;
	void gser(double gamser, double a, double x, double gln)
	{
		int itmax, n;
		double ap, sum, del, eps;
        itmax = 100;
        eps = 0.0000003;
        gln = gammln(a);
        if (x <= 0.0)
		{
            if (x < 0.0)
            {
                System.out.println("pause in gser");
                System.exit(1);
            }
            gamser = 0.0;
			gser_gamser = gamser; gser_a = a; gser_x = x; gser_gln = gln;
            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 (Math.abs(del) < (Math.abs(sum) * eps))
			{
                gamser = sum * Math.exp(-x + a * Math.log(x) - gln);
				break;
			}
        }
		gser_gamser = gamser; gser_a = a; gser_x = x; gser_gln = gln;
	}

	int expedv_idum;
	double expdev(int idum)
	{
		double temp;
        double dum;
		dum = 0.0;
        while (dum == 0)
		{
            dum = ran1(idum);
			idum = ran1_idum;
		}	
        temp = -Math.log(dum);
		expedv_idum = idum;
		return temp;
	}

	int ran1_idum;
	static double r[] = new double[98];
    static int ix1, ix2, ix3;
    double ran1(int idum)
	{
        double rm1, rm2 ,t;
		int m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j;
		boolean iff; iff=true;
        m1 = 259200; ia1 = 7141; ic1 = 54773; rm1 = 0.0000038580247;
        m2 = 134456; ia2 = 8121; ic2 = 28411; rm2 = 0.0000074373773;
        m3 = 243000; ia3 = 4561; ic3 = 51349;     
        if ((idum < 0) || (iff == false))
		{
            iff = true;
            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)) 
        {
            System.out.println( "abnormal exit");
            System.exit(1);
        }
        t = r[j];
        r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
        ran1_idum = idum;
        return t;
	}
}