## Copyright (C) 2016 Marco Caliari, Nir Krakauer, Andrew Lachlan
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

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";
  endif

  if (! isscalar (x) || ! isscalar (a))
    [err, x, a] = common_size (x, a);
    if (err > 0)
      error ("ngammainc: x, a must be of common size or scalars");
    endif
  endif

  #initialize output array
  y = nan (size(x));

  ##different x, a combinations are handled by different subfunctions
  i_done = false (size(x)); #keeps track of which combinations have already been calculated

  #case 1: x == 0
  ii = (x == 0);
  if any (ii(:))
    y(ii) = ngammainc_x0 (tail);
    i_done(ii) = true;
  endif

  #case 2: a == 0
  ii = ~i_done & (a == 0);
  if any (ii(:))
    y(ii) = ngammainc_a0 (x(ii), tail);
    i_done(ii) = true;
  endif

  #case 3: a == 1
  ii = ~i_done & (a == 1);
  if any (ii(:))
    y(ii) = ngammainc_a1 (x(ii), tail);
    i_done(ii) = true;
  endif

  #case 4: positive integer a; exp (x) and a! both under 1/eps
  ii = ~i_done & (a > 1) & (a == fix (a)) & (x <= 36) & (a <= 18);
  if any (ii(:))
    y(ii) = ngammainc_an (x(ii), a(ii), tail);
    i_done(ii) = true;
  endif

  #case 5: x relatively small
  ii = ~i_done & (x + 0.25 < a | x < 0);
  if any (ii(:))
    y(ii) = ngammainc_s (x(ii), a(ii), tail);
    i_done(ii) = true;
  endif

  #case 6: x positive and large relative to a
  ii = ~i_done;
  if any (ii(:))
    y(ii) = ngammainc_l (x(ii), a(ii), tail);
    i_done(ii) = true;
  endif

endfunction




  ##subfunctions to handle each case:

  #x == 0
  function y = ngammainc_x0 (tail)
    y = 1;
    if (strcmpi (tail,"upper") || strcmpi (tail,"scaledupper"))
      y = 0;
    endif
  endfunction

  #a == 0
  function y = ngammainc_a0 (x, tail)
    if strcmpi (tail,"lower")
      y = 1;
    elseif strcmpi (tail,"scaledlower")
      y = exp (x);
    else
      y = 0;
    endif
  endfunction

  #a == 1
  function y = ngammainc_a1 (x, tail)
    if strcmpi (tail,"lower")
      y = 1 - exp (-x);
    elseif strcmpi (tail,"scaledlower")
      if (abs (x) < 1/2)
        y = expm1 (x) / x;
      else
        y = (exp (x) - 1) / x;
      endif
    elseif strcmpi (tail,"upper")
      y = exp (-x);
    else strcmpi (tail,"scaledupper")
      y = 1 / x;
    endif
  endfunction

  #positive integer a; exp (x) and a! both under 1/eps
  ## uses closed-form expressions for nonnegative integer a
  ## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html
  function y = ngammainc_an (x, a, tail)
    y = t = ones (size(x));
    i = 1;
    while any(a(:) > i)
      jj = (a > i);
      t(jj) .*= (x(jj) / i);
      y(jj) += t(jj);
      i++;
    endwhile
    if strcmpi (tail,"upper")
      y .*= exp (-x);
    elseif strcmpi (tail,"lower")
      y = 1 - exp (-x) .* y;
    elseif strcmpi (tail,"scaledupper")
      y .*= D (a,x) .* exp (-x);
    elseif strcmpi (tail,"lower")
      y = D (a,x) .* (1 - exp (-x) .* y);
    endif
  endfunction

  #x + 0.25 < a | x < 0
    ## Numerical Recipes in Fortran 77 (6.2.5)
    ## series     
  function y = ngammainc_s (x, a, tail)
    if (strcmpi (tail,"scaledlower") || strcmpi (tail,"scaledupper"))
      y = ones (size(x));
      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);
    endif
    n = 1;
    while any (abs(term(:)) > abs(y(:)) * eps)
      ## y can be zero from the beginning (gammainc (1,1000))
      jj = abs(term) > abs(y) * eps;
      n += 1;
      y(jj) += term(jj);
      term(jj) .*= x(jj) ./ (a(jj) + n);
    endwhile
    if strcmpi (tail,"upper")
      y = 1 - y;
    elseif strcmpi (tail,"scaledupper")
      y = 1 ./ D (a,x) - y;
    endif
  endfunction

  #x positive and large relative to a
    ## NRF77 (6.2.7)
    ## Gamma (a,x)/Gamma (a)
    ## Lentz's algorithm
  function y = ngammainc_l (x, a, tail)
    tiny = 2^(-100);
    y = tiny;
    C0 = y;
    D0 = 0;
    bj = x + 1 - 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 any(abs ((Deltaj(:) - 1) ./ y(:)) > eps)
      Dj = bj + aj .* D0;
      Cj = bj + aj ./ C0;
      Dj = 1 ./ Dj;
      Deltaj = Cj .* Dj;
      y .*= Deltaj;
      D0 = Dj;
      C0 = Cj;
      bj += 2;
      aj = -j * (j - a);
      j += 1;
    endwhile
    if strcmpi (tail,"upper")
      y .*= D (a,x);
    elseif strcmpi (tail,"lower")
      y = 1 - y .* D (a,x);
    elseif strcmpi (tail,"scaledlower")
      y = 1 ./ D (a,x) - y;
    endif
endfunction

function y = D (a,x)
  y = zeros (size(x));
  j = x < 0;
  if any (j(:))
    ## to prevent imaginary parts when x negative and a integer
    y(j)  = x(j).^a(j) .* exp (-x(j) - gammaln (a(j) + 1));
  endif
  j = x > 0;
  if any (j(:))
    ## 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(j) = exp (a(j) .* log (x(j)) - x(j) - gammaln (a(j) + 1));
  endif
endfunction




%!test
%! assert (ngammainc([0 0],0),[1 1]);
%!test
%! assert (ngammainc (0,[0; 0],"upper"),[0; 0]);
%!test
%!error (ngammainc([0 0],[0; 0]));
%!test
%! assert (ngammainc ([1 10],[1000 0]),[0 1]);
%!test
%! assert (ngammainc ([3 2 36],[2 3 18], "upper"),[4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]);
%!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 (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 (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);