public static double incompleteGammaComplement(double a,double x){
  double ans, ax, c, yc, r, t, y, z;
  double pk, pkm1, pkm2, qk, qkm1, qkm2;
  if (x <= 0 || a <= 0) {
    return 1.0;
  }
  if (x < 1.0 || x < a) {
    return 1.0 - incompleteGamma(a,x);
  }
  ax=a * Math.log(x) - x - lnGamma(a);
  if (ax < -MAXLOG) {
    return 0.0;
  }
  ax=Math.exp(ax);
  y=1.0 - a;
  z=x + y + 1.0;
  c=0.0;
  pkm2=1.0;
  qkm2=x;
  pkm1=x + 1.0;
  qkm1=z * x;
  ans=pkm1 / qkm1;
  do {
    c+=1.0;
    y+=1.0;
    z+=2.0;
    yc=y * c;
    pk=pkm1 * z - pkm2 * yc;
    qk=qkm1 * z - qkm2 * yc;
    if (qk != 0) {
      r=pk / qk;
      t=Math.abs((ans - r) / r);
      ans=r;
    }
 else {
      t=1.0;
    }
    pkm2=pkm1;
    pkm1=pk;
    qkm2=qkm1;
    qkm1=qk;
    if (Math.abs(pk) > big) {
      pkm2*=biginv;
      pkm1*=biginv;
      qkm2*=biginv;
      qkm1*=biginv;
    }
  }
 while (t > MACHEP);
  return ans * ax;
}

public static double incompleteGammaComplement ( double a , double x ) { double ans , ax , c , yc , r , t , y , z ; double pk , pkm1 , pkm2 , qk , qkm1 , qkm2 ; if ( x <= 0 || a <= 0 ) { return 1.0 ; } if ( x < 1.0 || x < a ) { return 1.0 - incompleteGamma ( a , x ) ; } ax = a * Math . log ( x ) - x - lnGamma ( a ) ; if ( ax < - MAXLOG ) { return 0.0 ; } ax = Math . exp ( ax ) ; y = 1.0 - a ; z = x + y + 1.0 ; c = 0.0 ; pkm2 = 1.0 ; qkm2 = x ; pkm1 = x + 1.0 ; qkm1 = z * x ; ans = pkm1 / qkm1 ; do { c += 1.0 ; y += 1.0 ; z += 2.0 ; yc = y * c ; pk = pkm1 * z - pkm2 * yc ; qk = qkm1 * z - qkm2 * yc ; if ( qk != 0 ) { r = pk / qk ; t = Math . abs ( ( ans - r ) / r ) ; ans = r ; } else { t = 1.0 ; } pkm2 = pkm1 ; pkm1 = pk ; qkm2 = qkm1 ; qkm1 = qk ; if ( Math . abs ( pk ) > big ) { pkm2 *= biginv ; pkm1 *= biginv ; qkm2 *= biginv ; qkm1 *= biginv ; } } while ( t > MACHEP ) ; return ans * ax ; }
returns the complemented incomplete gamma function .
