GNU Scientific Library - Bugs: bug #50459, Non termination of the incomplete...
You are not allowed to post comments on this tracker with your current authentication level.
bug #50459: Non termination of the incomplete gamma function due to floating-point rounding errors
Submitter: | Fonenantsoa Maurica <fmaurica> | ||
Submitted: | 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 |
No files currently attached
Depends on the following items: None found
Items that depend on this one: None found
Carbon-Copy List
Follows 1 latest change.
Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|
2017-03-04 | fmaurica | Carbon-Copy | - | Added fmaurica |
Powered by Savane 3.13-f8d8.
Corresponding source code
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.