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 <[EMAIL PROTECTED]>/*
* gslbug.c
* compile: gcc -O2 -Wall -Wextra -o gslbug gslbug.c -lgsl -lgslcblas
*
* A small demo for some bugs in the GSL (V1.8)
* "Trigonometric Restriction" routines
*/
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_errno.h>
void
my_handler(const char *reason,
const char *file,
int line,
int gsl_errno) {
printf("GSL error %d (%s:%d %s) ignored\n",
gsl_errno, file, line, reason);
}
int
main(void) {
int i;
double val, res;
gsl_sf_result result;
gsl_error_handler_t *old_handler;
const double TwoPi = 2*M_PI;
printf("\n\tDemonstrate bugs in the GSL(-1.8) routines\n"
"\t\tgsl_sf_angle_restrict_symm and gsl_sf_angle_restrict_pos\n");
old_handler = gsl_set_error_handler(my_handler);
/*
* Part I.A. Missing check for huge negative arguments
*/
printf("\n\tPart I.A. Missing check for huge negative arguments:\n");
val = 1.0e+14;
for (i=0; i<2; i++) {
printf("gsl_sf_angle_restrict_symm(%g) => ", val);
printf("%g\n", gsl_sf_angle_restrict_symm(val));
printf("gsl_sf_angle_restrict_pos(%g) => ", val);
printf("%g\n", gsl_sf_angle_restrict_pos(val));
printf("gsl_sf_angle_restrict_symm(%g) => ", -val);
printf("%g\n", gsl_sf_angle_restrict_symm(-val));
printf("gsl_sf_angle_restrict_pos(%g) => ", -val);
printf("%g\n", gsl_sf_angle_restrict_pos(-val));
val *= 10.0;
}
/*
* Part I.B. gsl_sf_angle_restrict_pos(x)<0
*/
printf("\n\tPart I.B. gsl_sf_angle_restrict_pos(x)<0:\n");
for (i=-3; i<6; i+=3)
printf("gsl_sf_angle_restrict_pos(TwoPi%+d*DBL_EPSILON) = %.17g\n",
i, gsl_sf_angle_restrict_pos(TwoPi + i*DBL_EPSILON));
/*
* Part II.1. Wrong error estimates
*/
printf("\n\tPart II.1. Wrong error estimates:\n");
val = 1.0e-13;
for (i=0; i<3; i++) {
gsl_sf_angle_restrict_symm_err_e(val, &result);
printf("gsl_sf_angle_restrict_symm_err_e(%.15e) => %.15e +- %g\n",
val, result.val, result.err);
gsl_sf_angle_restrict_symm_err_e(-val, &result);
printf("gsl_sf_angle_restrict_symm_err_e(%.15e) => %.15e +- %g\n",
-val, result.val, result.err);
gsl_sf_angle_restrict_pos_err_e(val, &result);
printf("gsl_sf_angle_restrict_pos_err_e(%.15e) => %.15e +- %g\n",
val, result.val, result.err);
gsl_sf_angle_restrict_pos_err_e(-val, &result);
printf("gsl_sf_angle_restrict_pos_err_e(%.15e) => %.15e +- %g\n",
-val, result.val, result.err);
val /= 10.0;
}
printf("\n\t More error estimates:\n");
for (i=-3; i<9; i+=3) {
val = M_PI + i*DBL_EPSILON;
gsl_sf_angle_restrict_symm_err_e(val, &result);
printf("gsl_sf_angle_restrict_symm_err_e(Pi%+d*DBL_EPSILON) => %.16f +-
%g\n",
i, result.val, result.err);
}
printf("\n");
for (i=-3; i<9; i+=3) {
val = TwoPi + i*DBL_EPSILON;
gsl_sf_angle_restrict_pos_err_e(val, &result);
printf("gsl_sf_angle_restrict_pos_err_e(TwoPi%+d*DBL_EPSILON) => %.16f +-
%g\n",
i, result.val, result.err);
}
/*
* Part II.2. Huge (positive or negative) results
*/
printf("\n\tPart II.2. Huge (positive or negative) results:\n");
val = +1.0e+300;
do {
res = gsl_sf_angle_restrict_pos(val);
printf("gsl_sf_angle_restrict_pos(%g) = %g\n", val, res);
val = res;
} while (fabs(val) > 1.0e+200);
printf("\n");
return 0;
}
diff -ur gsl-1.8.orig/specfunc/trig.c gsl-1.8/specfunc/trig.c
--- gsl-1.8.orig/specfunc/trig.c 2005-06-26 15:25:35.000000000 +0200
+++ gsl-1.8/specfunc/trig.c 2006-08-21 22:05:31.000000000 +0200
@@ -544,11 +544,11 @@
if(r > M_PI) r -= TwoPi;
result->val = r;
- if(theta > 0.0625/GSL_DBL_EPSILON) {
+ if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
result->err = fabs(result->val);
GSL_ERROR ("error", GSL_ELOSS);
}
- else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
+ else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
}
@@ -568,14 +568,16 @@
const double TwoPi = 2*(P1 + P2 + P3);
const double y = 2*floor(theta/TwoPi);
+ double r = ((theta - y*P1) - y*P2) - y*P3;
- result->val = ((theta - y*P1) - y*P2) - y*P3;
+ if (r < 0.0) r += TwoPi; /* may happen due to FP rounding */
+ result->val = r;
- if(theta > 0.0625/GSL_DBL_EPSILON) {
+ if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
result->err = fabs(result->val);
GSL_ERROR ("error", GSL_ELOSS);
}
- else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
+ else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
return GSL_SUCCESS;
}
diff -ur gsl-1.8.orig/specfunc/trig.c gsl-1.8/specfunc/trig.c
--- gsl-1.8.orig/specfunc/trig.c 2005-06-26 15:25:35.000000000 +0200
+++ gsl-1.8/specfunc/trig.c 2006-08-23 14:01:08.000000000 +0200
@@ -538,24 +538,25 @@
const double P3 = 4 * 2.6951514290790594840552e-15;
const double TwoPi = 2*(P1 + P2 + P3);
- const double y = 2*floor(theta/TwoPi);
- double r = ((theta - y*P1) - y*P2) - y*P3;
+ double r = theta;
+ /* dont modify arg if already in allowed range */
+ while(fabs(r) > M_PI) {
+ const double y = 2*ceil(r/TwoPi - 0.5);
+ r = ((r - y*P1) - y*P2) - y*P3;
+ /* for |theta|=1.0e+300, we might get |r|>1.0e+283 */
+ }
+ /* now |r| < Pi + epsilon */
if(r > M_PI) r -= TwoPi;
+ else if(r <= -M_PI) r += TwoPi;
+
result->val = r;
+ result->err = GSL_DBL_EPSILON*fabs(theta-r);
- if(theta > 0.0625/GSL_DBL_EPSILON) {
- result->err = fabs(result->val);
+ if(result->err > 0.0625) {
GSL_ERROR ("error", GSL_ELOSS);
}
- else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
- result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
- else {
- result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
+ else return GSL_SUCCESS;
}
@@ -567,22 +568,25 @@
const double P3 = 4 * 2.69515142907905952645e-15;
const double TwoPi = 2*(P1 + P2 + P3);
- const double y = 2*floor(theta/TwoPi);
+ double r = theta;
- result->val = ((theta - y*P1) - y*P2) - y*P3;
+ /* dont modify arg if already in allowed range */
+ while(fabs(r - M_PI) > M_PI) {
+ const double y = 2*floor(r/TwoPi);
+ r = ((r - y*P1) - y*P2) - y*P3;
+ /* for |theta|=1.0e+300, we might get |r|>1.0e+283 */
+ }
+ /* now |r - Pi| < Pi + epsilon */
+ if(r >= TwoPi) r -= TwoPi;
+ else if(r < 0.0) r += TwoPi;
+
+ result->val = r;
+ result->err = GSL_DBL_EPSILON*fabs(theta-r);
- if(theta > 0.0625/GSL_DBL_EPSILON) {
- result->err = fabs(result->val);
+ if(result->err > 0.0625) {
GSL_ERROR ("error", GSL_ELOSS);
}
- else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
- result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
- else {
- result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
+ else return GSL_SUCCESS;
}
diff -ur gsl-1.8.orig/specfunc/trig.c gsl-1.8/specfunc/trig.c
--- gsl-1.8.orig/specfunc/trig.c 2005-06-26 15:25:35.000000000 +0200
+++ gsl-1.8/specfunc/trig.c 2006-08-23 15:27:05.000000000 +0200
@@ -538,24 +538,28 @@
const double P3 = 4 * 2.6951514290790594840552e-15;
const double TwoPi = 2*(P1 + P2 + P3);
- const double y = 2*floor(theta/TwoPi);
- double r = ((theta - y*P1) - y*P2) - y*P3;
+ double r = theta, e;
+ /* dont modify arg if already in allowed range */
+ while(fabs(r) > M_PI) {
+ const double y = 2*ceil(r/TwoPi - 0.5);
+ r = ((r - y*P1) - y*P2) - y*P3;
+ /* for |theta|=1.0e+300, we might get |r|>1.0e+283 */
+ }
+ /* now |r| < Pi + epsilon */
if(r > M_PI) r -= TwoPi;
+ else if(r <= -M_PI) r += TwoPi;
+
result->val = r;
+ e = fabs(theta-r)*(GSL_DBL_EPSILON/2);
+
+ if(fabs(r) + e >= M_PI) result->err = TwoPi;
+ else result->err = e;
- if(theta > 0.0625/GSL_DBL_EPSILON) {
- result->err = fabs(result->val);
+ if(e > 0.03125) {
GSL_ERROR ("error", GSL_ELOSS);
}
- else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
- result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
- else {
- result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
+ else return GSL_SUCCESS;
}
@@ -567,22 +571,28 @@
const double P3 = 4 * 2.69515142907905952645e-15;
const double TwoPi = 2*(P1 + P2 + P3);
- const double y = 2*floor(theta/TwoPi);
+ double r = theta, e;
- result->val = ((theta - y*P1) - y*P2) - y*P3;
+ /* dont modify arg if already in allowed range */
+ while(fabs(r - M_PI) > M_PI) {
+ const double y = 2*floor(r/TwoPi);
+ r = ((r - y*P1) - y*P2) - y*P3;
+ /* for |theta|=1.0e+300, we might get |r|>1.0e+283 */
+ }
+ /* now |r - Pi| < Pi + epsilon */
+ if(r >= TwoPi) r -= TwoPi;
+ else if(r < 0.0) r += TwoPi;
+
+ result->val = r;
+ e = fabs(theta-r)*(GSL_DBL_EPSILON/2);
+
+ if(fabs(r - M_PI) + e >= M_PI) result->err = TwoPi;
+ else result->err = e;
- if(theta > 0.0625/GSL_DBL_EPSILON) {
- result->err = fabs(result->val);
+ if(e > 0.03125) {
GSL_ERROR ("error", GSL_ELOSS);
}
- else if(theta > 0.0625/GSL_SQRT_DBL_EPSILON) {
- result->err = GSL_SQRT_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
- else {
- result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
- return GSL_SUCCESS;
- }
+ else return GSL_SUCCESS;
}
diff -ur gsl-1.8.orig/acconfig.h gsl-1.8/acconfig.h
--- gsl-1.8.orig/acconfig.h 2006-03-17 16:53:51.000000000 +0100
+++ gsl-1.8/acconfig.h 2006-03-17 16:53:51.000000000 +0100
@@ -39,10 +39,10 @@
#undef HAVE_OPENBSD_IEEE_INTERFACE
#undef HAVE_DARWIN_IEEE_INTERFACE
-/* Define this is IEEE comparisons work correctly (e.g. NaN != NaN) */
+/* Define this if IEEE comparisons work correctly (e.g. NaN != NaN) */
#undef HAVE_IEEE_COMPARISONS
-/* Define this is IEEE denormalized numbers are available */
+/* Define this if IEEE denormalized numbers are available */
#undef HAVE_IEEE_DENORMALS
/* Define a rounding function which moves extended precision values
diff -ur gsl-1.8.orig/config.h.in gsl-1.8/config.h.in
--- gsl-1.8.orig/config.h.in 2006-03-31 20:24:31.000000000 +0200
+++ gsl-1.8/config.h.in 2006-03-31 20:24:31.000000000 +0200
@@ -189,10 +189,10 @@
#undef HAVE_OPENBSD_IEEE_INTERFACE
#undef HAVE_DARWIN_IEEE_INTERFACE
-/* Define this is IEEE comparisons work correctly (e.g. NaN != NaN) */
+/* Define this if IEEE comparisons work correctly (e.g. NaN != NaN) */
#undef HAVE_IEEE_COMPARISONS
-/* Define this is IEEE denormalized numbers are available */
+/* Define this if IEEE denormalized numbers are available */
#undef HAVE_IEEE_DENORMALS
/* Define a rounding function which moves extended precision values
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
