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

[Peter Breitenlohner]

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