[Bug-gsl] A recommend repair of gsl_sf_angle_restrict_symm_err_e

[易昕]

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!