function y = newgammainc(x,a,tail)
% incomplete gamma function
%
% P(a,x) = 1/Gamma(a) * int_x^\infty exp(-t)*t^(a-1)dt
%
% a>=0, x real
%
% tail is either upper or lower 
% Numerical Recipes in Fortran77, \S 6.2
% P(a,x) = gamma(a,x)/Gamma(a)
% 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a)
if (nargin == 2)
  tail = 'lower';
end
if (abs(x) <= a+1 || x <= 0)
% some corner cases
  if (x == 0)
    y = 0;
  elseif (a == 0)
    y = 1;	 
  elseif (a == 1)
    y = 1-exp(-x);
  elseif (abs(x) >= 300) % to be fixed
% (6.2.5)
% series gamma(a,x)/Gamma(a) for large numbers
% take exp(log(each term of the series))
% not very accurate, Kahan summation does not help
    term = exp(a*log(x)-x-gammaln(a+1));
    y = 0;
    n = 0;
    while (abs(term) > eps*abs(y))
      y = y+term;
      n = n+1;
      term = exp((a+n)*log(x)-x-gammaln(a+1+n));
    end
  else
% (6.2.5)
% series     
    if ~strcmp(tail(1:5),'scale')
      gammaa = gamma(a);
      gamma1n = a*gammaa;
      xn = exp(-x)*x^a;
      % formula xn = exp(-x+a*log(x)) does not work with x<0, e.g. x=-10,a=10
      % but for positive x it may help
      y = xn/gamma1n;
      gamma1n = (a+1)*gamma1n;
      xn = xn*x;
      n = 1;
      while (abs(xn/gamma1n) > eps*abs(y))
	y = y+xn/gamma1n;
	gamma1n = (a+n+1)*gamma1n;
	xn = xn*x;
	n = n+1;
      end
    else
% a scaled version was required
% series, but compute a*gamma(a,x)/(exp(-x)*x^a) and simplify 
% Gamma(a)/Gamma(a+1+n) by 1/(a+n)!
      xn = 1;
      y = xn;
      rgammainv = (a+1);
      xn = xn*x;
      n = 1;
      while (abs(xn/rgammainv) > eps*abs(y))
	y = y+xn/rgammainv;
	rgammainv = (a+n+1)*rgammainv;
	xn = xn*x;
	n = n+1;
      end
    end  
  end
  if strcmp(tail,'upper')
     y = 1-y;
  elseif strcmp(tail,'scaledupper')
     y = exp(-a*log(x)+x+gammaln(a+1))-y;
  end
else
% (6.2.7)
% Gamma(a,x)/Gamma(a)
% Lentz's algorithm
% here x is positive
% scaled versions not implemented yet
% Could it be that continued fraction (6.2.6), although slower, is
% more accurate?
    tiny = 2^(-100);
    b0 = 0;
    f0 = tiny;
    C0 = f0;
    D0 = 0;
    bj = x+1-a;
    aj = 1;
    Deltaj = 0;
    j = 1;
    while (abs(Deltaj-1) > eps)
      Dj = bj+aj*D0;
      Cj = bj+aj/C0;
      Dj = 1/Dj;
      Deltaj = Cj*Dj;
      fj = f0*Deltaj;
      D0 = Dj;
      C0 = Cj;
      f0 = fj;
      bj = x+2*j+1-a;
      aj = -j*(j-a);
      j = j+1;
    end
% try to change next line to
% y = exp(-x-gammaln(a)+a*log(x))*fj;
% one more test will fail!
    y = exp(a*log(x)-x-gammaln(a))*fj;
    if strcmp(tail,'lower')
      y = 1-y;
    end
end
%!test
%! assert(newgammainc(0,0),0)
%!test
%! assert(newgammainc(0,0,'upper'),1)
%!test
%! assert(newgammainc(10,10),0.5420702855281477916858351,-2*eps)
%!test
%! % Here Matlab fails
%! assert(newgammainc(-10,10),3.112658265341493126871617e7,-2*eps)
%!test
%! assert(newgammainc(60,6,'upper'),6.18022358081160257327264261e-20,-13*eps)
%!test
%! assert(newgammainc(-10,10.1,'upper'),...
%!        -2.9582761911890713293e7-1i*9.612022339061679758e6,-27*eps)
%!test
%! % Here Matlab is better
%! assert(newgammainc(751,750,'upper'),0.4805914320558831327179457887,-1800*eps)
%!test
%! % Here Matlab fails
%! assert(newgammainc(-100,1,'upper'),exp(100),-eps)
%!test
%! assert(newgammainc(10,0),1,-eps)
%!test
%! assert(newgammainc(200,200,'upper'),0.49059658199276367497217454,-eps)
