bug-gsl
[Top][All Lists]

## [Bug-gsl] A recommend repair of gsl_sf_angle_restrict_symm_err_e

 From: 易昕 Subject: [Bug-gsl] A recommend repair of gsl_sf_angle_restrict_symm_err_e Date: Wed, 26 Jun 2019 10:20:03 +0800

```The program try to reduce input into the [-pi,pi]，I found inaccuracy when
the input is slight smaller than -2*k*pi.

For example, with the input -226.19467105845234

gsl_sf_angle_restrict_symm_err_e(-226.19467105845234) = 1.27704742478e-11
while I use the higher precision the right result should be
1.27701649912052602689124836895e-11
the relative error is 2.421711882e-5

A recommend repair is:
replace:
if(r >  M_PI) { r = (((r-2*P1)-2*P2)-2*P3); }  /* r-TwoPi */
else if (r < -M_PI) r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */

with:
if(r >  M_PI) { r = (((theta-(y+2)*P1)-(y+2)*P2)-(y+2)*P3); }  /* r-TwoPi */
else if (r < -M_PI) r = (((theta+(2-y)*P1)+(2-y)*P2)+(2-y)*P3); /*
r+TwoPi */

In this way to avoid the inaccuracy introduced in the previous
sentence: double r = ((theta - y*P1) - y*P2) - y*P3;

PS:
int gsl_sf_angle_restrict_symm_err_e(const double theta, gsl_sf_result *
result)
{
/* synthetic extended precision constants */
const double P1 = 4 * 7.8539812564849853515625e-01;
const double P2 = 4 * 3.7748947079307981766760e-08;
const double P3 = 4 * 2.6951514290790594840552e-15;
const double TwoPi = 2*(P1 + P2 + P3);

const double y = GSL_SIGN(theta) * 2 * floor(fabs(theta)/TwoPi);
double r = ((theta - y*P1) - y*P2) - y*P3;

if(r >  M_PI) { r = (((r-2*P1)-2*P2)-2*P3); }  /* r-TwoPi */
else if (r < -M_PI) r = (((r+2*P1)+2*P2)+2*P3); /* r+TwoPi */

*/* To avoid round error in the * double r = ((theta - y*P1) - y*P2) - y*P3;

*if(r >  M_PI) { r = (((theta-(y+2)*P1)-(y+2)*P2)-(y+2)*P3); }  /* r-TwoPi
*/  else if (r < -M_PI) r = (((theta+(2-y)*P1)+(2-y)*P2)+(2-y)*P3); /*
r+TwoPi */*
*}*

**/*
result->val = r;

if(fabs(theta) > 0.0625/GSL_DBL_EPSILON) {
result->val = GSL_NAN;
result->err = GSL_NAN;
GSL_ERROR ("error", GSL_ELOSS);
}
else if(fabs(theta) > 0.0625/GSL_SQRT_DBL_EPSILON) {
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val - theta);
return GSL_SUCCESS;
}
else {
double delta = fabs(result->val - theta);
result->err = 2.0 * GSL_DBL_EPSILON * ((delta < M_PI) ? delta : M_PI);
return GSL_SUCCESS;
}
}

Xin Yi
Best wishes!

```