From 84057895201ee244894e8847c933df6db5aa25a3 Mon Sep 17 00:00:00 2001 From: Ulrich Drepper Date: Tue, 25 Aug 2009 17:28:50 -0700 Subject: [PATCH] Add 64-bit optimized s_round. --- ChangeLog | 4 ++ sysdeps/ieee754/dbl-64/wordsize-64/s_round.c | 77 ++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 sysdeps/ieee754/dbl-64/wordsize-64/s_round.c diff --git a/ChangeLog b/ChangeLog index 4b1faea..910e188 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,7 @@ +2009-08-25 Ulrich Drepper + + * sysdeps/ieee754/dbl-64/wordsize-64/s_round.c: New file. + 2009-08-25 Joseph Myers * math/s_fdiml.c (__fdiml): Use fpclassify instead of fpclassifyl. 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 0000000..5bd8579 --- /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 , 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 + +#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 -- 2.7.4