musl - musl - an implementation of the standard library for Linux-based systems

index : musl
musl - an implementation of the standard library for Linux-based systems
summary refs log tree commit diff
path: root/src/math/rint.c
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2014年10月29日 00:25:50 +0100
committerRich Felker <dalias@aerifal.cx>2014年10月31日 11:35:40 -0400
commit79ca86094d70f43252b683c3a3ccb572d462cf28 (patch)
tree3396f9fcd33336c7cecb24399e716c407077acab /src/math/rint.c
parent2da3ab1382ca8e39eb1e4428103764a81fba73d3 (diff)
downloadmusl-79ca86094d70f43252b683c3a3ccb572d462cf28.tar.gz
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.
Diffstat (limited to 'src/math/rint.c')
-rw-r--r--src/math/rint.c 12
1 files changed, 10 insertions, 2 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 <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;
generated by cgit v1.2.1 (git 2.18.0) at 2025年10月04日 00:39:41 +0000

AltStyle によって変換されたページ (->オリジナル) /