
From:  Peter Breitenlohner 
Subject:  [Buggsl] Bug in gsl1.8 trigonometric restriction functions + patches 
Date:  Mon, 11 Sep 2006 16:42:11 +0200 (CEST) 
Hi, while looking at the GSL (gsl1.8) source code I noticed two bugs plus a few problems in the `Trigonometric Restriction Functions': As demonstrated by the attached small test program `gslbug.c': ==================== I.A. The functions `gsl_sf_angle_restrict_symm' & Co rightly detect a loss of accuracy (GSL_ELOSS) for excessive positive arguments, but fail to do so for equally excessive negative arguments. I.B. The result of the function `gsl_sf_angle_restrict_pos' may be negative due to subtle details of the floating point arithmetic.  The first attached `gsl1.8mod1patch' fixes problem A and the most obvious part of problem B. ==================== II.1. The error estimates in the functions `gsl_sf_angle_restrict_symm_err_e' and `gsl_sf_angle_restrict_pos_err_e' (why not mentioned in the manual?) are all wrong. The fact that the result happens to be very small implies by no means that the error is equally small (i.e., proportional to the result). A reasonable estimate would be result>err = GSL_DBL_EPSILON * fabs(theta  result>val) provided an argument already in the allowed range is not modified at all. II.2. The not so obvious part of problem B (above) is, that with a user specified error handler and sufficiently large arguments (i.e., >> 1.0e+15 in absolute value), the functions `gsl_sf_angle_restrict_symm' & Co may well exceed the allowed range (i.e., (\pi,\pi] for `gsl_sf_angle_restrict_symm*' resp. [0, 2\pi) for `gsl_sf_angle_restrict_pos*') by large amounts in either direction. That is a violation of the specification, at least in principle, and requires a somewhat more extensive modification.  The second, alternative, attached `gsl1.8mod2patch' addresses all these problems and ensures not to modify an argument already in the allowed range. ==================== The definition of an error for discontinuous functions such as *_symm* and *_pos* is of course problematic. One might argue that for arguments close to one of the discontinuities N*Pi (with N odd for *_symm* or N even for *_pos*), the error can be as large as TwoPi, but can never exceed TwoPi). This would yield the following code fragments for gsl_sf_angle_restrict_pos_err_e: result>val = r; e = GSL_DBL_EPSILON*fabs(thetar); if(fabs(r  M_PI) + e >= M_PI) result>err = TwoPi; else result>err = e; if(e > 0.0625) { GSL_ERROR ("error", GSL_ELOSS); } else return GSL_SUCCESS; and for gsl_sf_angle_restrict_symm_e: result>val = r; e = GSL_DBL_EPSILON*fabs(thetar); if(fabs(r) + e >= M_PI) result>err = TwoPi; else result>err = e; if(e > 0.0625) { GSL_ERROR ("error", GSL_ELOSS); } else return GSL_SUCCESS;  A third, alternative, `gsl1.8mod3patch' implements this notion of errors. ==================== In addition there is small, insignificant typo that you might want to correct as per the attached `patch01typo'. ==================== with kind regards, Peter Breitenlohner <address@hidden>
gslbug.c
Description: Text document
gsl1.8mod1patch
Description: Text document
gsl1.8mod2patch
Description: Text document
gsl1.8mod3patch
Description: Text document
patch01typo
Description: Text document
[Prev in Thread]  Current Thread  [Next in Thread] 