From 79ca86094d70f43252b683c3a3ccb572d462cf28 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Wed, 29 Oct 2014 00:25:50 +0100 Subject: [PATCH] fix rint.c and rintf.c when FLT_EVAL_METHOD!=0 The old code used the rounding idiom incorrectly: y = (double)(x + 0x1p52) - 0x1p52; the cast is useless if FLT_EVAL_METHOD==0 and causes a second rounding if FLT_EVAL_METHOD==2 which can give incorrect result in nearest rounding mode, so the correct idiom is to add/sub a power-of-2 according to the characteristics of double_t. This did not cause actual bug because only i386 is affected where rint is implemented in asm. Other rounding functions use a similar idiom, but they give correct results because they only rely on getting a neighboring integer result and the rounding direction is fixed up separately independently of the current rounding mode. However they should be fixed to use the idiom correctly too. --- src/math/rint.c | 12 ++++++++++-- src/math/rintf.c | 14 ++++++++++++-- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/src/math/rint.c b/src/math/rint.c index 81f4e622..fbba390e 100644 --- a/src/math/rint.c +++ b/src/math/rint.c @@ -1,6 +1,14 @@ +#include #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 +static const double_t toint = 1/EPS; + double rint(double x) { union {double f; uint64_t i;} u = {x}; @@ -11,9 +19,9 @@ double rint(double x) if (e >= 0x3ff+52) return x; if (s) - y = (double)(x - 0x1p52) + 0x1p52; + y = x - toint + toint; else - y = (double)(x + 0x1p52) - 0x1p52; + y = x + toint - toint; if (y == 0) return s ? -0.0 : 0; return y; diff --git a/src/math/rintf.c b/src/math/rintf.c index 9cfc2a26..9047688d 100644 --- a/src/math/rintf.c +++ b/src/math/rintf.c @@ -1,6 +1,16 @@ +#include #include #include +#if FLT_EVAL_METHOD==0 +#define EPS FLT_EPSILON +#elif FLT_EVAL_METHOD==1 +#define EPS DBL_EPSILON +#elif FLT_EVAL_METHOD==2 +#define EPS LDBL_EPSILON +#endif +static const float_t toint = 1/EPS; + float rintf(float x) { union {float f; uint32_t i;} u = {x}; @@ -11,9 +21,9 @@ float rintf(float x) if (e >= 0x7f+23) return x; if (s) - y = (float)(x - 0x1p23f) + 0x1p23f; + y = x - toint + toint; else - y = (float)(x + 0x1p23f) - 0x1p23f; + y = x + toint - toint; if (y == 0) return s ? -0.0f : 0.0f; return y; -- 2.25.1