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

@@ -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;

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;