On Tuesday, August 24, 2010, Andrew Benson wrote:
> Under GSL 1.14 (compiled with either gcc version 4.4.4 or gcc version
> 4.1.2) the gsl_sf_hyperg_2F1 function crashes with the following
> arguments:
>
> gsl_sf_hyperg_2F1(-10.34,2.05,3.05,0.1725)
>
> The expected answer (computed using Maple for example) is 0.3104735522
>
> I've attached an example case:
>
> $ g++ test.cpp -o test.x -lgsl -lgslcblas
> $ test.x
> gsl: hyperg_2F1.c:755: ERROR: error
> Default GSL error handler invoked.
> Abort
>
> The crash seems to occur only when the first argument is less than or equal
> to -10. For example,
>
> gsl_sf_hyperg_2F1(-10.0,2.05,3.05,0.1725)
>
> crashes, but
>
> gsl_sf_hyperg_2F1(-9.99999999999,2.05,3.05,0.1725)
>
> evaluates correctly to 0.321419346301974773.
Some further investigation shows that this combination of parameters is
falling through all of the various methods for computing the result in
gsl_sf_hyperg_2F1_e and hitting the "/* We give up. */" block at the end.
This seems to be related to http://savannah.gnu.org/bugs/?24812 - if I make
the changes suggested there (patch attached) then my test code works
successfully:
$ test.x
= 0.310473552213130111
I haven't checked that this doesn't break anything else (beyond running "make
check" on the newly compiled code).
--- specfunc/hyperg_2F1.c 2010-03-10 02:57:14.000000000 -0800
+++ specfunc/hyperg_2F1_new.c 2010-08-24 13:56:20.000000000 -0700
@@ -734,14 +734,14 @@ gsl_sf_hyperg_2F1_e(double a, double b,
return hyperg_2F1_luke(a, b, c, x, result);
}
- if(GSL_MAX_DBL(fabs(a),1.0)*fabs(bp)*fabs(x) < 2.0*fabs(c)) {
+ if(GSL_MAX_DBL(fabs(ap),1.0)*fabs(bp)*fabs(x) < 2.0*fabs(c)) {
/* If c is large enough or x is small enough,
* we can attempt the series anyway.
*/
return hyperg_2F1_series(a, b, c, x, result);
}
- if(fabs(bp*bp*x*x) < 0.001*fabs(bp) && fabs(a) < 10.0) {
+ if(fabs(bp*bp*x*x) < 0.001*fabs(bp) && fabs(ap) < 10.0) {
/* The famous but nearly worthless "large b" asymptotic.
*/
int stat = gsl_sf_hyperg_1F1_e(a, c, bp*x, result);
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
