# Gamma(a,x) is incomplete Gamma(a);
# gamma(n,a,x) is the nth order derivative of gamma(a,x);
# gamma(1/2,x_) := -sqrt(pi)*erf(sqrt(x));
# gamma(1/2,-x_) := -sqrt(pi)*i*erfi(sqrt(x));
# gamma(0,0) := inf;
# gamma(0,x_) = -Ei(-x)=En(1,x);


gamma(0,a_,x_):=gamma(a,x);
gamma(1,a_,y_):= d(gamma(a,y),x);
gamma(-1,a_,x):=x*gamma(a,x)-gamma(a+1,x);

gamma(a_,0):= gamma(a);
gamma(a_,inf):= 0;
gamma(0,x_) := -Ei(-x);
gamma(1,x_):= exp(-x);

gamma(-1/2):= -2*sqrt(pi);
gamma(1/2):= sqrt(pi);
gamma(3/2):= sqrt(pi)/2;
gamma(5/2):= 3/4*sqrt(pi);
gamma(zero):= infinity;
gamma(-zero):= -infinity;