public static double incompleteBetaFraction1(double a,double b,double x){
  double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
  double k1, k2, k3, k4, k5, k6, k7, k8;
  double r, t, ans, thresh;
  int n;
  k1=a;
  k2=a + b;
  k3=a;
  k4=a + 1.0;
  k5=1.0;
  k6=b - 1.0;
  k7=k4;
  k8=a + 2.0;
  pkm2=0.0;
  qkm2=1.0;
  pkm1=1.0;
  qkm1=1.0;
  ans=1.0;
  r=1.0;
  n=0;
  thresh=3.0 * MACHEP;
  do {
    xk=-(x * k1 * k2) / (k3 * k4);
    pk=pkm1 + pkm2 * xk;
    qk=qkm1 + qkm2 * xk;
    pkm2=pkm1;
    pkm1=pk;
    qkm2=qkm1;
    qkm1=qk;
    xk=(x * k5 * k6) / (k7 * k8);
    pk=pkm1 + pkm2 * xk;
    qk=qkm1 + qkm2 * xk;
    pkm2=pkm1;
    pkm1=pk;
    qkm2=qkm1;
    qkm1=qk;
    if (qk != 0) {
      r=pk / qk;
    }
    if (r != 0) {
      t=Math.abs((ans - r) / r);
      ans=r;
    }
 else {
      t=1.0;
    }
    if (t < thresh) {
      return ans;
    }
    k1+=1.0;
    k2+=1.0;
    k3+=2.0;
    k4+=2.0;
    k5+=1.0;
    k6-=1.0;
    k7+=2.0;
    k8+=2.0;
    if ((Math.abs(qk) + Math.abs(pk)) > big) {
      pkm2*=biginv;
      pkm1*=biginv;
      qkm2*=biginv;
      qkm1*=biginv;
    }
    if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
      pkm2*=big;
      pkm1*=big;
      qkm2*=big;
      qkm1*=big;
    }
  }
 while (++n < 300);
  return ans;
}

public static double incompleteBetaFraction1 ( double a , double b , double x ) { double xk , pk , pkm1 , pkm2 , qk , qkm1 , qkm2 ; double k1 , k2 , k3 , k4 , k5 , k6 , k7 , k8 ; double r , t , ans , thresh ; int n ; k1 = a ; k2 = a + b ; k3 = a ; k4 = a + 1.0 ; k5 = 1.0 ; k6 = b - 1.0 ; k7 = k4 ; k8 = a + 2.0 ; pkm2 = 0.0 ; qkm2 = 1.0 ; pkm1 = 1.0 ; qkm1 = 1.0 ; ans = 1.0 ; r = 1.0 ; n = 0 ; thresh = 3.0 * MACHEP ; do { xk = - ( x * k1 * k2 ) / ( k3 * k4 ) ; pk = pkm1 + pkm2 * xk ; qk = qkm1 + qkm2 * xk ; pkm2 = pkm1 ; pkm1 = pk ; qkm2 = qkm1 ; qkm1 = qk ; xk = ( x * k5 * k6 ) / ( k7 * k8 ) ; pk = pkm1 + pkm2 * xk ; qk = qkm1 + qkm2 * xk ; pkm2 = pkm1 ; pkm1 = pk ; qkm2 = qkm1 ; qkm1 = qk ; if ( qk != 0 ) { r = pk / qk ; } if ( r != 0 ) { t = Math . abs ( ( ans - r ) / r ) ; ans = r ; } else { t = 1.0 ; } if ( t < thresh ) { return ans ; } k1 += 1.0 ; k2 += 1.0 ; k3 += 2.0 ; k4 += 2.0 ; k5 += 1.0 ; k6 -= 1.0 ; k7 += 2.0 ; k8 += 2.0 ; if ( ( Math . abs ( qk ) + Math . abs ( pk ) ) > big ) { pkm2 *= biginv ; pkm1 *= biginv ; qkm2 *= biginv ; qkm1 *= biginv ; } if ( ( Math . abs ( qk ) < biginv ) || ( Math . abs ( pk ) < biginv ) ) { pkm2 *= big ; pkm1 *= big ; qkm2 *= big ; qkm1 *= big ; } } while ( ++ n < 300 ) ; return ans ; }
continued fraction expansion # 1 for incomplete beta integral .
