Mercurial > hg > octave-kai > gnulib-hg
changeset 12994:d63de5fca804
Fix values returned by sinl, cosl.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Sun, 14 Mar 2010 23:19:29 +0100 |
parents | e6ab76eb5d57 |
children | 6991d7777e06 |
files | ChangeLog lib/sincosl.c lib/trigl.h |
diffstat | 3 files changed, 43 insertions(+), 15 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2010-03-14 Bruno Haible <bruno@clisp.org> + + 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. + 2010-03-14 Bruno Haible <bruno@clisp.org> Fix compilation error when modules 'posix_spawn[p]' are not used.
--- a/lib/sincosl.c +++ b/lib/sincosl.c @@ -136,11 +136,12 @@ 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 + 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(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */ + 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,11 +159,14 @@ 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_SIN_LO] + - (sincosl_table[index + SINCOSL_SIN_HI] * cos_l_m1) + - (sincosl_table[index + SINCOSL_COS_HI] * sin_l)); + 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,11 +199,12 @@ 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 + 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(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */ + 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,10 +218,14 @@ 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)); + /* 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; } }
--- a/lib/trigl.h +++ b/lib/trigl.h @@ -18,7 +18,18 @@ 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); -