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