## Copyright (C) 2016 Marco Caliari
## Copyright (C) 2016 Nir Krakauer
##
## This file is part of Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {} {} gammainc (@var{x}, @var{a})
## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"lower\")
## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"upper\")
## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"scaledlower\")
## @deftypefnx {} {} gammainc (@var{x}, @var{a}, \"scaledupper\")
## Compute the normalized incomplete gamma function.
## 
## This is defined as
## @tex
## $$
##  \\gamma (x, a) = {1 \\over {\\Gamma (a)}}\\displaystyle{\\int_0^x t^{a-1} e^{-t} dt}
## $$
## @end tex
## @ifnottex
## 
## @example
## @group
##                                 x
##                        1       /
## gammainc (x, a) = ---------    | exp (-t) t^(a-1) dt
##                   gamma (a)    /
##                             t=0
## @end group
## @end example
## 
## @end ifnottex
## with the limiting value of 1 as @var{x} approaches infinity.
## The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and
## @nospell{Stegun} (6.5.1).
## 
## If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned
## for each element of @var{x} and vice versa.
## 
## If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and
## @var{a} must agree, and @code{gammainc} is applied element-by-element.
## 
## By default the incomplete gamma function integrated from 0 to @var{x} is
## computed.  If @qcode{\"upper\"} is given then the complementary function
## integrated from @var{x} to infinity is calculated.
##
## Algorithms used include the following.
## 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.
##
## @seealso{gamma, gammaln}
## @end deftypefn

## P(a,x) = gamma(a,x)/Gamma(a), upper
## 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower

function y = gammainc (x,a,tail)
  if (nargin == 2)
    tail = "lower";
  endif

  if (! isscalar (x) || ! isscalar (a))
    [err, x, a] = common_size (x, a);
    if (err > 0)
      error ("gammainc: 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)); # Track which elements have been calculated.

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

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

  ## Case 3: a == 1
  ii = ~i_done & (a == 1);
  if any (ii(:))
    y(ii) = gammainc_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) = gammainc_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) = gammainc_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) = gammainc_l (x(ii), a(ii), tail);
    i_done(ii) = true;
  endif

endfunction




  ##subfunctions to handle each case:

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

  ## a == 0
  function y = gammainc_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 = gammainc_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
      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 = gammainc_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 (x,a) .* exp (-x);
    elseif strcmpi (tail,"lower")
      y = D (x,a) .* (1 - exp (-x) .* y);
    endif
  endfunction

  ## x + 0.25 < a | x < 0
    ## Numerical Recipes in Fortran 77 (6.2.5)
    ## series     
  function y = gammainc_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 (x,a);
      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 (x,a) - y;
    endif
  endfunction

  ## x positive and large relative to a
    ## NRF77 (6.2.7)
    ## Gamma (a,x)/Gamma (a)
    ## Lentz's algorithm
  function y = gammainc_l (x, a, tail)
    n = numel (x);
    y = nan (size(x));
    for i = 1:n
      y(i) = __gammainc_lentz(x(i), a(i));
    endfor
    #tiny = 2^(-100);
    #y = tiny*ones(size(x));
    #Cj = y;
    #Dj = zeros(size(x));
    #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 = zeros(size(x));
    #j = 1;
    #while any (abs ((Deltaj(:) - 1) ./ y(:)) > eps)
    #  jj = abs ((Deltaj - 1) ./ y) > eps;
    #  Cj(jj) = bj(jj) + aj(jj) ./ Cj(jj);
    #  Dj(jj) = 1 ./ (bj(jj) + aj(jj) .* Dj(jj));
    #  Deltaj(jj) = Cj(jj) .* Dj(jj);
    #  y(jj) .*= Deltaj(jj);
    #  bj(jj) += 2;
    #  aj(jj) = -j * (j - a(jj));
    #  j++;
    #endwhile
    if strcmpi (tail,"upper")
      y .*= D (x,a);
    elseif strcmpi (tail,"lower")
      y = 1 - y .* D (x,a);
    elseif strcmpi (tail,"scaledlower")
      y = 1 ./ D (x,a) - y;
    endif
