d9r1f.java

解具有导数的多元函数的无约束极小值点或具有导数的非线性方程组的近似解

public class d9r1F
{
	static int iset;
	static double gset;
	int gasdev_idum;
	double gasdev(int idum)
	{		
        double t, v1, v2, fac, r;
        if (iset == 0.0)
		{
            do
            {
				v1 = 2.0 * ran1(idum) - 1.0;
				idum = ran1_idum;
				v2 = 2.0 * ran1(idum) - 1.0;
				idum = ran1_idum;
                r = v1 * v1 + v2 * v2;
			}while ((r >= 1.0) || (r == 0.0));
			fac = Math.sqrt(-2.0 * Math.log(r) / r);
			gset = v1 * fac;
			t = v2 * fac;
			iset = 1;
			gasdev_idum = idum;
			return t;
		}
        else
		{
			t = gset;
			gasdev_idum = idum;
			iset = 0;
			return t;
        }
	}
	
	int ran1_idum;
	static double r[] = new double[98];
    static int ix1, ix2, ix3;
	static int iff;
	double ran1(int idum)
	{
        double rm1, rm2 ,t;
		int m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j;		 
        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 == 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)) 
        {
            System.out.println( "abnormal exit");
            System.exit(1);
        }
        t = r[j];
        r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
        ran1_idum = idum;
        return t;
	}
	double fit_a;
	double fit_b;
	double fit_siga;
	double fit_sigb;
	double fit_chi2;
	double fit_q;
	void fit(double x[], double y[], int ndata, double sig[], int mwt, double a,
    double b, double siga, double sigb, double chi2, double q)
	{
		int i;
		double sigdat, t, sxoss, wt, ss, sx = 0.0;
		double sy = 0.0;
		double st2 = 0.0;
		b = 0.0;
		if (mwt != 0)
		{
            ss = 0.0;
            for (i = 1; i <= ndata; i++)
            {
                wt = 1.0 / (sig[i] * sig[i]);
                ss = ss + wt;
                sx = sx + x[i] * wt;
                sy = sy + y[i] * wt;
            }
		}
		else
		{
            for (i = 1; i <= ndata; i++)
            {
                sx = sx + x[i];
                sy = sy + y[i];
            }
            ss = ndata;
		}
		sxoss = sx / ss;
		if (mwt != 0)
		{
            for (i = 1; i <= ndata; i++)
			{
                t = (x[i] - sxoss) / sig[i];
                st2 = st2 + t * t;
                b = b + t * y[i] / sig[i];
            }
		}
		else
		{
            for (i = 1; i <= ndata; i++)
			{
                t = x[i] - sxoss;
                st2 = st2 + t * t;
                b = b + t * y[i];
            }
		}
		b = b / st2;
		a = (sy - sx * b) / ss;
		siga = Math.sqrt((1.0 + sx * sx / (ss * st2)) / ss);
		sigb = Math.sqrt(1.0 / st2);
		chi2 = 0.0;
		if (mwt == 0)
		{
            for (i = 1; i <= ndata; i++)
			{
                chi2 = chi2 + Math.pow((y[i] - a - b * x[i]), 2);
            }
            q = 1.0;
            sigdat = Math.sqrt(chi2 / (ndata - 2));
            siga = siga * sigdat;
            sigb = sigb * sigdat;
		}
		else
		{
            for (i = 1; i <= ndata; i++)
			{
                chi2 = chi2 + Math.pow(((y[i] - a - b * x[i]) / sig[i]) , 2);
			}
            q = gammq(0.5 * (ndata - 2), 0.5 * chi2);
		}
		fit_a = a; fit_b = b; fit_siga = siga; fit_sigb = sigb; fit_chi2 = chi2; fit_q = q;
	}

	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;
	}

	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;
	}
}