## 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 = D(a,x)
##  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?
if (x == 0)
  y = 0;
else
  if ((x > 0) && (a > athresh)) 
    lnGa = log (2 * pi *a) / 2+1 / (12 * a) - 1 / (360 * a ^ 3) + ...
           1 / (1260 * a ^ 5) - 1 / (1680 * a ^ 7) + 1 / (1188 * a ^ 9)- ...
           691 / (87360 * a ^ 11) + 1 / (156 * a ^ 13) - ...
           3617 / (122400 * a ^ 15) + 43867 / (244188 * a ^ 17) - ...
           174611 / (125400 * a ^ 19);
    if (a >= x)
      lns = log1p ((a - x) / x);
    else
      lns = -log1p ((x - a) / a);
    endif
    y = exp ((a - x) - a * lns - lnGa);
  else ## standard formula for a not so large
    y = exp(a * log (x) - x - gammaln (a + 1));
  endif
  if ((x < 0) && (a == fix (a)))
% remove spurious imaginary part
    y = real (y);
  endif
endif  
%!test
%! assert(D(10,10),0.12511003572113329898476497894772544708420990097708600,-7*eps)
%%!test
%%! assert(D(100,10),4.864649182067610436458701223522755117545181424119392e-63,-115*eps)
%!test
%! assert(D(100,100),0.0398609968091471352339206494591387451322894644246,-3*eps)
%!test
%! assert(D(200,200),0.0281977276859208217987020623393630120002075706645,-3*eps)
%!test
%! assert(D(200,201),0.0281275549310563187334442096017465724513046028959,-3*eps)
%!test
%! assert(D(1000,1000),0.01261461134872149971803693647457875764715489427,-4*eps)
%!test
%! assert(D(10,-10),6.069903492836947893782490257188118487200593203967e7,-8*eps)
%!test
%! assert(D(10,-10.1),7.410115412614356278962877535200296949490313010e7,-20*eps)
%!test
%! assert(D(9,-10),-6.069903492836947893782490257188118487200593203967e7,-8*eps)
%!test
%! assert(D(10.1,-10),5.741783267595397392122436851864044460198744897912388e7+1i*1.865618475142280818086513704349404376045316209289692e7,-25*eps)
