#include <math.h>
double gamma(double x);

int incog(double a,double x,double &gin,double &gim,double &gip)
{
    double xam,r,s,ga,t0;
    int k;

    if ((a < 0.0) || (x < 0)) return 1;
    xam = -x+a*log(x);
    if ((xam > 700) || ( a > 170.0)) return 1;
    if (x == 0.0) {
        gin = 0.0;
        gim = gamma(a);
        gip = 0.0;
        return 0;
    }
    if (x <= 1.0+a) {
        s = 1.0/a;
        r = s;
         for (k=1;k<=60;k++) {
            r *= x/(a+k);
            s += r;
            if (fabs(r/s) < 1e-15) break;
        }
        gin = exp(xam)*s;
        ga = gamma(a);
        gip = gin/ga;
        gim = ga-gin;
    }
    else {
        t0 = 0.0;
        for (k=60;k>=1;k--) {
            t0 = (k-a)/(1.0+k/(x+t0));
        }
        gim = exp(xam)/(x+t0);
        ga = gamma(a);
        gin = ga-gim;
        gip = 1.0-gim/ga;
    }
    return 0;
}