diff options
| author | Szabolcs Nagy <nsz@port70.net> | 2014-10-29 00:25:50 +0100 |
|---|---|---|
| committer | Rich Felker <dalias@aerifal.cx> | 2014-10-31 11:35:40 -0400 |
| commit | 79ca86094d70f43252b683c3a3ccb572d462cf28 (patch) | |
| tree | 3396f9fcd33336c7cecb24399e716c407077acab | |
| parent | 2da3ab1382ca8e39eb1e4428103764a81fba73d3 (diff) | |
| download | musl-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.
| -rw-r--r-- | src/math/rint.c | 12 | ||||
| -rw-r--r-- | 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 <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 9cfc2a26..9047688d 100644 --- 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; |