diff options
Diffstat (limited to 'sysdeps/libm-ieee754')
28 files changed, 1396 insertions, 626 deletions
diff --git a/sysdeps/libm-ieee754/s_clog10.c b/sysdeps/libm-ieee754/s_clog10.c new file mode 100644 index 0000000000..5a9de75733 --- /dev/null +++ b/sysdeps/libm-ieee754/s_clog10.c @@ -0,0 +1,65 @@ +/* Compute complex base 10 logarithm. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. 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_private.h" + + +__complex__ double +__clog10 (__complex__ double x) +{ + __complex__ double result; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls == FP_ZERO && icls == FP_ZERO) + { + /* Real and imaginary part are 0.0. */ + __imag__ result = signbit (__real__ x) ? M_PI : 0.0; + __imag__ result = __copysign (__imag__ result, __imag__ x); + /* Yes, the following line raises an exception. */ + __real__ result = -1.0 / fabs (__real__ x); + } + else if (rcls != FP_NAN && icls != FP_NAN) + { + /* Neither real nor imaginary part is NaN. */ + __real__ result = __ieee754_log10 (__ieee754_hypot (__real__ x, + __imag__ x)); + __imag__ result = __ieee754_atan2 (__imag__ x, __real__ x); + } + else + { + __imag__ result = __nan (""); + if (rcls == FP_INFINITE || icls == FP_INFINITE) + /* Real or imaginary part is infinite. */ + __real__ result = HUGE_VAL; + else + __real__ result = __nan (""); + } + + return result; +} +weak_alias (__clog10, clog10) +#ifdef NO_LONG_DOUBLE +strong_alias (__clog10, __clog10l) +weak_alias (__clog10, clog10l) +#endif diff --git a/sysdeps/libm-ieee754/s_clog10f.c b/sysdeps/libm-ieee754/s_clog10f.c new file mode 100644 index 0000000000..4e2fa83edf --- /dev/null +++ b/sysdeps/libm-ieee754/s_clog10f.c @@ -0,0 +1,61 @@ +/* Compute complex base 10 logarithm. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. 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_private.h" + + +__complex__ float +__clog10f (__complex__ float x) +{ + __complex__ float result; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls == FP_ZERO && icls == FP_ZERO) + { + /* Real and imaginary part are 0.0. */ + __imag__ result = signbit (__real__ x) ? M_PI : 0.0; + __imag__ result = __copysignf (__imag__ result, __imag__ x); + /* Yes, the following line raises an exception. */ + __real__ result = -1.0 / fabsf (__real__ x); + } + else if (rcls != FP_NAN && icls != FP_NAN) + { + /* Neither real nor imaginary part is NaN. */ + __real__ result = __ieee754_log10f (__ieee754_hypotf (__real__ x, + __imag__ x)); + __imag__ result = __ieee754_atan2f (__imag__ x, __real__ x); + } + else + { + __imag__ result = __nanf (""); + if (rcls == FP_INFINITE || icls == FP_INFINITE) + /* Real or imaginary part is infinite. */ + __real__ result = HUGE_VALF; + else + __real__ result = __nanf (""); + } + + return result; +} +weak_alias (__clog10f, clog10f) diff --git a/sysdeps/libm-ieee754/s_clog10l.c b/sysdeps/libm-ieee754/s_clog10l.c new file mode 100644 index 0000000000..bf7d394c42 --- /dev/null +++ b/sysdeps/libm-ieee754/s_clog10l.c @@ -0,0 +1,61 @@ +/* Compute complex natural logarithm. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. 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_private.h" + + +__complex__ long double +__clog10l (__complex__ long double x) +{ + __complex__ long double result; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls == FP_ZERO && icls == FP_ZERO) + { + /* Real and imaginary part are 0.0. */ + __imag__ result = signbit (__real__ x) ? M_PI : 0.0; + __imag__ result = __copysignl (__imag__ result, __imag__ x); + /* Yes, the following line raises an exception. */ + __real__ result = -1.0 / fabsl (__real__ x); + } + else if (rcls != FP_NAN && icls != FP_NAN) + { + /* Neither real nor imaginary part is NaN. */ + __real__ result = __ieee754_log10l (__ieee754_hypotl (__real__ x, + __imag__ x)); + __imag__ result = __ieee754_atan2l (__imag__ x, __real__ x); + } + else + { + __imag__ result = __nanl (""); + if (rcls == FP_INFINITE || icls == FP_INFINITE) + /* Real or imaginary part is infinite. */ + __real__ result = HUGE_VALL; + else + __real__ result = __nanl (""); + } + + return result; +} +weak_alias (__clog10l, clog10l) diff --git a/sysdeps/libm-ieee754/s_fma.c b/sysdeps/libm-ieee754/s_fma.c new file mode 100644 index 0000000000..53974dc40c --- /dev/null +++ b/sysdeps/libm-ieee754/s_fma.c @@ -0,0 +1,28 @@ +/* Compute x * y + z as ternary operation. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +double +__fma (double x, double y, double z) +{ + return (x * y) + z; +} +weak_alias (__fma, fma) diff --git a/sysdeps/libm-ieee754/s_fmaf.c b/sysdeps/libm-ieee754/s_fmaf.c new file mode 100644 index 0000000000..44e5e15def --- /dev/null +++ b/sysdeps/libm-ieee754/s_fmaf.c @@ -0,0 +1,28 @@ +/* Compute x * y + z as ternary operation. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +float +__fmaf (float x, float y, float z) +{ + return (x * y) + z; +} +weak_alias (__fmaf, fmaf) diff --git a/sysdeps/libm-ieee754/s_fmal.c b/sysdeps/libm-ieee754/s_fmal.c new file mode 100644 index 0000000000..f73751f4db --- /dev/null +++ b/sysdeps/libm-ieee754/s_fmal.c @@ -0,0 +1,28 @@ +/* Compute x * y + z as ternary operation. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +long double +__fmal (long double x, long double y, long double z) +{ + return (x * y) + z; +} +weak_alias (__fmal, fmal) diff --git a/sysdeps/libm-ieee754/s_llrint.c b/sysdeps/libm-ieee754/s_llrint.c index faae106ece..a50ba96952 100644 --- a/sysdeps/libm-ieee754/s_llrint.c +++ b/sysdeps/libm-ieee754/s_llrint.c @@ -23,9 +23,7 @@ #include "math_private.h" -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ -static long double two52[2] = +static const long double two52[2] = { 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ @@ -33,210 +31,65 @@ static long double two52[2] = long long int -__llrint (long double x) +__llrint (double x) { - int32_t j0,sx; - u_int32_t i0, i1, i; - long double t, w; + int32_t j0; + u_int32_t i1, i0; long long int result; + volatile double w; + double t; + int sx; EXTRACT_WORDS (i0, i1, x); - - sx = i0 >> 31; j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + sx = i0 >> 31; + i0 &= 0xfffff; + i0 |= 0x100000; if (j0 < 20) { - if (j0 < 0) - { - if (((i0 & 0x7fffffff) | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xfffe0000; - i0 |= ((i1 | -i1) >> 12) & 0x80000; - SET_HIGH_WORD (x, i0); - w = two52[sx] + x; - t = w - two52[sx]; - GET_HIGH_WORD (i0, t); - if ((i0 & 0x7fffffff) >= 0x3fff0000) - result = sx ? -1 : 1; - else - result = 0; - } - } + if (j0 < -1) + return 0; else { - u_int32_t i = 0x000fffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 19) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000 >> j0); - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - } - } - - result = ((i0 >> (20 - j0)) & 0xfffff) | (0x00100000 >> (20 - j0)); - if (sx) - result = -result; + w = two52[sx] + x; + t = w - two52[sx]; + EXTRACT_WORDS (i0, i1, t); + i0 = i & 0xfffff; + i0 |= 0x100000; + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + result = i0 >> (20 - j0); } } - else if ((unsigned int) j0 < sizeof (long long int) * 8 && j0 < 53) - { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 20)); - } - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - - result = i0 | 0x00100000; - if (j0 > 20) - { - result <<= j0 - 20; - result |= i1 >> (52 - j0); - } - if (sx) - result = -result; - } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long long int) x; - - return result; -} - -#else -static long double two63[2] = -{ - 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ - -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ -}; - - -long long int -__llrint (long double x) -{ - int32_t se,j0,sx; - u_int32_t i0, i1, i; - long long int result; - long double w, t; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - - if (j0 < 31) + else if (j0 < (int32_t) (8 * sizeof (long long int))) { - if (j0 < 0) - { - if (((se & 0x7fff) | i0 | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xe0000000; - i0 |= (i1 | -i1) & 0x80000000; - SET_LDOUBLE_MSW (x, i0); - w = two63[sx] + x; - t = w - two63[sx]; - GET_LDOUBLE_EXP (i0, t); - if ((i0 & 0x7fff) >= 0x3fff) - result = sx ? -1 : 1; - else - result = 0; - } - } + if (j0 >= 52) + result = ((long long int) i0 << (j0 - 20)) | (i1 << (j0 - 52)); else { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 31) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000000 >> j0); - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - } - } - - - result = i0 >> (31 - j0); + w = two52[sx] + x; + t = w - two52[sx]; + EXTRACT_WORDS (i0, i1, t); + i0 = i & 0xfffff; + i0 |= 0x100000; + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + result = ((long long int) i0 << (j0 - 20)) | (j >> (52 - j0)); } } - else if ((unsigned int) j0 < sizeof (long long int) * 8 && j0 < 64) + else { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 31)); - } - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - j0 = (se & 0x7fff) - 0x3fff; - - result = i0; - if (j0 > 31) - { - result <<= j0 - 31; - result |= i1 >> (63 - j0); - } + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long long int) x; - return result; + return sx ? -result : result; } -#endif weak_alias (__llrint, llrint) +#ifdef NO_LONG_DOUBLE +strong_alias (__llrint, __llrintl) +weak_alias (__llrint, llrintl) +#endif diff --git a/sysdeps/libm-ieee754/s_llrintf.c b/sysdeps/libm-ieee754/s_llrintf.c new file mode 100644 index 0000000000..f73afd0fd2 --- /dev/null +++ b/sysdeps/libm-ieee754/s_llrintf.c @@ -0,0 +1,78 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +static const double two23[2] = +{ + 8.3886080000e+06, /* 0x4B000000 */ + -8.3886080000e+06, /* 0xCB000000 */ +}; + + +long long int +__llrintf (float x) +{ + int32_t j0; + u_int32_t i, i0; + volatile float w; + float t; + long long int result; + int sx; + + GET_FLOAT_WORD (i, x); + + sx = i0 >> 31; + j0 = ((i0 >> 23) & 0xff) - 0x7f; + i0 = i & 0x7fffff; + i0 |= 0x800000; + + if (j0 < (int32_t) (sizeof (long long int) * 8)) + { + if (j0 < -1) + return 0; + else if (j0 >= 23) + result = (long long int) i0 << (j0 - 23); + else + { + w = two23[sx] + x; + t = w - two23[sx]; + GET_FLOAT_WORD (i, t); + i0 = i & 0x7fffff; + i0 |= 0x800000; + j0 = ((i0 >> 23) & 0xff) - 0x7f; + + result = i0 >> (23 - j0); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; + } + + return sx ? -result : result; +} + +weak_alias (__llrintf, llrintf) diff --git a/sysdeps/libm-ieee754/s_llrintl.c b/sysdeps/libm-ieee754/s_llrintl.c new file mode 100644 index 0000000000..29fc3a300d --- /dev/null +++ b/sysdeps/libm-ieee754/s_llrintl.c @@ -0,0 +1,86 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +static const long double two63[2] = +{ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + + +long long int +__llrintl (long double x) +{ + int32_t se,j0; + u_int32_t i0, i1, i; + long long int result; + volatile long double w; + long double t; + int sx; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < 31) + { + if (j0 < -1) + return 0; + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + result = i0 >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long long int))) + { + if (j0 >= 63) + result = ((long long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + result = ((long long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; + } + + return sx ? -result : result; +} + +weak_alias (__llrintl, llrintl) diff --git a/sysdeps/libm-ieee754/s_llround.c b/sysdeps/libm-ieee754/s_llround.c index aee0e31fc0..1deb630e96 100644 --- a/sysdeps/libm-ieee754/s_llround.c +++ b/sysdeps/libm-ieee754/s_llround.c @@ -1,4 +1,4 @@ -/* Round long double value to long long int. +/* Round double value to long long int. Copyright (C) 1997 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -23,157 +23,56 @@ #include "math_private.h" -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ - long long int -__llround (long double x) +__llround (double x) { int32_t j0; u_int32_t i1, i0; long long int result; + int sign; EXTRACT_WORDS (i0, i1, x); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + sign = (i0 & 0x80000000) != 0 ? -1 : 1; + i0 &= 0xfffff; + i0 |= 0x100000; + if (j0 < 20) { if (j0 < 0) - result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); + return j0 < -1 ? 0 : sign; else { - u_int32_t i = 0xfffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long long int) ((i0 & 0xfffff) | 0x100000) >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x80000 >> j0); - if (j < i0) - result = (long long int) 0x80000 >> (20 - j0); - else - result = (j | 0x100000) >> (20 - j0); - } + i0 += 0x80000 >> j0; + + result = i0 >> (20 - j0); } } - else if (j0 >= 8 * sizeof (long long int) || j0 > 51) - { - /* The number is too large. It is left implementation defined - what happens. */ - result = (long long int) x; - } - else + else if (j0 < (int32_t) (8 * sizeof (long long int))) { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) + if (j0 >= 52) + result = ((long long int) i0 << (j0 - 20)) | (i1 << (j0 - 52)); + else { - /* x is not integral. */ u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); if (j < i1) - { - j = i0 + 1; - if ((j & 0xfffff) == 0) - { - if (sizeof (long long int) <= 4) - /* Overflow. */ - result = (long long int) x; - else - result = 1ll << (j0 + 1); - } - else - result = ((long long int) ((i0 & 0xfffff) | 0x100000) - << (j0 - 31)); - } - else - { - result = ((long long int) ((i0 & 0xfffff) | 0x100000) - << (j0 - 31)); - if (sizeof (long long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - if (sizeof (long long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } - } + ++i0; - return i0 & 0x80000000 ? -result : result; -} -#else -long long int -__llround (long double x) -{ - int32_t j0; - u_int32_t se, i1, i0; - long long int result; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - j0 = (se & 0x7fff) - 0x3fff; - if (j0 < 31) - { - if (j0 < 0) - result = j0 < -1 ? 0 : 1; - else - { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long long int) i0 >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x40000000 >> j0); - if (j < i0) - result = 0x80000000l >> (30 - j0); - else - result = j >> (31 - j0); - } + result = ((long long int) i0 << (j0 - 20)) | (j >> (52 - j0)); } } - else if ((unsigned int) j0 >= 8 * sizeof (long long int) || j0 > 62) + else { /* The number is too large. It is left implementation defined what happens. */ - result = (long long int) x; - } - else - { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); - if (j < i1) - { - j = i0 + 1; - if (j == 0) - { - if (sizeof (long long int) <= 4) - /* Overflow. */ - result = (long long int) x; - else - result = 1ll << (j0 + 1); - } - else - result = (long long int) i0 << (j0 - 31); - } - else - { - result = (long long int) i0 << (j0 - 31); - if (sizeof (long long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long long int) i0 << (j0 - 31); - if (sizeof (long long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } + return (long long int) x; } - return se & 0x8000 ? -result : result; + return sign * result; } -#endif + weak_alias (__llround, llround) +#ifdef NO_LONG_DOUBLE +strong_alias (__llround, __llroundl) +weak_alias (__llround, llroundl) +#endif diff --git a/sysdeps/libm-ieee754/s_llroundf.c b/sysdeps/libm-ieee754/s_llroundf.c new file mode 100644 index 0000000000..87ed1a3117 --- /dev/null +++ b/sysdeps/libm-ieee754/s_llroundf.c @@ -0,0 +1,63 @@ +/* Round float value to long long int. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long long int +__llroundf (float x) +{ + int32_t j0; + u_int32_t i; + long long int result; + int sign; + + GET_FLOAT_WORD (i, x); + j0 = ((i >> 23) & 0xff) - 0x7f; + sign = (i & 0x80000000) != 0 ? -1 : 1; + i &= 0x7fffff; + i |= 0x800000; + + if (j0 < (int32_t) (8 * sizeof (long long int))) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else if (j0 >= 23) + result = (long int) i << (j0 - 23); + else + { + i += 0x400000 >> j0; + + result = i >> (23 - j0); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; + } + + return sign * result; +} + +weak_alias (__llroundf, llroundf) diff --git a/sysdeps/libm-ieee754/s_llroundl.c b/sysdeps/libm-ieee754/s_llroundl.c new file mode 100644 index 0000000000..959b1e3fe9 --- /dev/null +++ b/sysdeps/libm-ieee754/s_llroundl.c @@ -0,0 +1,78 @@ +/* Round long double value to long long int. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long long int +__llroundl (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long long int result; + int sign; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + sign = (se & 0x8000) != 0 ? -1 : 1; + + if (j0 < 31) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else + { + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + { + j >>= 1; + j |= 0x80000000; + ++j0; + } + + result = j >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long long int))) + { + if (j0 >= 63) + result = ((long long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); + else + { + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + ++i0; + + result = ((long long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; + } + + return sign * result; +} + +weak_alias (__llroundl, llroundl) diff --git a/sysdeps/libm-ieee754/s_lrint.c b/sysdeps/libm-ieee754/s_lrint.c index 6779f974bd..0ccd06a061 100644 --- a/sysdeps/libm-ieee754/s_lrint.c +++ b/sysdeps/libm-ieee754/s_lrint.c @@ -23,9 +23,7 @@ #include "math_private.h" -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ -static long double two52[2] = +static const double two52[2] = { 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ @@ -33,210 +31,65 @@ static long double two52[2] = long int -__lrint (long double x) +__lrint (double x) { - int32_t j0,sx; + int32_t j0; u_int32_t i0,i1,i; - long double t, w; + volatile double w; + double t; long int result; + int sx; EXTRACT_WORDS (i0, i1, x); - - sx = i0 >> 31; j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + sx = i0 >> 31; + i0 &= 0xfffff; + i0 |= 0x100000; if (j0 < 20) { - if (j0 < 0) - { - if (((i0 & 0x7fffffff) | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xfffe0000; - i0 |= ((i1 | -i1) >> 12) & 0x80000; - SET_HIGH_WORD (x, i0); - w = two52[sx] + x; - t = w - two52[sx]; - GET_HIGH_WORD (i0, t); - if ((i0 & 0x7fffffff) >= 0x3fff0000) - result = sx ? -1 : 1; - else - result = 0; - } - } + if (j0 < -1) + return 0; else { - u_int32_t i = 0x000fffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 19) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000 >> j0); - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - } - } - - result = ((i0 >> (20 - j0)) & 0xfffff) | (0x00100000 >> (20 - j0)); - if (sx) - result = -result; + w = two52[sx] + x; + t = w - two52[sx]; + EXTRACT_WORDS (i0, i1, t); + i0 = i & 0xfffff; + i0 |= 0x100000; + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + result = i0 >> (20 - j0); } } - else if ((unsigned int) j0 < sizeof (long int) * 8 && j0 < 53) + else if (j0 < (int32_t) (8 * sizeof (long int))) { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 20)); - } - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - - result = i0 | 0x00100000; - if (j0 > 20) - { - result <<= j0 - 20; - result |= i1 >> (52 - j0); - } - if (sx) - result = -result; - } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long int) x; - - return result; -} - -#else -static long double two63[2] = -{ - 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ - -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ -}; - - -long int -__lrint (long double x) -{ - int32_t se,j0,sx; - u_int32_t i0,i1,i; - long int result; - long double w, t; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - - if (j0 < 31) - { - if (j0 < 0) - { - if (((se & 0x7fff) | i0 | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xe0000000; - i0 |= (i1 | -i1) & 0x80000000; - SET_LDOUBLE_MSW (x, i0); - w = two63[sx] + x; - t = w - two63[sx]; - GET_LDOUBLE_EXP (i0, t); - if ((i0 & 0x7fff) >= 0x3fff) - result = sx ? -1 : 1; - else - result = 0; - } - } + if (j0 >= 52) + result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52)); else { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 31) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000000 >> j0); - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - } - } - - - result = i0 >> (31 - j0); + w = two52[sx] + x; + t = w - two52[sx]; + EXTRACT_WORDS (i0, i1, t); + i0 = i & 0xfffff; + i0 |= 0x100000; + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + result = ((long int) i0 << (j0 - 20)) | (j >> (52 - j0)); } } - else if ((unsigned int) j0 < sizeof (long int) * 8 && j0 < 64) + else { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 31)); - } - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - j0 = (se & 0x7fff) - 0x3fff; - - result = i0; - if (j0 > 31) - { - result <<= j0 - 31; - result |= i1 >> (63 - j0); - } + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long int) x; - return result; + return sx ? -result : result; } -#endif weak_alias (__lrint, lrint) +#ifdef NO_LONG_DOUBLE +strong_alias (__lrint, __lrintl) +weak_alias (__lrint, lrintl) +#endif diff --git a/sysdeps/libm-ieee754/s_lrintf.c b/sysdeps/libm-ieee754/s_lrintf.c new file mode 100644 index 0000000000..6d74b362f5 --- /dev/null +++ b/sysdeps/libm-ieee754/s_lrintf.c @@ -0,0 +1,78 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +static const double two23[2] = +{ + 8.3886080000e+06, /* 0x4B000000 */ + -8.3886080000e+06, /* 0xCB000000 */ +}; + + +long int +__lrintf (float x) +{ + int32_t j0; + u_int32_t i, i0; + volatile float w; + float t; + long int result; + int sx; + + GET_FLOAT_WORD (i, x); + + sx = i0 >> 31; + j0 = ((i0 >> 23) & 0xff) - 0x7f; + i0 = i & 0x7fffff; + i0 |= 0x800000; + + if (j0 < (int32_t) (sizeof (long int) * 8)) + { + if (j0 < -1) + return 0; + else if (j0 >= 23) + result = (long int) i0 << (j0 - 23); + else + { + w = two23[sx] + x; + t = w - two23[sx]; + GET_FLOAT_WORD (i, t); + i0 = i & 0x7fffff; + i0 |= 0x800000; + j0 = ((i0 >> 23) & 0xff) - 0x7f; + + result = i0 >> (23 - j0); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sx ? -result : result; +} + +weak_alias (__lrintf, lrintf) diff --git a/sysdeps/libm-ieee754/s_lrintl.c b/sysdeps/libm-ieee754/s_lrintl.c new file mode 100644 index 0000000000..8607851669 --- /dev/null +++ b/sysdeps/libm-ieee754/s_lrintl.c @@ -0,0 +1,86 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +static const long double two63[2] = +{ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + + +long int +__lrintl (long double x) +{ + int32_t se,j0; + u_int32_t i0,i1,i; + long int result; + volatile long double w; + long double t; + int sx; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < 31) + { + if (j0 < -1) + return 0; + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + result = i0 >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long int))) + { + if (j0 >= 63) + result = ((long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + result = ((long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sx ? -result : result; +} + +weak_alias (__lrintl, lrintl) diff --git a/sysdeps/libm-ieee754/s_lround.c b/sysdeps/libm-ieee754/s_lround.c index 974dbde050..a6468ba78a 100644 --- a/sysdeps/libm-ieee754/s_lround.c +++ b/sysdeps/libm-ieee754/s_lround.c @@ -1,4 +1,4 @@ -/* Round long double value to long int. +/* Round double value to long int. Copyright (C) 1997 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -23,13 +23,8 @@ #include "math_private.h" -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ - -/* This code does not presently work. */ - long int -__lround (long double x) +__lround (double x) { int32_t j0; u_int32_t i1, i0; @@ -40,138 +35,44 @@ __lround (long double x) j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; sign = (i0 & 0x80000000) != 0 ? -1 : 1; i0 &= 0xfffff; + i0 |= 0x100000; + if (j0 < 20) { if (j0 < 0) - result = j0 < -1 ? 0 : sign; + return j0 < -1 ? 0 : sign; else { - u_int32_t i = 0xfffff >> j0; - i0 |= 0x100000; - if (((i0 & i) | i1) == 0) - result = i0 >> j0; - else - /* X is not integral. */ - result = (i0 + (0x80000 >> j0)) >> (20 - j0); + i0 += 0x80000 >> j0; + + result = i0 >> (20 - j0); } } - else if (j0 >= 8 * sizeof (long int) || j0 > 51) - { - /* The number is too large. It is left implementation defined - what happens. */ - return (long int) x; - } - else + else if (j0 < (int32_t) (8 * sizeof (long int))) { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) + if (j0 >= 52) + result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52)); + else { - /* x is not integral. */ u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); if (j < i1) - { - j = i0 + 1; - if ((j & 0xfffff) == 0) - { - if (sizeof (long int) <= 4) - /* Overflow. */ - result = (long int) x; - else - result = 1l << (j0 + 1); - } - else - result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - } - else - { - result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } - } - - return sign * result; -} -#else -long int -__lround (long double x) -{ - int32_t j0; - u_int32_t se, i1, i0; - long int result; + ++i0; - GET_LDOUBLE_WORDS (se, i0, i1, x); - j0 = (se & 0x7fff) - 0x3fff; - if (j0 < 31) - { - if (j0 < 0) - result = j0 < -1 ? 0 : 1; - else - { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long int) i0 >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x40000000 >> j0); - if (j < i0) - result = 0x80000000l >> (30 - j0); - else - result = j >> (31 - j0); - } + result = ((long int) i0 << (j0 - 20)) | (j >> (52 - j0)); } } - else if ((unsigned int) j0 >= 8 * sizeof (long int) || j0 > 62) + else { /* The number is too large. It is left implementation defined what happens. */ - result = (long int) x; - } - else - { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); - if (j < i1) - { - j = i0 + 1; - if (j == 0) - { - if (sizeof (long int) <= 4) - /* Overflow. */ - result = (long int) x; - else - result = 1l << (j0 + 1); - } - else - result = (long int) i0 << (j0 - 31); - } - else - { - result = (long int) i0 << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long int) i0 << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } + return (long int) x; } - return se & 0x8000 ? -result : result; + return sign * result; } -#endif + weak_alias (__lround, lround) +#ifdef NO_LONG_DOUBLE +strong_alias (__lround, __lroundl) +weak_alias (__lround, lroundl) +#endif diff --git a/sysdeps/libm-ieee754/s_lroundf.c b/sysdeps/libm-ieee754/s_lroundf.c new file mode 100644 index 0000000000..0d3889c24a --- /dev/null +++ b/sysdeps/libm-ieee754/s_lroundf.c @@ -0,0 +1,63 @@ +/* Round float value to long int. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long int +__lroundf (float x) +{ + int32_t j0; + u_int32_t i; + long int result; + int sign; + + GET_FLOAT_WORD (i, x); + j0 = ((i >> 23) & 0xff) - 0x7f; + sign = (i & 0x80000000) != 0 ? -1 : 1; + i &= 0x7fffff; + i |= 0x800000; + + if (j0 < (int32_t) (8 * sizeof (long int))) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else if (j0 >= 23) + result = (long int) i << (j0 - 23); + else + { + i += 0x400000 >> j0; + + result = i >> (23 - j0); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sign * result; +} + +weak_alias (__lroundf, lroundf) diff --git a/sysdeps/libm-ieee754/s_lroundl.c b/sysdeps/libm-ieee754/s_lroundl.c new file mode 100644 index 0000000000..d8b7db4b5c --- /dev/null +++ b/sysdeps/libm-ieee754/s_lroundl.c @@ -0,0 +1,78 @@ +/* Round long double value to long int. + Copyright (C) 1997 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 Library General Public License as + published by the Free Software Foundation; either version 2 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long int +__lroundl (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long int result; + int sign; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + sign = (se & 0x8000) != 0 ? -1 : 1; + + if (j0 < 31) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else + { + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + { + j >>= 1; + j |= 0x80000000; + ++j0; + } + + result = j >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long int))) + { + if (j0 >= 63) + result = ((long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); + else + { + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + ++i0; + + result = ((long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sign * result; +} + +weak_alias (__lroundl, lroundl) diff --git a/sysdeps/libm-ieee754/s_nextafterl.c b/sysdeps/libm-ieee754/s_nextafterl.c index 1d3fc86c44..d7bdf026e7 100644 --- a/sysdeps/libm-ieee754/s_nextafterl.c +++ b/sysdeps/libm-ieee754/s_nextafterl.c @@ -97,3 +97,5 @@ static char rcsid[] = "$NetBSD: $"; return x; } weak_alias (__nextafterl, nextafterl) +strong_alias (__nextafterl, __nextafterxl) +weak_alias (__nextafterl, nextafterxl) diff --git a/sysdeps/libm-ieee754/s_nextafterx.c b/sysdeps/libm-ieee754/s_nextafterx.c new file mode 100644 index 0000000000..745ceaf34d --- /dev/null +++ b/sysdeps/libm-ieee754/s_nextafterx.c @@ -0,0 +1,98 @@ +/* s_nextafterx.c + * Conversion freom s_nextafter.c by Ulrich Drepper, Cygnus Support, + * drepper@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* IEEE functions + * nextafterx(x,y) + * return the next machine floating-point number of x in the + * direction toward y. + * Special cases: + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + double __nextafterx(double x, long double y) +#else + double __nextafterx(x,y) + double x; + long double y; +#endif +{ + int32_t hx,ix,iy; + u_int32_t lx,hy,ly,esy; + + EXTRACT_WORDS(hx,lx,x); + GET_LDOUBLE_WORDS(esy,hy,ly,y); + ix = hx&0x7fffffff; /* |x| */ + iy = esy&0x7fff; /* |y| */ + + if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ + ((iy>=0x7fff)&&((hy|ly)|-(hy|ly))!=0)) /* y is nan */ + return x+y; + if((long double) x==y) return x; /* x=y, return x */ + if((ix|lx)==0) { /* x == 0 */ + double x2; + INSERT_WORDS(x,esy&0x8000?0x80000000:0,1);/* return +-minsub */ + x2 = x*x; + if(x2==x) return x2; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if (esy>=0x8000||((ix>>20)&0x7ff)>iy + || (((ix>>20)&0x7ff)==iy + && (((hx<<11)|(lx>>21))>(hy&0x7fffffff) + || (((hx<<11)|(lx>>21))==(hy&0x7fffffff) + && (lx<<11)>ly)))) { /* x > y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x < y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } else { /* x < 0 */ + if (esy<0x8000||((ix>>20)&0x7ff)>iy + || (((ix>>20)&0x7ff)==iy + && (((hx<<11)|(lx>>21))>(hy&0x7fffffff) + || (((hx<<11)|(lx>>21))==(hy&0x7fffffff) + && (lx<<11)>ly)))) {/* x < y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x > y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } + hy = hx&0x7ff00000; + if(hy>=0x7ff00000) return x+x; /* overflow */ + if(hy<0x00100000) { /* underflow */ + double x2 = x*x; + if(x2!=x) { /* raise underflow flag */ + INSERT_WORDS(x2,hx,lx); + return x2; + } + } + INSERT_WORDS(x,hx,lx); + return x; +} +weak_alias (__nextafterx, nextafterx) +#ifdef NO_LONG_DOUBLE +strong_alias (__nextafterx, __nextafterxl) +weak_alias (__nextafterx, nextafterxl) +#endif diff --git a/sysdeps/libm-ieee754/s_nextafterxf.c b/sysdeps/libm-ieee754/s_nextafterxf.c new file mode 100644 index 0000000000..7e870805e2 --- /dev/null +++ b/sysdeps/libm-ieee754/s_nextafterxf.c @@ -0,0 +1,78 @@ +/* s_nextafterxf.c -- float version of s_nextafter.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + float __nextafterxf(float x, long double y) +#else + float __nextafterxf(x,y) + float x; + log double y; +#endif +{ + int32_t hx,ix,iy; + u_int32_t hy,ly,esy; + + GET_FLOAT_WORD(hx,x); + GET_LDOUBLE_WORDS(esy,hy,ly,y); + ix = hx&0x7fffffff; /* |x| */ + iy = esy&0x7fff; /* |y| */ + + if((ix>0x7f800000) || /* x is nan */ + (iy>=0x7fff&&((hy|ly)|-(hy|ly))!=0)) /* y is nan */ + return x+y; + if((long double) x==y) return y; /* x=y, return y */ + if(ix==0) { /* x == 0 */ + float x2; + SET_FLOAT_WORD(x,(esy&0x8000?0x80000000:0)|1);/* return +-minsub*/ + x2 = x*x; + if(x2==x) return x2; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if(esy>=0x8000||((ix>>23)&0xff)>iy + || (((ix>>23)&0xff)==iy + && ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x > y, x -= ulp */ + hx -= 1; + } else { /* x < y, x += ulp */ + hx += 1; + } + } else { /* x < 0 */ + if(esy<0x8000||((ix>>23)&0xff)>iy + || (((ix>>23)&0xff)==iy + && ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x < y, x -= ulp */ + hx -= 1; + } else { /* x > y, x += ulp */ + hx += 1; + } + } + hy = hx&0x7f800000; + if(hy>=0x7f800000) return x+x; /* overflow */ + if(hy<0x00800000) { /* underflow */ + float x2 = x*x; + if(x2!=x) { /* raise underflow flag */ + SET_FLOAT_WORD(x2,hx); + return x2; + } + } + SET_FLOAT_WORD(x,hx); + return x; +} +weak_alias (__nextafterxf, nextafterxf) diff --git a/sysdeps/libm-ieee754/s_nextafterxl.c b/sysdeps/libm-ieee754/s_nextafterxl.c new file mode 100644 index 0000000000..73c3610fc1 --- /dev/null +++ b/sysdeps/libm-ieee754/s_nextafterxl.c @@ -0,0 +1 @@ +/* This function is the same as nextafterl so we use an alias there. */ diff --git a/sysdeps/libm-ieee754/s_scalbln.c b/sysdeps/libm-ieee754/s_scalbln.c new file mode 100644 index 0000000000..aa6134f093 --- /dev/null +++ b/sysdeps/libm-ieee754/s_scalbln.c @@ -0,0 +1,70 @@ +/* @(#)s_scalbn.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $"; +#endif + +/* + * scalbn (double x, int n) + * scalbn(x,n) returns x* 2**n computed by exponent + * manipulation rather than by actually performing an + * exponentiation or a multiplication. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ +twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ +huge = 1.0e+300, +tiny = 1.0e-300; + +#ifdef __STDC__ + double __scalbln (double x, long int n) +#else + double __scalbln (x,n) + double x; long int n; +#endif +{ + int32_t k,hx,lx; + EXTRACT_WORDS(hx,lx,x); + k = (hx&0x7ff00000)>>20; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ + x *= two54; + GET_HIGH_WORD(hx,x); + k = ((hx&0x7ff00000)>>20) - 54; + } + if (k==0x7ff) return x+x; /* NaN or Inf */ + k = k+n; + if (n> 50000 || k > 0x7fe) + return huge*__copysign(huge,x); /* overflow */ + if (n< -50000) return tiny*__copysign(tiny,x); /*underflow*/ + if (k > 0) /* normal result */ + {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} + if (k <= -54) + return tiny*__copysign(tiny,x); /*underflow*/ + k += 54; /* subnormal result */ + SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); + return x*twom54; +} +weak_alias (__scalbln, scalbln) +#ifdef NO_LONG_DOUBLE +strong_alias (__scalbln, __scalblnl) +weak_alias (__scalbln, scalblnl) +#endif diff --git a/sysdeps/libm-ieee754/s_scalblnf.c b/sysdeps/libm-ieee754/s_scalblnf.c new file mode 100644 index 0000000000..4ed48733cf --- /dev/null +++ b/sysdeps/libm-ieee754/s_scalblnf.c @@ -0,0 +1,63 @@ +/* s_scalbnf.c -- float version of s_scalbn.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_scalbnf.c,v 1.4 1995/05/10 20:48:10 jtc Exp $"; +#endif + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +two25 = 3.355443200e+07, /* 0x4c000000 */ +twom25 = 2.9802322388e-08, /* 0x33000000 */ +huge = 1.0e+30, +tiny = 1.0e-30; + +#ifdef __STDC__ + float __scalblnf (float x, long int n) +#else + float __scalblnf (x,n) + float x; long int n; +#endif +{ + int32_t k,ix; + GET_FLOAT_WORD(ix,x); + k = (ix&0x7f800000)>>23; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((ix&0x7fffffff)==0) return x; /* +-0 */ + x *= two25; + GET_FLOAT_WORD(ix,x); + k = ((ix&0x7f800000)>>23) - 25; + } + if (k==0xff) return x+x; /* NaN or Inf */ + k = k+n; + if (n> 50000 || k > 0xfe) + return huge*copysignf(huge,x); /* overflow */ + if (n< -50000) + return tiny*copysignf(tiny,x); /*underflow*/ + if (k > 0) /* normal result */ + {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;} + if (k <= -25) + return tiny*copysignf(tiny,x); /*underflow*/ + k += 25; /* subnormal result */ + SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); + return x*twom25; +} +weak_alias (__scalblnf, scalblnf) diff --git a/sysdeps/libm-ieee754/s_scalblnl.c b/sysdeps/libm-ieee754/s_scalblnl.c new file mode 100644 index 0000000000..aafb73e688 --- /dev/null +++ b/sysdeps/libm-ieee754/s_scalblnl.c @@ -0,0 +1,71 @@ +/* s_scalbnl.c -- long double version of s_scalbn.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* + * scalbnl (long double x, int n) + * scalbnl(x,n) returns x* 2**n computed by exponent + * manipulation rather than by actually performing an + * exponentiation or a multiplication. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +two63 = 4.50359962737049600000e+15, +twom63 = 1.08420217248550443400e-19; +huge = 1.0e+4900L, +tiny = 1.0e-4900L; + +#ifdef __STDC__ + long double __scalblnl (long double x, long int n) +#else + long double __scalblnl (x,n) + long double x; long int n; +#endif +{ + int32_t k,es,hx,lx; + GET_LDOUBLE_WORDS(es,hx,lx,x); + k = es&0x7fff; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ + x *= two63; + GET_LDOUBLE_EXP(es,x); + k = (hx&0x7fff) - 63; + } + if (k==0x7fff) return x+x; /* NaN or Inf */ + k = k+n; + if (n> 50000 || k > 0x7ffe) + return huge*__copysignl(huge,x); /* overflow */ + if (n< -50000) + return tiny*__copysignl(tiny,x); + if (k > 0) /* normal result */ + {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;} + if (k <= -63) + return tiny*__copysignl(tiny,x); /*underflow*/ + k += 54; /* subnormal result */ + SET_LDOUBLE_EXP(x,(es&0x8000)|k); + return x*twom63; +} +weak_alias (__scalblnl, scalblnl) diff --git a/sysdeps/libm-ieee754/s_scalbn.c b/sysdeps/libm-ieee754/s_scalbn.c index 0ce0ffdf2f..3dbfe8fef0 100644 --- a/sysdeps/libm-ieee754/s_scalbn.c +++ b/sysdeps/libm-ieee754/s_scalbn.c @@ -35,10 +35,10 @@ huge = 1.0e+300, tiny = 1.0e-300; #ifdef __STDC__ - double __scalbn (double x, long int n) + double __scalbn (double x, int n) #else double __scalbn (x,n) - double x; long int n; + double x; int n; #endif { int32_t k,hx,lx; diff --git a/sysdeps/libm-ieee754/s_scalbnf.c b/sysdeps/libm-ieee754/s_scalbnf.c index 4799c825c9..11b77bd6c2 100644 --- a/sysdeps/libm-ieee754/s_scalbnf.c +++ b/sysdeps/libm-ieee754/s_scalbnf.c @@ -31,10 +31,10 @@ huge = 1.0e+30, tiny = 1.0e-30; #ifdef __STDC__ - float __scalbnf (float x, long int n) + float __scalbnf (float x, int n) #else float __scalbnf (x,n) - float x; long int n; + float x; int n; #endif { int32_t k,ix; diff --git a/sysdeps/libm-ieee754/s_scalbnl.c b/sysdeps/libm-ieee754/s_scalbnl.c index 3e80d85136..4e1d08ef90 100644 --- a/sysdeps/libm-ieee754/s_scalbnl.c +++ b/sysdeps/libm-ieee754/s_scalbnl.c @@ -39,10 +39,10 @@ huge = 1.0e+4900L, tiny = 1.0e-4900L; #ifdef __STDC__ - long double __scalbnl (long double x, long int n) + long double __scalbnl (long double x, int n) #else long double __scalbnl (x,n) - long double x; long int n; + long double x; int n; #endif { int32_t k,es,hx,lx; |