GNU Scientific Library  Bugs: bug #41837, bugs in gsl_sf_hyperg_U
You are not allowed to post comments on this tracker with your current authentication level.
bug #41837: bugs in gsl_sf_hyperg_U
Submitter:  Patrick Alken <psa>  
Submitted:  Tue 11 Mar 2014 07:52:19 PM UTC  
Category:  Runtime error  Severity:  3  Normal 
Operating System:  Status:  None  
Assigned to:  None  Open/Closed:  Open 
Release: 
Fri 14 Mar 2014 06:20:06 PM UTC, comment #4: 
Raymond Rogers <rrogers> 
Tue 11 Mar 2014 07:54:53 PM UTC, comment #3: Also see bug #30510 
Patrick Alken <psa> 
Tue 11 Mar 2014 07:54:28 PM UTC, comment #2: Sorry for so many emails... I also tried the GSL SHELL software and it

Patrick Alken <psa> 
Tue 11 Mar 2014 07:54:03 PM UTC, comment #1: Sorry for the typo but for example 1) listed in the previous email:
gsl_sf_hyperg_U(1, 500, 0.1) gives the correct value: ~1.99e3.
The problem happens when calculating gsl_sf_hyperg_U(1, 500, 40): gsl:

Patrick Alken <psa> 
Tue 11 Mar 2014 07:52:19 PM UTC, original submission:
Hi,
I have been using the gsl_sf_hyperg_U function for a while and I found
1) gsl_sf_hyperg_U(1, 500, 0.1) gives incorrect value: 9.3246e+03, which
2) Many domain error bugs. For example, gsl_sf_hyperg_U(1, 6.67, 1) should

Patrick Alken <psa> 
Depends on the following items: None found
Items that depend on this one: None found
Follows 1 latest change.
Date  Changed by  Updated Field  Previous Value  =>  Replaced by 

20140311  psa  Attached File    Added gsl_sf_hyperg_U_bug_report.cpp, #30863 
Caveat: I am not experienced enough to be a maintainer for this code and have just looked into this particular problem.
It needs a more professional eye.
This was analyzed with ddd/gdb on Linux Mint.
I have looked into this and the problem seems to be in:
hyperg_U.c
In this section:
gsl_sf_result sum;
int stat_sum = hyperg_U_finite_sum(N, a, b, x, xeps, &sum);
int stat_inf;
/* Evaluate infinite sum. */
if(fabs(xeps1.0) > 0.5 ) {
stat_inf = hyperg_U_infinite_sum_stable(N, a, bint, b, beps, x, xeps, sum, result);
} else if (1+ab < 0 && 1+ab==floor(1+ab) && beps != 0) {
stat_inf = hyperg_U_infinite_sum_simple(N, a, bint, b, beps, x, xeps, sum, result);
} else {
stat_inf = hyperg_U_infinite_sum_improved(N, a, bint, b, beps, x, xeps, sum, result);
}
What happens is that when ab+1=n then a finite condition 1b+1= n is satisfied and the finite sum succeeds; as it supposed to.
But the program continues to the hyperg_U_infinite_sum_simple where that condition is a no..no and bombs in gamma.c deep inside so it never gets back to check if the finite sum worked.
What should happen is the success of the finite sum should skip the latter tests
Recommendation: If the finite sum succeeded skip the infinite sum tests (and execution).