question: 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 ; }
target:  Continued fraction expansion #1 for incomplete beta integral.