public static double betacf(double a, double b, double x)
        {
                double MAXIT = 100;
                double EPS = 3e-7;
                double FPMIN = 1e-30;
                int m, m2;
                double aa, c, d, del, h, qab, qam, qap;

                qab = a+b;
                qap = a+1;
                qam = a-1;
                c = 1;
                d = 1-qab*x/qap;
                if (Math.abs(d) < FPMIN) d = FPMIN;
                d = 1/d;
                h = d;
                for (m = 1; m <= MAXIT; m++)
                {
                        m2 = 2*m;
                        aa = m*(b-m)*x/((qam+m2)*(a+m2));
                        d = 1+aa*d;
                        if (Math.abs(d) < FPMIN) d = FPMIN;
                        c = 1+aa/c;
                        if (Math.abs(c) < FPMIN) c = FPMIN;
                        d = 1/d;
                        h *= d*c;
                        aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
                        d = 1+aa*d;
                        if (Math.abs(d) < FPMIN) d = FPMIN;
                        c = 1+aa/c;
                        if (Math.abs(c) < FPMIN) c = FPMIN;
                        d = 1/d;
                        del = d*c;
                        h *= del;
                        if (Math.abs(del-1) < EPS) break;
                }
                if (m > MAXIT)
                {
                        System.out.println("a or b too big, or MAXIT too small in betacf");
                        return Double.NaN;
                }
                else
                        return h;
        }