diff options
| -rw-r--r-- | src/math/llrint.c | 12 | ||||
| -rw-r--r-- | src/math/llrintf.c | 12 | ||||
| -rw-r--r-- | src/math/llrintl.c | 28 | ||||
| -rw-r--r-- | src/math/lrint.c | 87 | ||||
| -rw-r--r-- | src/math/lrintf.c | 12 | ||||
| -rw-r--r-- | src/math/lrintl.c | 28 |
6 files changed, 99 insertions, 80 deletions
diff --git a/src/math/llrint.c b/src/math/llrint.c index c0a40721..ee783b8e 100644 --- a/src/math/llrint.c +++ b/src/math/llrint.c @@ -1,8 +1,8 @@ -#define type double -#define roundit rint -#define dtype long long -#define fn llrint - -#include "lrint.c" +#include <math.h> +/* assumes LLONG_MAX > 2^53, see comments in lrint.c */ +long long llrint(double x) +{ + return rint(x); +} diff --git a/src/math/llrintf.c b/src/math/llrintf.c index f06a3c27..e41b6d41 100644 --- a/src/math/llrintf.c +++ b/src/math/llrintf.c @@ -1,6 +1,8 @@ -#define type float -#define roundit rintf -#define dtype long long -#define fn llrintf +#include <math.h> -#include "lrint.c" +/* assumes LLONG_MAX > 2^24, see comments in lrint.c */ + +long long llrintf(float x) +{ + return rintf(x); +} diff --git a/src/math/llrintl.c b/src/math/llrintl.c index ec2cf86b..f1cc47ed 100644 --- a/src/math/llrintl.c +++ b/src/math/llrintl.c @@ -1,5 +1,7 @@ -#include <math.h> -#include <float.h> +#include <limits.h> +#include <fenv.h> +#include "libm.h" + #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long long llrintl(long double x) @@ -7,10 +9,22 @@ long long llrintl(long double x) return llrint(x); } #else -#define type long double -#define roundit rintl -#define dtype long long -#define fn llrintl +/* +see comments in lrint.c -#include "lrint.c" +Note that if LLONG_MAX == 0x7fffffffffffffff && LDBL_MANT_DIG == 64 +then x == 2**63 - 0.5 is the only input that overflows and +raises inexact (with tonearest or upward rounding mode) +*/ +long long llrintl(long double x) +{ + int e; + + e = fetestexcept(FE_INEXACT); + x = rintl(x); + if (!e && (x > LLONG_MAX || x < LLONG_MIN)) + feclearexcept(FE_INEXACT); + /* conversion */ + return x; +} #endif diff --git a/src/math/lrint.c b/src/math/lrint.c index 9754fa74..feba28d0 100644 --- a/src/math/lrint.c +++ b/src/math/lrint.c @@ -1,58 +1,45 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_lrint.c */ -/*- - * Copyright (c) 2005 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 <limits.h> #include <fenv.h> #include "libm.h" -#ifndef type -#define type double -#define roundit rint -#define dtype long -#define fn lrint -#endif - /* - * C99 says we should not raise a spurious inexact exception when an - * invalid exception is raised. Unfortunately, the set of inputs - * that overflows depends on the rounding mode when 'dtype' has more - * significant bits than 'type'. Hence, we bend over backwards for the - * sake of correctness; an MD implementation could be more efficient. - */ -dtype fn(type x) +If the result cannot be represented (overflow, nan), then +lrint raises the invalid exception. + +Otherwise if the input was not an integer then the inexact +exception is raised. + +C99 is a bit vague about whether inexact exception is +allowed to be raised when invalid is raised. +(F.9 explicitly allows spurious inexact exceptions, F.9.6.5 +does not make it clear if that rule applies to lrint, but +IEEE 754r 7.8 seems to forbid spurious inexact exception in +the ineger conversion functions) + +So we try to make sure that no spurious inexact exception is +raised in case of an overflow. + +If the bit size of long > precision of double, then there +cannot be inexact rounding in case the result overflows, +otherwise LONG_MAX and LONG_MIN can be represented exactly +as a double. +*/ + +#if LONG_MAX < 1U<<53 +long lrint(double x) { - fenv_t env; - dtype d; + int e; - feholdexcept(&env); - d = (dtype)roundit(x); -#if defined(FE_INVALID) && defined(FE_INEXACT) - if (fetestexcept(FE_INVALID)) + e = fetestexcept(FE_INEXACT); + x = rint(x); + if (!e && (x > LONG_MAX || x < LONG_MIN)) feclearexcept(FE_INEXACT); -#endif - feupdateenv(&env); - return d; + /* conversion */ + return x; } +#else +long lrint(double x) +{ + return rint(x); +} +#endif diff --git a/src/math/lrintf.c b/src/math/lrintf.c index caed7ca5..34d1081c 100644 --- a/src/math/lrintf.c +++ b/src/math/lrintf.c @@ -1,6 +1,8 @@ -#define type float -#define roundit rintf -#define dtype long -#define fn lrintf +#include <math.h> -#include "lrint.c" +/* assumes LONG_MAX > 2^24, see comments in lrint.c */ + +long lrintf(float x) +{ + return rintf(x); +} diff --git a/src/math/lrintl.c b/src/math/lrintl.c index c076326f..0e579bc5 100644 --- a/src/math/lrintl.c +++ b/src/math/lrintl.c @@ -1,5 +1,7 @@ -#include <math.h> -#include <float.h> +#include <limits.h> +#include <fenv.h> +#include "libm.h" + #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long lrintl(long double x) @@ -7,10 +9,22 @@ long lrintl(long double x) return lrint(x); } #else -#define type long double -#define roundit rintl -#define dtype long -#define fn lrintl +/* +see comments in lrint.c -#include "lrint.c" +Note that if LONG_MAX == 0x7fffffffffffffff && LDBL_MANT_DIG == 64 +then x == 2**63 - 0.5 is the only input that overflows and +raises inexact (with tonearest or upward rounding mode) +*/ +long lrintl(long double x) +{ + int e; + + e = fetestexcept(FE_INEXACT); + x = rintl(x); + if (!e && (x > LONG_MAX || x < LONG_MIN)) + feclearexcept(FE_INEXACT); + /* conversion */ + return x; +} #endif |