double betai(double a,double b,double x)
{
	double bt,aaa,temp;
    if (x < 0.0 || x > 1.0) cout<<"bad argument x in betai";
    if (x == 0.0 || x == 1.0)
	{
        bt = 0.0;
	}
    else
	{
        aaa = gammln(a + b) - gammln(a) - gammln(b);
        bt = exp(aaa + a * log(x) + b * log(1.0 - x));
	}
    if (x < ((a + 1.0) / (a + b + 2.0)))
	{
        temp = bt * betacf(a, b, x) / a;
        _c_exit();
	}
    else
	{
        temp = 1.0 - bt * betacf(b, a, 1.0 - x) / b;
        _c_exit();
	}
	return temp;
}