[Octave-bug-tracker] [bug #47800] gammainc(x, a, "upper") rounds down to zero if output is below eps

[Marco Caliari]

Follow-up Comment #12, bug #47800 (project octave):

Hi Nir. If I understand correctly, you (Gautschi) can avoid overflow/underflow
by computing the series for the scaled version and scaling at the end. The
Legendre continued fraction is instead formula (6.2.6) of Numerical Recipes. I
used (6.2.7) (claimed to be faster), but only for abs(x) > a+1. Anyway, I
found this interesting paper 

https://www.researchgate.net/publication/231929375_A_Set_of_Algorithms_for_the_Incomplete_Gamma_Functions

In section 5.1 it is clear that either the series or the continued fraction
method is used. But it is interesting the first part of the paper, in which
the accurate computation of D(a,x)=x^a*exp(-x)/Gamma(a+1) is discussed. I
already notice that that standard way exp(a*log(x)-x-gammaln(a+1)) is not
always accurate (and may explain the loss of 3-4 digits).
To summarize, I would guess (or hope) that the series/continued fraction
method + accurate computation of D(a,x) would manage all the cases (instead of
currently two series and two continued fractions).  

By the way, "improvement of special functions": what a nice project for next
GSOC!

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?47800>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/