[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [bug #50459] Non termination of the incomplete gamma function
From: |
Fonenantsoa Maurica |
Subject: |
[Bug-gsl] [bug #50459] Non termination of the incomplete gamma function due to floating-point rounding errors |
Date: |
Sat, 4 Mar 2017 13:56:46 -0500 (EST) |
User-agent: |
Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:51.0) Gecko/20100101 Firefox/51.0 |
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/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [bug #50459] Non termination of the incomplete gamma function due to floating-point rounding errors,
Fonenantsoa Maurica <=