public static double betacf( double a, double b, double x )
	{
		int MAXIT=100;
		double EPS = getEpsilon();
		double FPMIN = Double.MIN_VALUE / EPS;
		
		int m, m2;
		double aa, c, d, del, h, qab, qam, qap;
		
		qab = a+b;
		qap = a + 1.0;
		qam = a - 1.0;
		
		c=1.0;
		d=1.0-qab*x/qap;
		
		if ( Math.abs(d) < FPMIN ) d=FPMIN;
		d = 1.0/d;
		h=d;
		
		for( m=1;m<=MAXIT;m++) {
			m2 = 2*m;
			aa = m*(b-m)*x / ((qam+m2)*(a+m2));
			d = 1.0 + aa * d;
			if ( Math.abs(d) < FPMIN ) d = FPMIN;
			c= 1.0 + aa / c;
			if ( Math.abs(c) < FPMIN) c = FPMIN;
			d= 1.0 / d;
			h *= d * c;
			aa = -(a+m) * (qab + m ) * x / ((a+m2) * ( qap + m2));
			d=1.0 + aa * d;
			if ( Math.abs(d) < FPMIN) d = FPMIN;
			c = 1.0 + aa / c;
			if ( Math.abs(c) < FPMIN) c= FPMIN;
			d=1.0/d;
			del=d*c;
			h *= del;
			if ( Math.abs( del - 1.0 ) <= EPS ) break;
		}
		
		if ( m > MAXIT ) 
		{
		    
			throw new RuntimeException("a or b too big, of MAXIT too small");
		}
		
		return h;
	}