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