[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: cosl test failure
From: |
Bruno Haible |
Subject: |
Re: cosl test failure |
Date: |
Sun, 14 Mar 2010 23:18:31 +0100 |
User-agent: |
KMail/1.9.9 |
Hi Paolo,
On 2010-01-18, in
<http://lists.gnu.org/archive/html/bug-gnulib/2010-01/msg00256.html>,
I observed that gnulib's cosl replacement function does not have the
necessary accuracy. This was due to a wrong formula: The term
sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1
was being subtracted, when it should have been added.
Also, in this formula and the other one for sinl(), you omitted one of
the five summands; I cannot see a reason why.
This fixes it.
2010-03-14 Bruno Haible <address@hidden>
Fix values returned by sinl, cosl.
* lib/trigl.h: Add specification comments.
* lib/sincosl.c (kernel_sinl, kernel_cosl): Fix comments and formula
that combines the values from the precomputed table with the values of
the Chebyshev polynomials.
*** lib/trigl.h.orig Sun Mar 14 23:11:58 2010
--- lib/trigl.h Sun Mar 14 22:14:54 2010
***************
*** 18,24 ****
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. */
extern int ieee754_rem_pio2l (long double x, long double *y);
extern long double kernel_sinl (long double x, long double y, int iy);
- extern long double kernel_cosl (long double x, long double y);
--- 18,35 ----
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>. */
+ /* Decompose x into x = k * π/2 + r
+ where k is an integer and abs(r) <= π/4.
+ Store r in y[0] and y[1] (main part in y[0], small additional part in
+ y[1], r = y[0] + y[1]).
+ Return k. */
extern int ieee754_rem_pio2l (long double x, long double *y);
+
+ /* Compute and return sinl (x + y), where x is the main part and y is the
+ small additional part of a floating-point number.
+ iy is 0 when y is known to be 0.0, otherwise iy is 1. */
extern long double kernel_sinl (long double x, long double y, int iy);
+ /* Compute and return cosl (x + y), where x is the main part and y is the
+ small additional part of a floating-point number. */
+ extern long double kernel_cosl (long double x, long double y);
*** lib/sincosl.c.orig Sun Mar 14 23:11:58 2010
--- lib/sincosl.c Sun Mar 14 23:04:53 2010
***************
*** 136,146 ****
else
{
/* So that we don't have to use too large polynomial, we find
! l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83
! possible values for h. We look up cosl(h) and sinl(h) in
pre-computed tables, compute cosl(l) and sinl(l) using a
Chebyshev polynomial of degree 10(11) and compute
! sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */
x -= 0.1484375L;
index = (int) (x * 128L + 0.5L);
h = index / 128.0L;
--- 136,147 ----
else
{
/* So that we don't have to use too large polynomial, we find
! k and l such that x = k + l, where fabsl(l) <= 1.0/256 with 83
! possible values for k. We look up cosl(k) and sinl(k) in
pre-computed tables, compute cosl(l) and sinl(l) using a
Chebyshev polynomial of degree 10(11) and compute
! sinl(k+l) = sinl(k)cosl(l) + cosl(k)sinl(l).
! Furthermore write k = 0.1484375 + h. */
x -= 0.1484375L;
index = (int) (x * 128L + 0.5L);
h = index / 128.0L;
***************
*** 158,168 ****
z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
index *= 4;
z =
! sincosl_table[index + SINCOSL_SIN_HI] +
! (sincosl_table[index + SINCOSL_SIN_LO] +
! (sincosl_table[index + SINCOSL_SIN_HI] * cos_l_m1) +
! (sincosl_table[index + SINCOSL_COS_HI] * sin_l));
return z * sign;
}
}
--- 159,172 ----
z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
index *= 4;
+ /* We rely on this expression not being "contracted" by the compiler
+ (cf. ISO C 99 section 6.5 paragraph 8). */
z =
! sincosl_table[index + SINCOSL_SIN_HI]
! + (sincosl_table[index + SINCOSL_COS_HI] * sin_l
! + (sincosl_table[index + SINCOSL_SIN_HI] * cos_l_m1
! + (sincosl_table[index + SINCOSL_SIN_LO] * (1 + cos_l_m1)
! + sincosl_table[index + SINCOSL_COS_LO] * sin_l)));
return z * sign;
}
}
***************
*** 195,205 ****
else
{
/* So that we don't have to use too large polynomial, we find
! l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83
! possible values for h. We look up cosl(h) and sinl(h) in
pre-computed tables, compute cosl(l) and sinl(l) using a
Chebyshev polynomial of degree 10(11) and compute
! sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */
x -= 0.1484375L;
index = (int) (x * 128L + 0.5L);
h = index / 128.0L;
--- 199,210 ----
else
{
/* So that we don't have to use too large polynomial, we find
! k and l such that x = k + l, where fabsl(l) <= 1.0/256 with 83
! possible values for k. We look up cosl(k) and sinl(k) in
pre-computed tables, compute cosl(l) and sinl(l) using a
Chebyshev polynomial of degree 10(11) and compute
! cosl(k+l) = cosl(k)cosl(l) - sinl(k)sinl(l).
! Furthermore write k = 0.1484375 + h. */
x -= 0.1484375L;
index = (int) (x * 128L + 0.5L);
h = index / 128.0L;
***************
*** 213,222 ****
z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
index *= 4;
! z = sincosl_table [index + SINCOSL_COS_HI]
! + (sincosl_table [index + SINCOSL_COS_LO]
! - (sincosl_table [index + SINCOSL_SIN_HI] * sin_l)
! - (sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1));
return z;
}
}
--- 218,231 ----
z * (SCOS1 + z * (SCOS2 + z * (SCOS3 + z * (SCOS4 + z * SCOS5))));
index *= 4;
! /* We rely on this expression not being "contracted" by the compiler
! (cf. ISO C 99 section 6.5 paragraph 8). */
! z =
! sincosl_table [index + SINCOSL_COS_HI]
! - (sincosl_table [index + SINCOSL_SIN_HI] * sin_l
! - (sincosl_table [index + SINCOSL_COS_HI] * cos_l_m1
! + (sincosl_table [index + SINCOSL_COS_LO] * (1 + cos_l_m1)
! - sincosl_table [index + SINCOSL_SIN_LO] * sin_l)));
return z;
}
}
- Re: cosl test failure,
Bruno Haible <=