private static double gcf(double a,double x,double gln){
  double gold=0.0, g, fac=1.0, b1=1.0;
  double b0=0.0, anf, ana, an, a1, a0=1.0;
  double p;
  boolean done=false;
  a1=x;
  p=0.0;
  for (an=1.0; (!done) && (an <= ITMAX); an+=1.0) {
    ana=an - a;
    a0=(a1 + a0 * ana) * fac;
    b0=(b1 + b0 * ana) * fac;
    anf=an * fac;
    a1=x * a0 + anf * a1;
    b1=x * b0 + anf * b1;
    if (a1 != 0.0) {
      fac=1.0 / a1;
      g=b1 * fac;
      if (Math.abs((g - gold) / g) < EPSILON) {
        p=Math.exp(-x + a * Math.log(x) - gln) * g;
        done=true;
      }
      gold=g;
    }
  }
  return p;
}

private static double gcf ( double a , double x , double gln ) { double gold = 0.0 , g , fac = 1.0 , b1 = 1.0 ; double b0 = 0.0 , anf , ana , an , a1 , a0 = 1.0 ; double p ; boolean done = false ; a1 = x ; p = 0.0 ; for ( an = 1.0 ; ( ! done ) && ( an <= ITMAX ) ; an += 1.0 ) { ana = an - a ; a0 = ( a1 + a0 * ana ) * fac ; b0 = ( b1 + b0 * ana ) * fac ; anf = an * fac ; a1 = x * a0 + anf * a1 ; b1 = x * b0 + anf * b1 ; if ( a1 != 0.0 ) { fac = 1.0 / a1 ; g = b1 * fac ; if ( Math . abs ( ( g - gold ) / g ) < EPSILON ) { p = Math . exp ( - x + a * Math . log ( x ) - gln ) * g ; done = true ; } gold = g ; } } return p ; }
compute gamma cdf by its series representation
