[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[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!
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] A recommend repair of gsl_sf_angle_restrict_symm_err_e,
易昕 <=