endfunction


function y = D(x,a)
##  Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large
##
## L. Knusel, Computation of the Chi-square and Poisson distribution,
## SIAM J. Sci. Stat. Comput., 7(3), 1986
## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41
  athresh = 10; ## FIXME: can be better tuned?
  y = zeros(size(x));
  i_done = false (size (x));
  i_done(x == 0) = true;
  ii = ~i_done & (x > 0) & (a > athresh) & (a >= x);
  if any (ii(:))
    lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ...
           1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ...
           1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ...
           691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ...
           3617 ./ (122400 * a(ii) .^ 15) + ...
           43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19);
    lns = log1p ((a(ii) - x(ii)) ./ x(ii));
    y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa);
    i_done(ii) = true;
  endif
  ii = ~i_done & (x > 0) & (a > athresh) & (a < x);
  if any (ii(:))
    lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ...
           1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ...
           1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ...
           691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ...
           3617 ./ (122400 * a(ii) .^ 15) + ...
           43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19);
    lns = -log1p ((x(ii) - a(ii)) ./ a(ii));
    y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa);
    i_done(ii) = true;
  endif
  ii = ~i_done & ((x <= 0) | (a <= athresh));
  if any (ii(:)) ## standard formula for a not so large
    y(ii) = exp (a(ii) .* log (x(ii)) - x(ii) - gammaln (a(ii) + 1));
    i_done(ii) = true;
  endif
  ii = (x < 0) & (a == fix (a));
  if any(ii(:)) ## remove spurious imaginary part
    y(ii) = real (y(ii));
  endif
endfunction

%!test
%! assert (gammainc([0 0],0),[1 1]);
%!test
%! assert (gammainc (0,[0; 0],"upper"),[0; 0]);
%!test
%!error (gammainc([0 0],[0; 0]));
%!test
%! assert (gammainc ([1 10],[1000 0]),[0 1]);
%%!test
%%! assert (gammainc ([3 2 36],[2 3 18], "upper"),[4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]);
%!test
%! assert (gammainc (10,10),1 - (5719087 / 567) * exp (-10),-eps);
%!test
%! assert (gammainc (10,10, "upper"),(5719087 / 567) * exp (-10),-eps);
%!test
%! assert (gammainc (-10,10),3.112658265341493126871617e7,-2 * eps);
%!test
%! % Here Matlab fails
%! assert (isreal (gammainc (-10,10)),true);
%!test
%! assert (gammainc (60,6,"upper"),6.18022358081160257327264261e-20,-10 * eps);
%!test
%! assert (gammainc (-10,10.1,"upper"),...
%!        -2.9582761911890713293e7-1i * 9.612022339061679758e6,-30 * eps);
%!test
%! % Here Matlab is better
%! assert (gammainc (751,750,"upper"),0.4805914320558831327179457887,-12 * eps);
%!test
%! % Here Matlab fails
%! assert (gammainc (-100,1,"upper"),exp (100),-eps);
%!test
%! assert (gammainc (200,200,"upper"),0.49059658199276367497217454,-4 * eps);
%!test
%! assert (gammainc (200,200),0.509403418007236325027825459574527043,-3 * eps);
%!test
%! assert (gammainc (200,200,"scaledupper"),17.3984438553791505135122900,-eps);
%!test
%! assert (gammainc (200,200,"scaledlower"),18.065406676779221643065,-6 * eps);
%!test
%! assert (gammainc (200,201,"upper"),0.5187943096786844967708766027,-2 * eps);
%!test
%! assert (gammainc(201,200,"upper"),0.46249244908276709524913736667,-7 * eps);
%!test
%! assert (gammainc (200,201,"scaledupper"),18.4904360746560462660798514,-eps);
%!test
%! % here we are very good (no D (x,a)) involved
%! assert (gammainc(1000,1000.5,"scaledlower"),39.48467539583672271,-2* eps);
%!test
%! assert (gammainc (709,1000,"upper"), 0.99999999999999999999999954358, -eps);
