diff options
| -rw-r--r-- | ChangeLog | 14 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/k_cosl.c | 20 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/k_sincosl.c | 21 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/k_sinl.c | 20 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_ceill.c | 31 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_cprojl.c | 54 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c | 87 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_ctanl.c | 88 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_floorl.c | 21 |
9 files changed, 328 insertions, 28 deletions
diff --git a/ChangeLog b/ChangeLog index f8239cc861..7a74feac63 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,17 @@ +2006-01-30 Steven Munroe <sjmunroe@us.ibm.com> + + * sysdeps/ieee754/ldbl-128ibm/k_cosl.c (__kernel_cosl): Correct index + for __sincosl_table. + * sysdeps/ieee754/ldbl-128ibm/k_sincosl.c (__kernel_sincosl): Likewise. + * sysdeps/ieee754/ldbl-128ibm/k_sinl.c (__kernel_sinl): Likewise. + + * sysdeps/ieee754/ldbl-128ibm/s_ceill.c: Correct sign of 0.0. + * sysdeps/ieee754/ldbl-128ibm/s_floorl.c: Likewise. + + * sysdeps/ieee754/ldbl-128ibm/s_cprojl.c: New file. + * sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c: New file. + * sysdeps/ieee754/ldbl-128ibm/s_ctanl.c: New file. + 2006-01-29 Daniel Jacobowitz <dan@codesourcery.com> * sysdeps/unix/sysv/linux/mips/bits/errno.h (ECANCELED, diff --git a/sysdeps/ieee754/ldbl-128ibm/k_cosl.c b/sysdeps/ieee754/ldbl-128ibm/k_cosl.c index 3baf8b7b66..b442582b3f 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_cosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_cosl.c @@ -104,6 +104,24 @@ __kernel_cosl(long double x, long double y) pre-computed tables, compute cosl(l) and sinl(l) using a Chebyshev polynomial of degree 10(11) and compute cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l). */ + int six = tix; + tix = ((six - 0x3ff00000) >> 4) + 0x3fff0000; + index = 0x3ffe - (tix >> 16); + hix = (tix + (0x200 << index)) & (0xfffffc00 << index); + x = fabsl (x); + switch (index) + { + case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break; + case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break; + default: + case 2: index = (hix - 0x3ffc3000) >> 10; break; + } + hix = (hix << 4) & 0x3fffffff; +/* + The following should work for double but generates the wrong index. + For now the code above converts double to ieee extended to compute + the index back to double for the h value. + index = 0x3fe - (tix >> 20); hix = (tix + (0x200 << index)) & (0xfffffc00 << index); x = fabsl (x); @@ -114,7 +132,7 @@ __kernel_cosl(long double x, long double y) default: case 2: index = (hix - 0x3fc30000) >> 14; break; } - +*/ SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0); l = y - (h - x); z = l * l; diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c index 1cea6fef6e..cd2ce7ad1d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_sincosl.c @@ -131,6 +131,25 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c Chebyshev polynomial of degree 10(11) and compute sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l) and cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l). */ + int six = tix; + tix = ((six - 0x3ff00000) >> 4) + 0x3fff0000; + index = 0x3ffe - (tix >> 16); + hix = (tix + (0x200 << index)) & (0xfffffc00 << index); + x = fabsl (x); + switch (index) + { + case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break; + case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break; + default: + case 2: index = (hix - 0x3ffc3000) >> 10; break; + } + hix = (hix << 4) & 0x3fffffff; +/* + The following should work for double but generates the wrong index. + For now the code above converts double to ieee extended to compute + the index back to double for the h value. + + index = 0x3fe - (tix >> 20); hix = (tix + (0x2000 << index)) & (0xffffc000 << index); x = fabsl (x); @@ -141,7 +160,7 @@ __kernel_sincosl(long double x, long double y, long double *sinx, long double *c default: case 2: index = (hix - 0x3fc30000) >> 14; break; } - +*/ SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0); if (iy) l = y - (h - x); diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c index c4dce26f9a..24cb551b6e 100644 --- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c +++ b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c @@ -104,6 +104,24 @@ __kernel_sinl(long double x, long double y, int iy) 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). */ + int six = tix; + tix = ((six - 0x3ff00000) >> 4) + 0x3fff0000; + index = 0x3ffe - (tix >> 16); + hix = (tix + (0x200 << index)) & (0xfffffc00 << index); + x = fabsl (x); + switch (index) + { + case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break; + case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break; + default: + case 2: index = (hix - 0x3ffc3000) >> 10; break; + } + hix = (hix << 4) & 0x3fffffff; +/* + The following should work for double but generates the wrong index. + For now the code above converts double to ieee extended to compute + the index back to double for the h value. + index = 0x3fe - (tix >> 20); hix = (tix + (0x2000 << index)) & (0xffffc000 << index); x = fabsl (x); @@ -114,7 +132,7 @@ __kernel_sinl(long double x, long double y, int iy) default: case 2: index = (hix - 0x3fc30000) >> 14; break; } - +*/ SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0); if (iy) l = y - (h - x); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c index a4cbbe1048..a606548873 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_ceill.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_ceill.c @@ -18,9 +18,6 @@ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. */ -/* This has been coded in assembler because GCC makes such a mess of it - when it's coded in C. */ - #include <math.h> #include <fenv.h> #include <math_ldbl_opt.h> @@ -44,19 +41,23 @@ __ceill (x) u.d = x; if (fabs (u.dd[0]) < TWO52) - { + { + double high = u.dd[0]; fesetround(FE_UPWARD); - if (u.dd[0] > 0.0) + if (high > 0.0) { - u.dd[0] += TWO52; - u.dd[0] -= TWO52; + high += TWO52; + high -= TWO52; + if (high == -0.0) high = 0.0; } - else if (u.dd[0] < 0.0) + else if (high < 0.0) { - u.dd[0] -= TWO52; - u.dd[0] += TWO52; + high -= TWO52; + high += TWO52; + if (high == 0.0) high = -0.0; } - u.dd[1] = 0.0; + u.dd[0] = high; + u.dd[1] = 0.0; fesetround(mode); } else if (fabs (u.dd[1]) < TWO52 && u.dd[1] != 0.0) @@ -82,10 +83,10 @@ __ceill (x) adjust for that. */ high = nextafter (u.dd[0], 0.0); low = u.dd[1] + (u.dd[0] - high); - } + } fesetround(FE_UPWARD); low += TWO52; - low -= TWO52; + low -= TWO52; fesetround(mode); } else if (u.dd[0] < 0.0) @@ -103,10 +104,10 @@ __ceill (x) adjust for that. */ high = nextafter (u.dd[0], 0.0); low = u.dd[1] + (u.dd[0] - high); - } + } fesetround(FE_UPWARD); low -= TWO52; - low += TWO52; + low += TWO52; fesetround(mode); } u.dd[0] = high + low; diff --git a/sysdeps/ieee754/ldbl-128ibm/s_cprojl.c b/sysdeps/ieee754/ldbl-128ibm/s_cprojl.c new file mode 100644 index 0000000000..2167db3d91 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/s_cprojl.c @@ -0,0 +1,54 @@ +/* Compute projection of complex long double value to Riemann sphere. + Copyright (C) 1997,1999,2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <complex.h> +#include <math.h> +#include <math_ldbl_opt.h> + +__complex__ long double +__cprojl (__complex__ long double x) +{ + __complex__ long double res; + + if (isnan (__real__ x) && isnan (__imag__ x)) + return x; + else if (!isfinite (__real__ x) || !isfinite (__imag__ x)) + { + __real__ res = INFINITY; + __imag__ res = __copysignl (0.0, __imag__ x); + } + else + { + long double den = (__real__ x * __real__ x + __imag__ x * __imag__ x + + 1.0); + + __real__ res = (2.0 * __real__ x) / den; + __imag__ res = (2.0 * __imag__ x) / den; + /* __gcc_qmul does not respect -0.0 so we need the following fixup. */ + if (__real__ x == 0.0) + __real__ res = __real__ x; + + if (__imag__ x == 0.0) + __imag__ res = __imag__ x; + } + + return res; +} +long_double_symbol (libm, __cprojl, cprojl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c new file mode 100644 index 0000000000..7d619039ca --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c @@ -0,0 +1,87 @@ +/* Complex hyperbole tangent for long double. IBM extended format version. + Copyright (C) 1997,2005,2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <complex.h> +#include <fenv.h> +#include <math.h> +#include <math_ldbl_opt.h> + +#include "math_private.h" + + +__complex__ long double +__ctanhl (__complex__ long double x) +{ + __complex__ long double res; + + if (!isfinite (__real__ x) || !isfinite (__imag__ x)) + { + if (__isinfl (__real__ x)) + { + __real__ res = __copysignl (1.0, __real__ x); + __imag__ res = __copysignl (0.0, __imag__ x); + } + else if (__imag__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + +#ifdef FE_INVALID + if (__isinfl (__imag__ x)) + feraiseexcept (FE_INVALID); +#endif + } + } + else + { + long double sin2ix, cos2ix; + long double den; + + __sincosl (2.0 * __imag__ x, &sin2ix, &cos2ix); + + den = (__ieee754_coshl (2.0 * __real__ x) + cos2ix); + + if (den == 0.0L) + { + __complex__ long double ez = __cexpl (x); + __complex__ long double emz = __cexpl (-x); + + res = (ez - emz) / (ez + emz); + } + else + { + __real__ res = __ieee754_sinhl (2.0 * __real__ x) / den; + __imag__ res = sin2ix / den; + } + /* __gcc_qmul does not respect -0.0 so we need the following fixup. */ + if ((__real__ res == 0.0) && (__real__ x == 0.0)) + __real__ res = __real__ x; + + if ((__real__ res == 0.0) && (__imag__ x == 0.0)) + __imag__ res = __imag__ x; + } + + return res; +} +long_double_symbol (libm, __ctanhl, ctanhl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c new file mode 100644 index 0000000000..00d2aa6b6e --- /dev/null +++ b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c @@ -0,0 +1,88 @@ +/* Complex tangent function for long double. IBM extended format version. + Copyright (C) 1997,2005,2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <complex.h> +#include <fenv.h> +#include <math.h> +#include <math_ldbl_opt.h> + +#include "math_private.h" + + +__complex__ long double +__ctanl (__complex__ long double x) +{ + __complex__ long double res; + + if (!isfinite (__real__ x) || !isfinite (__imag__ x)) + { + if (__isinfl (__imag__ x)) + { + __real__ res = __copysignl (0.0, __real__ x); + __imag__ res = __copysignl (1.0, __imag__ x); + } + else if (__real__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + +#ifdef FE_INVALID + if (__isinfl (__real__ x)) + feraiseexcept (FE_INVALID); +#endif + } + } + else + { + long double sin2rx, cos2rx; + long double den; + + __sincosl (2.0 * __real__ x, &sin2rx, &cos2rx); + + den = cos2rx + __ieee754_coshl (2.0 * __imag__ x); + + + if (den == 0.0) + { + __complex__ long double ez = __cexpl (1.0i * x); + __complex__ long double emz = __cexpl (-1.0i * x); + + res = (ez - emz) / (ez + emz) * -1.0i; + } + else + { + __real__ res = sin2rx / den; + __imag__ res = __ieee754_sinhl (2.0 * __imag__ x) / den; + } + /* __gcc_qmul does not respect -0.0 so we need the following fixup. */ + if ((__real__ res == 0.0) && (__real__ x == 0.0)) + __real__ res = __real__ x; + + if ((__real__ res == 0.0) && (__imag__ x == 0.0)) + __imag__ res = __imag__ x; + } + + return res; +} +long_double_symbol (libm, __ctanl, ctanl); diff --git a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c index e8ef6e849d..2be5e28698 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_floorl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_floorl.c @@ -18,9 +18,6 @@ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA. */ -/* This has been coded in assembler because GCC makes such a mess of it - when it's coded in C. */ - #include <math.h> #include <fenv.h> #include <math_ldbl_opt.h> @@ -44,18 +41,22 @@ __floorl (x) u.d = x; if (fabs (u.dd[0]) < TWO52) - { + { + double high = u.dd[0]; fesetround(FE_DOWNWARD); - if (u.dd[0] > 0.0) + if (high > 0.0) { - u.dd[0] += TWO52; - u.dd[0] -= TWO52; + high += TWO52; + high -= TWO52; + if (high == -0.0) high = 0.0; } - else if (u.dd[0] < 0.0) + else if (high < 0.0) { - u.dd[0] -= TWO52; - u.dd[0] += TWO52; + high -= TWO52; + high += TWO52; + if (high == 0.0) high = -0.0; } + u.dd[0] = high; u.dd[1] = 0.0; fesetround(mode); } |