diff options
Diffstat (limited to 'src/math/rintl.c')
| -rw-r--r-- | src/math/rintl.c | 86 |
1 files changed, 16 insertions, 70 deletions
diff --git a/src/math/rintl.c b/src/math/rintl.c index b13cfeb3..26725073 100644 --- a/src/math/rintl.c +++ b/src/math/rintl.c @@ -1,30 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_rintl.c */ -/*- - * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG> - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -33,53 +6,26 @@ long double rintl(long double x) return rint(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -#define BIAS (LDBL_MAX_EXP - 1) - -static const float -shift[2] = { #if LDBL_MANT_DIG == 64 - 0x1.0p63, -0x1.0p63 +#define TOINT 0x1p63 #elif LDBL_MANT_DIG == 113 - 0x1.0p112, -0x1.0p112 +#define TOINT 0x1p112 #endif -}; -static const float zero[2] = { 0.0, -0.0 }; - long double rintl(long double x) { - union IEEEl2bits u; - uint32_t expsign; - int ex, sign; - - u.e = x; - expsign = u.xbits.expsign; - ex = expsign & 0x7fff; - - if (ex >= BIAS + LDBL_MANT_DIG - 1) { - if (ex == BIAS + LDBL_MAX_EXP) - return x + x; /* Inf, NaN, or unsupported format */ - return x; /* finite and already an integer */ - } - sign = expsign >> 15; - - /* - * The following code assumes that intermediate results are - * evaluated in long double precision. If they are evaluated in - * greater precision, double rounding may occur, and if they are - * evaluated in less precision (as on i386), results will be - * wildly incorrect. - */ - x += shift[sign]; - x -= shift[sign]; - - /* - * If the result is +-0, then it must have the same sign as x, but - * the above calculation doesn't always give this. Fix up the sign. - */ - if (ex < BIAS && x == 0.0) - return zero[sign]; - - return x; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + int s = u.i.se >> 15; + long double y; + + if (e >= 0x3fff+LDBL_MANT_DIG-1) + return x; + if (s) + y = x - TOINT + TOINT; + else + y = x + TOINT - TOINT; + if (y == 0) + return 0*x; + return y; } #endif |