faster lrint and llrint functions

A faster workaround for spurious inexact exceptions
when the result cannot be represented. The old code
actually could be wrong, because gcc reordered the
integer conversion and the exception check.

diff --git a/src/math/llrint.c b/src/math/llrint.c
index c0a40721c45e2ba6760285312fed28772b166726..ee783b8e403475075715bfbdf3ee514b51fb00cf 100644 (file)
--- 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 f06a3c278eaf91e00f6dabd89ee7b7652fc27fb1..e41b6d41862304f09e6fe00c0cc788ade8940065 100644 (file)
--- 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 ec2cf86b5251969a7fcaeb0c896f521f2fe24735..f1cc47ed453c82d07e40554f11b9f5cfea53b83b 100644 (file)
--- 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 9754fa749f7cc5e13a67e2e0f28e2fe2cc6c6e44..feba28d072c153e50dac0fe207e833e2bc7361a4 100644 (file)
--- 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 caed7ca5e4f56a5019f4cb76190e1e971b33cf13..34d1081c0fbf1a1fe98cb94f158f7036d45605ec 100644 (file)
--- 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 c076326f70b5cd32d779b4ec3b7e9f42f9cc49e7..0e579bc51c9a4066e7a0b197614c6f9ad6a458b2 100644 (file)
--- 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