function y = ngammainc(x,a,tail)
% incomplete (upper) 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
% Regions roughly as in Gautschi, W. (1979). 
% A computational procedure for incomplete gamma functions. ACM
% Transactions on Mathematical Software (TOMS), 5(4), 466-481.
% P(a,x) = gamma(a,x)/Gamma(a), upper
% 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower
if (nargin == 2)
  tail = 'lower';
end
% some corner cases
if (x == 0)
  y = 1;
  if (strcmp(tail,'upper') || strcmp(tail,'scaledupper'))
    y = 0;
  end
%closed-form expressions for nonnegative integer a -- http://mathworld.wolfram.com/IncompleteGammaFunction.html
elseif (a == 0)
  if strcmp(tail,'lower')
    y = 1;
  elseif strcmp(tail,'scaledlower')
    y = exp(x);
  else
    y = 0;
  end
elseif (a == 1)
  if strcmp(tail,'lower')
    y = 1-exp(-x);
  elseif strcmp(tail,'scaledlower')
    if (abs(x) < 1/2)
      y = expm1(x)/x;
    else
      y = (exp(x)-1)/x;
    end
  elseif strcmp(tail,'upper')
    y = exp(-x);
  else strcmp(tail,'scaledupper')
    y = 1/x;
  end
elseif (a > 1) && (a == fix(a)) && (x <= 36) && (a <= 18) %positive integer a; exp(x) and a! both under 1/eps
  y = 1; t = 1;
  for i = 1:(a-1)
    t = t * (x/i);
    y = y + t;
  end
  if strcmp(tail,'upper')
    y = exp(-x) * y;
  elseif strcmp(tail,'lower')
    y = 1 - exp(-x) * y;
  elseif strcmp(tail,'scaledupper')
    y = D(a,x) * exp(-x) * y;
  elseif strcmp(tail,'lower')
    y = D(a,x) * (1 - exp(-x) * y);
  end
elseif (x + 0.25 < a || x < 0)
% NRF77 (6.2.5)
% series     
  if (strcmp(tail,'scaledlower') || strcmp(tail,'scaledupper'))
    y = 1;
    term = x/(a+1);
  else
% Of course it is possible to scale at the end. But some tests fail.
% And try gammainc(1,1000), it take 0 iterations if you scale now.
    y = D(a,x);
    term = y*x/(a+1);
  end
  n = 1;
  while (abs(term) > abs(y)*eps)
% y can be zero from the beginning (gammainc(1,1000))
    n = n+1;
    y = y+term;
    term = term*(x/(a+n));
  end
  if strcmp(tail,'upper')
    y = 1-y;
  elseif strcmp(tail,'scaledupper')
    y = 1/D(a,x)-y;
  end
else
% NRF77 (6.2.7)
% Gamma(a,x)/Gamma(a)
% Lentz's algorithm
% here x is positive
  tiny = 2^(-100);
  y = tiny;
  C0 = y;
  D0 = 0;
  bj = x+1-a;
  aj = a;
% if 'lower' or 'upper', it should be 
%
% aj = a*exp(a*log(x)-x-gammaln(a+1))
%
% and no scaling at the end. But it suffers cancellation
% Better start with aj = a always.
  aj = a;
  Deltaj = 0;
  j = 1;
  while (abs((Deltaj-1)/y) > eps)
    Dj = bj+aj*D0;
    Cj = bj+aj/C0;
    Dj = 1/Dj;
    Deltaj = Cj*Dj;
    y = y*Deltaj;
    D0 = Dj;
    C0 = Cj;
    bj = bj+2;
    aj = -j*(j-a);
    j = j+1;
  end
  if strcmp(tail,'upper')
    y = y*D(a,x);
  elseif strcmp(tail,'lower')
    y = 1-y*D(a,x);
  elseif strcmp(tail,'scaledlower')
    y = 1/D(a,x)-y;
  end
end
end

function y = D(a,x)
  if (x == 0)
    y = 0;
  elseif (x < 0)
% to prevent imaginary parts when x negative and a integer
     y  = x^a*exp(-x-gammaln(a+1));
  else
% Here, Temme's formula should be used
% N. M. Temme, A set of algorithms for the incomplete Gamma functions,
% Probab. Engrg. Inform. Sci., 8, 1994.
    y = exp(a*log(x)-x-gammaln(a+1));
  end
end
%!test
%! assert(ngammainc(0,0),1)
%!test
%! assert(ngammainc(0,0,'upper'),0)
%!test
%! assert(ngammainc(10,10),1-(5719087/567)*exp(-10),-eps)
%!test
%! assert(ngammainc(10,10, 'upper'),(5719087/567)*exp(-10),-eps)
%!test
%! assert(ngammainc(-10,10),3.112658265341493126871617e7,-11*eps)
%!test
%! % Here Matlab fails
%! assert(isreal(ngammainc(-10,10)),true)
%!test
%! assert(ngammainc(60,6,'upper'),6.18022358081160257327264261e-20,-10*eps)
%!test
%! assert(ngammainc(-10,10.1,'upper'),...
%!        -2.9582761911890713293e7-1i*9.612022339061679758e6,-30*eps)
%!test
%! % Here Matlab is better
%! assert(ngammainc(751,750,'upper'),0.4805914320558831327179457887,-512*eps)
%!test
%! % Here Matlab fails
%! assert(ngammainc(-100,1,'upper'),exp(100),-eps)
%!test
%! assert(ngammainc(10,0),1)
%!test
%! assert(ngammainc(200,200,'upper'),0.49059658199276367497217454,-39*eps)
%!test
%! assert(ngammainc(200,200),0.5094034180072363250278254595745270431549,-38*eps)
%!test
%! assert(ngammainc(200,200,'scaledupper'),17.3984438553791505135122900,-eps)
%!test
%! assert(ngammainc(200,200,'scaledlower'),18.0654066767792216430658593,-78*eps)
%!test
%! assert(ngammainc(200,201,'upper'),0.5187943096786844967708766027648,-808*eps)
%!test
%! assert(ngammainc(201,200,'upper'),0.4624924490827670952491373666731,-698*eps)
%!test
%! assert(ngammainc(1,1000),0)
%!test
%! assert(ngammainc(200,201,'scaledupper'),18.4904360746560462660798514700,-1678*eps)
%!test
%! % here we are very good (no D(a,x)) involved
%! assert(ngammainc(1000,1000.5,'scaledlower'),39.4846753958367227115137,-2*eps)
%!test
%! assert(ngammainc(709,1000,'upper'), 0.999999999999999999999999543586098796378631194583297890832, -eps)
%!test
%! assert(ngammainc(3,2, 'upper'),4/exp(3))
%!test
%! assert(ngammainc(2,3, 'upper'),5*exp(-2))
%!test
%! assert(ngammainc(36,18, 'upper'),(4369755579265807723/2977975)/exp(36))
