[Bug-gsl] [bug #50459] Non termination of the incomplete gamma function due to floating-point rounding errors

[Fonenantsoa Maurica]

URL:
  <http://savannah.gnu.org/bugs/?50459>

                 Summary: Non termination of the incomplete gamma function due
to floating-point rounding errors
                 Project: GNU Scientific Library
            Submitted by: fmaurica
            Submitted on: Sat 04 Mar 2017 06:56:45 PM UTC
                Category: None
                Severity: 3 - Normal
        Operating System: 
                  Status: None
             Assigned to: None
             Open/Closed: Open
                 Release: 2.3
         Discussion Lock: Any

    _______________________________________________________

Details:

Bug of non termination of the incomplete gamma function due to floating-point
rounding errors in GSL 2.3

The function
int gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result* result)
from gsl/gsl_sf_gamma.h
does not terminate for any value of a and x such that a < -2^52 && 0 < x <=
0.25

For example, the call to gsl_sf_gamma_inc_e in the following program:
#include <stdio.h>
#include <gsl/gsl_sf_gamma.h>
int main (void)
{
  double a = -1E25, x = 0.1;
  gsl_sf_result* result;
  printf ("%i\n", gsl_sf_gamma_inc_e(a, x, result));
  return 0;
}

does not terminate as explained below:
int gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result* result)
{
...
  else
  {
    ...
    double alpha = da; // alpha = 0
    ...
    do
    {
      ...
      alpha -= 1.0; // when alpha reaches -2^-52 (~4.5E15)
                    // then alpha - 1 is rounded to alpha itself
                    // thus alpha will not be updated anymore
    } while(alpha > a); // a = -1E25
  ...
  }
}

To guarantee termination for any possible input, the case a < -2^52 && 0 < x
<= 0.25 should be prevented.





    _______________________________________________________

Reply to this item at:

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

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