>From 902e0db0a7a9d4f11870e1f784e1d14b5491ce8c Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Mon, 16 Sep 2019 20:33:11 +0000 Subject: [PATCH] math: optimize lrint on 32bit targets lrint in (LONG_MAX, 1/DBL_EPSILON) and in (-1/DBL_EPSILON, LONG_MIN) is not trivial: rounding to int may be inexact, but the conversion to int may overflow and then the inexact flag must not be raised. (the overflow threshold is rounding mode dependent). this matters on 32bit targets (without single instruction lrint or rint), so the common case (when there is no overflow) is optimized by inlining the lrint logic, otherwise the old code is kept as a fallback. (on my laptop an i486 lrint call is asm:10ns, old c:31ns, new c:22ns) --- src/math/lrint.c | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/math/lrint.c b/src/math/lrint.c index bdca8b7c..08abf326 100644 --- a/src/math/lrint.c +++ b/src/math/lrint.c @@ -1,5 +1,6 @@ #include #include +#include #include "libm.h" /* @@ -26,7 +27,18 @@ as a double. */ #if LONG_MAX < 1U<<53 && defined(FE_INEXACT) -long lrint(double x) +#include +#include +#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 +#define EPS DBL_EPSILON +#elif FLT_EVAL_METHOD==2 +#define EPS LDBL_EPSILON +#endif +#ifdef __GNUC__ +/* avoid stack frame in lrint */ +__attribute__((noinline)) +#endif +static long lrint_slow(double x) { #pragma STDC FENV_ACCESS ON int e; @@ -38,6 +50,24 @@ long lrint(double x) /* conversion */ return x; } + +long lrint(double x) +{ + union {double f; uint64_t i;} u = {x}; + uint32_t abstop = u.i>>32 & 0x7fffffff; + + if (abstop < 0x41dfffff) { + /* |x| < 0x7ffffc00 */ + double_t y, t; + union {double f; uint64_t i;} toint = {1/EPS}; + toint.i |= u.i & (1ULL << 63); + t = toint.f; + y = x + t - t; + /* conversion */ + return y; + } + return lrint_slow(x); +} #else long lrint(double x) { -- 2.21.0