This file is the GSL bug tracking system. The CVS version of this file should be kept up-to-date. ---------------------------------------------------------------------- BUG#1 -- gsl_sf_hyperg_2F1_e fails for some arguments From: keith.briggs@bt.com Subject: gsl_sf_hyperg_2F1 bug report Date: Thu, 31 Jan 2002 12:30:04 -0000 gsl_sf_hyperg_2F1_e fails with arguments (1,13,14,0.999227196008978,&r). It should return 53.4645... . #include #include int main (void) { gsl_sf_result r; gsl_sf_hyperg_2F1_e (1,13,14,0.999227196008978,&r); printf("r = %g %g\n", r.val, r.err); } NOTES: The program overflows the maximum number of iterations in gsl_sf_hyperg_2F1, due to the presence of a nearby singularity at (c=a+b,x=1) so the sum is slowly convergent. The exact result is 53.46451441879150950530608621 as calculated by gp-pari using sumpos(k=0,gamma(a+k)*gamma(b+k)*gamma(c)*gamma(1)/ (gamma(c+k)*gamma(1+k)*gamma(a)*gamma(b))*x^k) The code needs to be extended to handle the case c=a+b. This is the main problem. The case c=a+b is special and needs to be computed differently. There is a special formula given for it in Abramowitz & Stegun 15.3.10 As reported by Lee Warren another set of arguments which fail are: #include #include int main (void) { gsl_sf_result r; gsl_sf_hyperg_2F1_e (-1, -1, -0.5, 1.5, &r); printf("r = %g %g\n", r.val, r.err); } The correct value is -2. See also, From: Olaf Wucknitz To: bug-gsl@gnu.org Subject: [Bug-gsl] gsl_sf_hyperg_2F1 Hi, I am having a problem with gsl_sf_hyperg_2F1. With the parameters (-0.5, 0.5, 1, x) the returned values show a jump at x=0.5. For x<0.5 the results seem to be correct, while for x>0.5 they aren't. The function gsl_sf_hyperg_2F1_e calls hyperg_2F1_series for x<0.5, but hyperg_2F1_reflect for x>0.5. The latter function checks for c-a-b being an integer (which it is in my case). If I change one of the parameters a,b,c by a small amount, the other branch of the function is taken and the results are correct again. Unfortunately I know too little about the numerics of 2F1 to suggest a patch at the moment. Regards, Olaf Wucknitz -- Joint Institute for VLBI in Europe wucknitz@jive.nl ---------------------------------------------------------------------- BUG#44 -- gamma_inc_P and gamma_inc_Q only satisfy P+Q=1 within errors The sum of gamma_inc_P and gamma_inc_Q doesn't always satisfy the identity P+Q=1 exactly (although it is correct within errors), due the slightly different branch conditions for the series and continued fraction expansions. These could be made identical so that P+Q=1 exactly. #include #include int main (void) { gsl_sf_result r1, r2; double a = 0.3, x = 1.0; gsl_sf_gamma_inc_P_e (a, x, &r1); gsl_sf_gamma_inc_Q_e (a, x, &r2); printf("%.18e\n", r1.val); printf("%.18e\n", r2.val); printf("%.18e\n", r1.val + r2.val); } $ ./a.out 9.156741562411074842e-01 8.432584375889111417e-02 9.999999999999985567e-01 ---------------------------------------------------------------------- BUG#49 - fdjacobian step size causes problems From: eknecronzontas To: help-gsl@gnu.org Subject: [Help-gsl] Bug (or not?) in multiroots/fdjacobian.c Date: Thu, 8 Jun 2006 13:38:49 -0700 (PDT) Since I can't quite decide if this is a bug, I'll post it here first... The stepsize for finite-differencing in gsl-1.8/multiroots/fdjacobian.c is specified by the lines > double xj = gsl_vector_get (x, j); > double dx = epsrel * fabs (xj); This is, of course mathematically correct. Nevertheless, the behavior is less than ideal if one of the elements of the vector 'x' is sufficiently small. This occurs often if one component of the root happens to be zero. If this occurs, then 'dx' can be so small that one of the columns of the Jacobian is identically zero. This immediately leads to NAN's which are not easy to trace. --------------------------------------------------------------------- BUG#50 - gsl_linalg_solve_symm_tridiag requires positive definite matrix A zero on the diagonal will cause NaNs even though a reasonable solution could be computed in principle. #include int main (void) { double d[] = { 0.00, 1.21, 0.80, 1.55, 0.76 } ; double e[] = { 0.82, 0.39, 0.09, 0.68 } ; double b[] = { 0.07, 0.62, 0.81, 0.11, 0.65} ; double x[] = { 0.00, 0.00, 0.00, 0.00, 0.00} ; gsl_vector_view dv = gsl_vector_view_array(d, 5); gsl_vector_view ev = gsl_vector_view_array(e, 4); gsl_vector_view bv = gsl_vector_view_array(b, 5); gsl_vector_view xv = gsl_vector_view_array(x, 5); gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector); gsl_vector_fprintf(stdout, &xv.vector, "% .5f"); d[0] += 1e-5; gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector); gsl_vector_fprintf(stdout, &xv.vector, "% .5f"); } $ ./a.out nan nan nan nan nan 0.13626 0.08536 1.03840 -0.60009 1.39219 ---------------------------------------------------------------------- BUG#52 - beta functions do not handle negative arguments The beta functions (complete and incomplete) are well defined for non-integer negative arguments in terms of the gamma function but return domain errors or nans. The relation of I_x(a,b,x) = (1/a) x^a 2F1(a,1-b,a+1,x)/B(a,b) could be documented even if it is not implemented. DONE for beta function, still needed for incomplete beta function. ---------------------------------------------------------------------- BUG#57 - incorrect rounding error in deriv functions? Needs a factor of 1/h^2 ? From: Rene Girard To: help-gsl@gnu.org Subject: [Help-gsl] Origin of 2nd round-off term "dy" in function central_deriv.c double dy = GSL_MAX(fabs(r3),fabs(r5))* fabs(x)*GSL_DBL_EPSILON I understand the rest of the function very well and I am able to derive the equation for the optimal stepsize "h_opt" in function "gsl_deriv_central.c"; however, I am trying to look for a derivation for the expression of "dy". ---------------------------------------------------------------------- Last assigned bug number = 57