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
elseif (a == 0)
  if strcmp(tail,'lower')
    y = 1;
  elseif strcmp(tail,'scaledlower')
    y = exp(x);
  else
    y = 0;
  end
%elseif (a == fix(a)) % a integer
%http://mathworld.wolfram.com/IncompleteGammaFunction.html
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 (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 possibble to scale at the end. Try gammainc(1,1000),
% it take 0 iterations if you scale now.
% moreover, some test fail
    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),0.5420702855281477916858351,-7*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,-eps)
%!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)


