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
%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)
% 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a)

%parse the input flag
if (nargin == 2)
  tail = 'lower';
end
if (numel(tail) >= 6) && strcmp(tail(1:6),'scaled')
  scaled = true; %return P or Q times scale = gamma(a+1)*exp(x)*x^-a
  tail = tail(7:end);
else
  scaled = false;
end


if (x > 1.5) && (x > 20) & (x < 300) && ((a - x) <= 0.25) %Gautschi, Sec. 4.3 [Legendre continued fraction for y = (x - a + 1) * x^-a * exp(x) * Gamma(a,x) -- works for large positive inputs, but unfortunately incurs a good bit of roundoff error]
  p = 0; q = (x - a - 1) * (x - a + 1); r = 4 * (x - a + 1); s = 1 - a; rho = 0; t = 1; y = t;
  while abs(t / y) > eps
    p += s;
    q += r;
    r += 8;
    s += 2;
    tau = p * (1 + rho);
    rho = tau / (q - tau);
    t *= rho;
    y += t;
  end
  if scaled
    y *= a / (x - a + 1); %scale*Q
    if strcmp(tail, 'lower')
      y = exp(x - a*log(x) + gammaln(a+1)) - y; %scale*P
    end
  else
    y *= exp(-x + a*log(x) - gammaln(a)) / (x - a + 1); %Q
    if strcmp(tail, 'lower')
      y = 1 - y; %P
    end
  end
  return
end

if (x >= 0) && ((a - x) > 0.25) %Gautschi, Sec. 5 [power series for y = gamma(a+1) * x^-a * exp(x) * P = scale*P]
  t = 1;
  y = t;
  n = 0;
  while abs(t / y) > eps
    n++;
    t *= x / (a + n);
    y += t;
  end
  if scaled
    if strcmp(tail, 'upper')
      y = exp(x - a*log(x) + gammaln(a+1)) - y; %scale*Q
    end
  else
    y *= exp(-x + a*log(x) - gammaln(a+1)); %P
    if strcmp(tail, 'upper')
      y = 1 - y; %Q
    end
  end
  return
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 ~scaled
      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')
    if ~scaled
      y = 1-y;
    else
      y = exp(-a*log(x)+x+gammaln(a+1))-y;
    end
  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,-3*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,-300*eps)
%!test
%! assert(newgammainc(200,201,'upper'),0.518794309678684496770876602764835968845235402065399405783,-1000*eps)
%!test
%! assert(newgammainc(201,200,'upper'),0.462492449082767095249137366673156681609762469721457129052,-1000*eps)