From 79ca86094d70f43252b683c3a3ccb572d462cf28 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Wed, 29 Oct 2014 00:25:50 +0100 Subject: 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 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) (limited to 'src/math/rint.c') 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; -- cgit v1.2.1