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);
 
    /* continued fraction */
    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;
  }