summaryrefslogtreecommitdiff
path: root/src/math/floor.c
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2014-10-29 00:34:37 +0100
committerRich Felker <dalias@aerifal.cx>2014-10-31 11:35:40 -0400
commit0ce946cf808274c2d6e5419b139e130c8ad4bd30 (patch)
treee6614e756dde0afbcd48e38916d0208fed93ece1 /src/math/floor.c
parent79ca86094d70f43252b683c3a3ccb572d462cf28 (diff)
downloadmusl-0ce946cf808274c2d6e5419b139e130c8ad4bd30.tar.gz
musl-0ce946cf808274c2d6e5419b139e130c8ad4bd30.tar.bz2
musl-0ce946cf808274c2d6e5419b139e130c8ad4bd30.tar.xz
musl-0ce946cf808274c2d6e5419b139e130c8ad4bd30.zip
math: use the rounding idiom consistently
the idiomatic rounding of x is n = x + toint - toint; where toint is either 1/EPSILON (x is non-negative) or 1.5/EPSILON (x may be negative and nearest rounding mode is assumed) and EPSILON is according to the evaluation precision (the type of toint is not very important, because single precision float can represent the 1/EPSILON of ieee binary128). in case of FLT_EVAL_METHOD!=0 this avoids a useless store to double or float precision, and the long double code became cleaner with 1/LDBL_EPSILON instead of ifdefs for toint. __rem_pio2f and __rem_pio2 functions slightly changed semantics: on i386 a double-rounding is avoided so close to half-way cases may get evaluated differently eg. as sin(pi/4-eps) instead of cos(pi/4+eps)
Diffstat (limited to 'src/math/floor.c')
-rw-r--r--src/math/floor.c11
1 files changed, 9 insertions, 2 deletions
diff --git a/src/math/floor.c b/src/math/floor.c
index ebc9fabe..14a31cd8 100644
--- a/src/math/floor.c
+++ b/src/math/floor.c
@@ -1,5 +1,12 @@
#include "libm.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 floor(double x)
{
union {double f; uint64_t i;} u = {x};
@@ -10,9 +17,9 @@ double floor(double x)
return x;
/* y = int(x) - x, where int(x) is an integer neighbor of x */
if (u.i >> 63)
- y = (double)(x - 0x1p52) + 0x1p52 - x;
+ y = x - toint + toint - x;
else
- y = (double)(x + 0x1p52) - 0x1p52 - x;
+ y = x + toint - toint - x;
/* special case because of non-nearest rounding modes */
if (e <= 0x3ff-1) {
FORCE_EVAL(y);