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 /src | |
| parent | 2da3ab1382ca8e39eb1e4428103764a81fba73d3 (diff) | |
| download | musl-79ca86094d70f43252b683c3a3ccb572d462cf28.tar.gz musl-79ca86094d70f43252b683c3a3ccb572d462cf28.tar.bz2 musl-79ca86094d70f43252b683c3a3ccb572d462cf28.tar.xz musl-79ca86094d70f43252b683c3a3ccb572d462cf28.zip | |
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')
| -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; |