bug-gsl
[Top][All Lists]

## [Bug-gsl] Bug in gsl-1.8 trigonometric restriction functions + patches

 From: Peter Breitenlohner Subject: [Bug-gsl] Bug in gsl-1.8 trigonometric restriction functions + patches Date: Mon, 11 Sep 2006 16:42:11 +0200 (CEST)

Hi,

while looking at the GSL (gsl-1.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 gsl-1.8-mod1-patch' 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 gsl-1.8-mod2-patch' 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(theta-r);

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(theta-r);

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, gsl-1.8-mod3-patch' implements this notion of errors.

====================

In addition there is small, insignificant typo that you might want to
correct as per the attached patch-01-typo'.

====================

with kind regards,
Peter Breitenlohner <address@hidden>`

gslbug.c
Description: Text document

gsl-1.8-mod1-patch
Description: Text document

gsl-1.8-mod2-patch
Description: Text document

gsl-1.8-mod3-patch
Description: Text document

patch-01-typo
Description: Text document