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;
   }