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.

diff --git a/src/math/rint.c b/src/math/rint.c
index 81f4e6223b748b59db415725a9b07a8969dfeac4..fbba390e7d723774d41ba2a9da74957c20643d63 100644 (file)
--- a/src/math/rint.c
+++ b/src/math/rint.c
@@ -1,6 +1,14 @@
+#include <float.h>
 #include <math.h>
 #include <stdint.h>
 
+#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 9cfc2a261f3c7f758b776859691c65e1bf6654bd..9047688d246a69c2698ce27399135e1db7ac2f65 100644 (file)
--- a/src/math/rintf.c
+++ b/src/math/rintf.c
@@ -1,6 +1,16 @@
+#include <float.h>
 #include <math.h>
 #include <stdint.h>
 
+#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;