diff options
author | Roland McGrath <roland@redhat.com> | 2009-08-27 13:41:23 -0700 |
---|---|---|
committer | Roland McGrath <roland@redhat.com> | 2009-08-27 13:41:23 -0700 |
commit | 8a344d4a3bd01d1704a405b741c99086261713cf (patch) | |
tree | 1d28740af54a3d9194373f7f561baf579c3a8f44 /sysdeps/ieee754/dbl-64 | |
parent | bc9dd948ffc9cf1992d3e77e753789ab897e3dc6 (diff) | |
parent | 2df4be8c3545f6158e468660a73f794573fb6f07 (diff) |
Merge commit 'origin/master' into fedora/master
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c | 67 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c | 30 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c | 39 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c | 1 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c | 80 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c | 86 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c | 79 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_round.c | 77 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c | 68 | ||||
-rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c | 56 |
10 files changed, 583 insertions, 0 deletions
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c new file mode 100644 index 0000000000..9123fdc7bd --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c @@ -0,0 +1,67 @@ +/* @(#)s_ceil.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. + * ==================================================== + */ + +/* + * ceil(x) + * Return x rounded toward -inf to integral value + * Method: + * Bit twiddling. + * Exception: + * Inexact flag raised if x not equal to ceil(x). + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const double huge = 1.0e300; +#else +static double huge = 1.0e300; +#endif + +#ifdef __STDC__ + double __ceil(double x) +#else + double __ceil(x) + double x; +#endif +{ + int64_t i0,i; + int32_t j0; + EXTRACT_WORDS64(i0,x); + j0 = ((i0>>52)&0x7ff)-0x3ff; + if(j0<=51) { + if(j0<0) { /* raise inexact if x != 0 */ + if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=INT64_C(0x8000000000000000);} + else if(i0!=0) { i0=INT64_C(0x3ff0000000000000);} + } + } else { + i = INT64_C(0x000fffffffffffff)>>j0; + if((i0&i)==0) return x; /* x is integral */ + if(huge+x>0.0) { /* raise inexact flag */ + if(i0>0) i0 += UINT64_C(0x0010000000000000)>>j0; + i0 &= (~i); + } + } + } else { + if(j0==0x400) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } + INSERT_WORDS64(x,i0); + return x; +} +weak_alias (__ceil, ceil) +#ifdef NO_LONG_DOUBLE +strong_alias (__ceil, __ceill) +weak_alias (__ceil, ceill) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c new file mode 100644 index 0000000000..c48e979a84 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c @@ -0,0 +1,30 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changed to return -1 for -Inf by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + +/* + * isinf(x) returns 1 is x is inf, -1 if x is -inf, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +int +__isinf (double x) +{ + int64_t ix; + EXTRACT_WORDS64(ix,x); + int64_t t = ix & UINT64_C(0x7fffffffffffffff); + t ^= UINT64_C(0x7ff0000000000000); + t |= -t; + return ~(t >> 63) & (ix >> 62); +} +hidden_def (__isinf) +weak_alias (__isinf, isinf) +#ifdef NO_LONG_DOUBLE +strong_alias (__isinf, __isinfl) +weak_alias (__isinf, isinfl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c new file mode 100644 index 0000000000..3b08c54dd4 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c @@ -0,0 +1,39 @@ +/* @(#)s_isnan.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. + * ==================================================== + */ + +/* + * isnan(x) returns 1 is x is nan, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + int __isnan(double x) +#else + int __isnan(x) + double x; +#endif +{ + int64_t hx; + EXTRACT_WORDS64(hx,x); + hx &= UINT64_C(0x7fffffffffffffff); + hx = UINT64_C(0x7ff0000000000000) - hx; + return (int)(((uint64_t)hx)>>63); +} +hidden_def (__isnan) +weak_alias (__isnan, isnan) +#ifdef NO_LONG_DOUBLE +strong_alias (__isnan, __isnanl) +weak_alias (__isnan, isnanl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c new file mode 100644 index 0000000000..f7c9ea568f --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c @@ -0,0 +1 @@ +/* The code is the same as lround. Use an alias, see l_round.c. */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c new file mode 100644 index 0000000000..30ea5db7cc --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c @@ -0,0 +1,80 @@ +/* Round double value to long int. + Copyright (C) 1997, 2004, 2009 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. */ + +#define llround __hidden_llround +#define __llround __hidden___llround + +#include <math.h> + +#include "math_private.h" + + +long int +__lround (double x) +{ + int32_t j0; + int64_t i0; + long int result; + int sign; + + EXTRACT_WORDS64 (i0, x); + j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; + sign = i0 < 0 ? -1 : 1; + i0 &= UINT64_C(0xfffffffffffff); + i0 |= UINT64_C(0x10000000000000); + + if (j0 < (int32_t) (8 * sizeof (long int)) - 1) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else if (j0 >= 52) + result = i0 << (j0 - 52); + else + { + i0 += UINT64_C(0x8000000000000) >> j0; + + result = i0 >> (52 - j0); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sign * result; +} + +weak_alias (__lround, lround) +#ifdef NO_LONG_DOUBLE +strong_alias (__lround, __lroundl) +weak_alias (__lround, lroundl) +#endif + +/* long long has the same width as long on 64-bit machines. */ +#undef llround +#undef __llround +strong_alias (__lround, __llround) +weak_alias (__lround, llround) +#ifdef NO_LONG_DOUBLE +strong_alias (__lround, __llroundl) +weak_alias (__lround, llroundl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c new file mode 100644 index 0000000000..cb49019ddb --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c @@ -0,0 +1,86 @@ +/* Adapted for use as nearbyint by Ulrich Drepper <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. + * ==================================================== + */ + +/* + * rint(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rint(x). + */ + +#include <fenv.h> +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +TWO52[2]={ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + +#ifdef __STDC__ + double __nearbyint(double x) +#else + double __nearbyint(x) + double x; +#endif +{ + fenv_t env; + int64_t i0,sx; + int32_t j0; + EXTRACT_WORDS64(i0,x); + sx = (i0>>63)&1; + j0 = ((i0>>52)&0x7ff)-0x3ff; + if(j0<52) { + if(j0<0) { + if((i0&UINT64_C(0x7fffffffffffffff))==0) return x; + uint64_t i = i0 & UINT64_C(0xfffffffffffff); + i0 &= UINT64_C(0xfffe000000000000); + i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000)); + INSERT_WORDS64(x,i0); + feholdexcept (&env); + double w = TWO52[sx]+x; + double t = w-TWO52[sx]; + fesetenv (&env); + EXTRACT_WORDS64(i0,t); + INSERT_WORDS64(t,(i0&UINT64_C(0x7fffffffffffffff))|(sx<<63)); + return t; + } else { + uint64_t i = UINT64_C(0x000fffffffffffff)>>j0; + if((i0&i)==0) return x; /* x is integral */ + i>>=1; + if((i0&i)!=0) + i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0); + } + } else { + if(j0==0x400) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } + INSERT_WORDS64(x,i0); + feholdexcept (&env); + double w = TWO52[sx]+x; + double t = w-TWO52[sx]; + fesetenv (&env); + return t; +} +weak_alias (__nearbyint, nearbyint) +#ifdef NO_LONG_DOUBLE +strong_alias (__nearbyint, __nearbyintl) +weak_alias (__nearbyint, nearbyintl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c new file mode 100644 index 0000000000..4a60aa3278 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c @@ -0,0 +1,79 @@ +/* @(#)s_rint.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. + * ==================================================== + */ + +/* + * rint(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rint(x). + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +TWO52[2]={ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + +#ifdef __STDC__ + double __rint(double x) +#else + double __rint(x) + double x; +#endif +{ + int64_t i0,sx; + int32_t j0; + EXTRACT_WORDS64(i0,x); + sx = (i0>>63)&1; + j0 = ((i0>>52)&0x7ff)-0x3ff; + if(j0<52) { + if(j0<0) { + if((i0 & UINT64_C(0x7fffffffffffffff))==0) return x; + uint64_t i = i0 & UINT64_C(0xfffffffffffff); + i0 &= UINT64_C(0xfffe000000000000); + i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000)); + INSERT_WORDS64(x,i0); + double w = TWO52[sx]+x; + double t = w-TWO52[sx]; + EXTRACT_WORDS64(i0,t); + INSERT_WORDS64(t,(i0&UINT64_C(0x7fffffffffffffff))|(sx<<63)); + return t; + } else { + uint64_t i = UINT64_C(0x000fffffffffffff)>>j0; + if((i0&i)==0) return x; /* x is integral */ + i>>=1; + if((i0&i)!=0) + i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0); + } + } else { + if(j0==0x400) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } + INSERT_WORDS64(x,i0); + double w = TWO52[sx]+x; + return w-TWO52[sx]; +} +weak_alias (__rint, rint) +#ifdef NO_LONG_DOUBLE +strong_alias (__rint, __rintl) +weak_alias (__rint, rintl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c new file mode 100644 index 0000000000..5bd857910c --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -0,0 +1,77 @@ +/* Round double to integer away from zero. + Copyright (C) 1997, 2009 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 <math.h> + +#include "math_private.h" + + +static const double huge = 1.0e300; + + +double +__round (double x) +{ + int64_t i0, j0; + + EXTRACT_WORDS64 (i0, x); + j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; + if (__builtin_expect (j0 < 52, 1)) + { + if (j0 < 0) + { + if (huge + x > 0.0) + { + i0 &= UINT64_C(0x8000000000000000); + if (j0 == -1) + i0 |= UINT64_C(0x3ff0000000000000); + } + } + else + { + uint64_t i = UINT64_C(0x000fffffffffffff) >> j0; + if ((i0 & i) == 0) + /* X is integral. */ + return x; + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + i0 += UINT64_C(0x0008000000000000) >> j0; + i0 &= ~i; + } + } + } + else + { + if (j0 == 0x400) + /* Inf or NaN. */ + return x + x; + else + return x; + } + + INSERT_WORDS64 (x, i0); + return x; +} +weak_alias (__round, round) +#ifdef NO_LONG_DOUBLE +strong_alias (__round, __roundl) +weak_alias (__round, roundl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c new file mode 100644 index 0000000000..25cf3b1699 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_scalbn.c @@ -0,0 +1,68 @@ +/* @(#)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. + * ==================================================== + */ + +/* + * 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 __scalbn (double x, int n) +#else + double __scalbn (x,n) + double x; int n; +#endif +{ + int64_t ix; + int64_t k; + EXTRACT_WORDS64(ix,x); + k = (ix >> 52) & 0x7ff; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((ix & UINT64_C(0xfffffffffffff))==0) return x; /* +-0 */ + x *= two54; + EXTRACT_WORDS64(ix,x); + k = ((ix >> 52) & 0x7ff) - 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 */ + {INSERT_WORDS64(x,(ix&UINT64_C(0x800fffffffffffff))|(k<<52)); + return x;} + if (k <= -54) + return tiny*__copysign(tiny,x); /*underflow*/ + k += 54; /* subnormal result */ + INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52)); + return x*twom54; +} +weak_alias (__scalbn, scalbn) +#ifdef NO_LONG_DOUBLE +strong_alias (__scalbn, __scalbnl) +weak_alias (__scalbn, scalbnl) +#endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c new file mode 100644 index 0000000000..9add5ade71 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c @@ -0,0 +1,56 @@ +/* Truncate argument to nearest integral value not larger than the argument. + Copyright (C) 1997, 1998, 2009 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 <math.h> + +#include "math_private.h" + + +double +__trunc (double x) +{ + int64_t i0, j0; + int64_t sx; + + EXTRACT_WORDS64 (i0, x); + sx = i0 & UINT64_C(0x8000000000000000); + j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; + if (j0 < 52) + { + if (j0 < 0) + /* The magnitude of the number is < 1 so the result is +-0. */ + INSERT_WORDS64 (x, sx); + else + INSERT_WORDS64 (x, sx | (i0 & ~(UINT64_C(0x000fffffffffffff) >> j0))); + } + else + { + if (j0 == 0x400) + /* x is inf or NaN. */ + return x + x; + } + + return x; +} +weak_alias (__trunc, trunc) +#ifdef NO_LONG_DOUBLE +strong_alias (__trunc, __truncl) +weak_alias (__trunc, truncl) +#endif |