double betacf(double a,double b,double x)
{
	double itmax,eps,am,bm,az,qab,qap,qam,bz,em,tem,ap,bp;
	double aap,bpp,aold,temp,d;
	int m;
    itmax = 100;
    eps = 0.0000003;
    am = 1.0;
    bm = 1.0;
    az = 1.0;
    qab = a + b;
    qap = a + 1.0;
    qam = a - 1.0;
    bz = 1.0 - qab * x / qap;
    for( m = 1;m<=itmax;m++)
	{
        em = m;
        tem = em + em;
        d = em * (b - m) * x / ((qam + tem) * (a + tem));
        ap = az + d * am;
        bp = bz + d * bm;
        d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem));
        aap = ap + d * az;
        bpp = bp + d * bz;
        aold = az;
        am = ap / bpp;
        bm = bp / bpp;
        az = aap / bpp;
        bz = 1.0;
        if( fabs(az - aold) < eps * fabs(az)) 
			goto yi;
    }
    cout<<"a or b too big, or itmax too small";
yi: temp = az;
	return temp;
}