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)

diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c
index 5fafc4a8aca547ebd671f40680618960a61705b0..a40db9fc05fdd028d90f11d8f9ce8f8db6992fcf 100644 (file)
--- a/src/math/__rem_pio2.c
+++ b/src/math/__rem_pio2.c
@@ -19,6 +19,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
+
 /*
  * invpio2:  53 bits of 2/pi
  * pio2_1:   first  33 bit of pi/2
@@ -29,6 +35,7 @@
  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  */
 static const double
+toint   = 1.5/EPS,
 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 pio2_1  = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
@@ -41,8 +48,8 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 int __rem_pio2(double x, double *y)
 {
        union {double f; uint64_t i;} u = {x};
-       double_t z,w,t,r;
-       double tx[3],ty[2],fn;
+       double_t z,w,t,r,fn;
+       double tx[3],ty[2];
        uint32_t ix;
        int sign, n, ex, ey, i;
 
@@ -111,8 +118,7 @@ int __rem_pio2(double x, double *y)
        if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
 medium:
                /* rint(x/(pi/2)), Assume round-to-nearest. */
-               fn = x*invpio2 + 0x1.8p52;
-               fn = fn - 0x1.8p52;
+               fn = x*invpio2 + toint - toint;
                n = (int32_t)fn;
                r = x - fn*pio2_1;
                w = fn*pio2_1t;  /* 1st round, good to 85 bits */

diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c
index 838e1fcac1921da9164a79450eafc460ad0e4d3e..f5163856512880fdf82597b169e0b75d6ae0257c 100644 (file)
--- a/src/math/__rem_pio2f.c
+++ b/src/math/__rem_pio2f.c
@@ -22,12 +22,19 @@
 
 #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
+
 /*
  * invpio2:  53 bits of 2/pi
  * pio2_1:   first 25 bits of pi/2
  * pio2_1t:  pi/2 - pio2_1
  */
 static const double
+toint   = 1.5/EPS,
 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 pio2_1  = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
 pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
@@ -35,7 +42,8 @@ pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
 int __rem_pio2f(float x, double *y)
 {
        union {float f; uint32_t i;} u = {x};
-       double tx[1],ty[1],fn;
+       double tx[1],ty[1];
+       double_t fn;
        uint32_t ix;
        int n, sign, e0;
 
@@ -43,8 +51,7 @@ int __rem_pio2f(float x, double *y)
        /* 25+53 bit pi is good enough for medium size */
        if (ix < 0x4dc90fdb) {  /* |x| ~< 2^28*(pi/2), medium size */
                /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-               fn = x*invpio2 + 0x1.8p52;
-               fn = fn - 0x1.8p52;
+               fn = x*invpio2 + toint - toint;
                n  = (int32_t)fn;
                *y = x - fn*pio2_1 - fn*pio2_1t;
                return n;

diff --git a/src/math/__rem_pio2l.c b/src/math/__rem_pio2l.c
index 8b15b7b2a3d802192772fa2ea34dd583cfa2995f..77255bd80ae480004dfe51490a2091f746d9a00f 100644 (file)
--- a/src/math/__rem_pio2l.c
+++ b/src/math/__rem_pio2l.c
@@ -20,10 +20,11 @@
  * use __rem_pio2_large() for large x
  */
 
+static const long double toint = 1.5/LDBL_EPSILON;
+
 #if LDBL_MANT_DIG == 64
 /* u ~< 0x1p25*pi/2 */
 #define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.m>>48) < ((0x3fff + 25)<<16 | 0x921f>>1 | 0x8000))
-#define TOINT 0x1.8p63
 #define QUOBITS(x) ((uint32_t)(int32_t)x & 0x7fffffff)
 #define ROUND1 22
 #define ROUND2 61
@@ -50,7 +51,6 @@ pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
 #elif LDBL_MANT_DIG == 113
 /* u ~< 0x1p45*pi/2 */
 #define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.top) < ((0x3fff + 45)<<16 | 0x921f))
-#define TOINT 0x1.8p112
 #define QUOBITS(x) ((uint32_t)(int64_t)x & 0x7fffffff)
 #define ROUND1 51
 #define ROUND2 119
@@ -77,7 +77,7 @@ int __rem_pio2l(long double x, long double *y)
        ex = u.i.se & 0x7fff;
        if (SMALL(u)) {
                /* rint(x/(pi/2)), Assume round-to-nearest. */
-               fn = x*invpio2 + TOINT - TOINT;
+               fn = x*invpio2 + toint - toint;
                n = QUOBITS(fn);
                r = x-fn*pio2_1;
                w = fn*pio2_1t;  /* 1st round good to 102/180 bits (ld80/ld128) */

diff --git a/src/math/ceil.c b/src/math/ceil.c
index 22dc224c5f38b1fdebbdba03d36bf0f408298baa..b13e6f2d63689818ff349f1327fb5d1bb25eb49d 100644 (file)
--- a/src/math/ceil.c
+++ b/src/math/ceil.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 ceil(double x)
 {
        union {double f; uint64_t i;} u = {x};
@@ -10,9 +17,9 @@ double ceil(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);

diff --git a/src/math/ceill.c b/src/math/ceill.c
index a2cb0a7f999607e8b55a3582a040497618553c65..60a83020dd1cb6f6764d0ddb5b64ba9476e71d15 100644 (file)
--- a/src/math/ceill.c
+++ b/src/math/ceill.c
@@ -6,11 +6,9 @@ long double ceill(long double x)
        return ceil(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double ceill(long double x)
 {
        union ldshape u = {x};
@@ -21,9 +19,9 @@ long double ceill(long double x)
                return x;
        /* y = int(x) - x, where int(x) is an integer neighbor of x */
        if (u.i.se >> 15)
-               y = x - TOINT + TOINT - x;
+               y = x - toint + toint - x;
        else
-               y = x + TOINT - TOINT - x;
+               y = x + toint - toint - x;
        /* special case because of non-nearest rounding modes */
        if (e <= 0x3fff-1) {
                FORCE_EVAL(y);

diff --git a/src/math/floor.c b/src/math/floor.c
index ebc9fabe908223d54522c73546ac314df59922a8..14a31cd8c4c549750205369a6fe9e71eec8f71c5 100644 (file)
--- 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);

diff --git a/src/math/floorl.c b/src/math/floorl.c
index 961f9e890b7362655a1c24938cf02e61e21e3d75..16aaec48eea8e4e28e747122d6179d95e279981f 100644 (file)
--- a/src/math/floorl.c
+++ b/src/math/floorl.c
@@ -6,11 +6,9 @@ long double floorl(long double x)
        return floor(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double floorl(long double x)
 {
        union ldshape u = {x};
@@ -21,9 +19,9 @@ long double floorl(long double x)
                return x;
        /* y = int(x) - x, where int(x) is an integer neighbor of x */
        if (u.i.se >> 15)
-               y = x - TOINT + TOINT - x;
+               y = x - toint + toint - x;
        else
-               y = x + TOINT - TOINT - x;
+               y = x + toint - toint - x;
        /* special case because of non-nearest rounding modes */
        if (e <= 0x3fff-1) {
                FORCE_EVAL(y);

diff --git a/src/math/modfl.c b/src/math/modfl.c
index 4b03a4be0aa462b915c9ad5c85144771a706fcdc..a47b1924f75d1125b712e3ccf5b0f39af27fd8f2 100644 (file)
--- a/src/math/modfl.c
+++ b/src/math/modfl.c
@@ -11,11 +11,9 @@ long double modfl(long double x, long double *iptr)
        return r;
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double modfl(long double x, long double *iptr)
 {
        union ldshape u = {x};
@@ -40,7 +38,7 @@ long double modfl(long double x, long double *iptr)
 
        /* raises spurious inexact */
        absx = s ? -x : x;
-       y = absx + TOINT - TOINT - absx;
+       y = absx + toint - toint - absx;
        if (y == 0) {
                *iptr = x;
                return s ? -0.0 : 0.0;

diff --git a/src/math/rintl.c b/src/math/rintl.c
index 267250737f0a8e3b3924ee7bb662bbce2bbf6a9e..374327db2282d450ceb79e12286db767eb452148 100644 (file)
--- a/src/math/rintl.c
+++ b/src/math/rintl.c
@@ -6,11 +6,9 @@ long double rintl(long double x)
        return rint(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double rintl(long double x)
 {
        union ldshape u = {x};
@@ -21,9 +19,9 @@ long double rintl(long double x)
        if (e >= 0x3fff+LDBL_MANT_DIG-1)
                return x;
        if (s)
-               y = x - TOINT + TOINT;
+               y = x - toint + toint;
        else
-               y = x + TOINT - TOINT;
+               y = x + toint - toint;
        if (y == 0)
                return 0*x;
        return y;

diff --git a/src/math/round.c b/src/math/round.c
index 4b38d1fdbb2adca6f1f89a8ff86dc29e9206f272..130d58d2571e77cc15a67b355010e11c5ebebe91 100644 (file)
--- a/src/math/round.c
+++ b/src/math/round.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 round(double x)
 {
        union {double f; uint64_t i;} u = {x};
@@ -12,10 +19,10 @@ double round(double x)
                x = -x;
        if (e < 0x3ff-1) {
                /* raise inexact if x!=0 */
-               FORCE_EVAL(x + 0x1p52);
+               FORCE_EVAL(x + toint);
                return 0*u.f;
        }
-       y = (double)(x + 0x1p52) - 0x1p52 - x;
+       y = x + toint - toint - x;
        if (y > 0.5)
                y = y + x - 1;
        else if (y <= -0.5)

diff --git a/src/math/roundf.c b/src/math/roundf.c
index c6b277979038e0abd5ddc566342e674c8d570810..e8210af5621aa392543616d308200ddb22523281 100644 (file)
--- a/src/math/roundf.c
+++ b/src/math/roundf.c
@@ -1,5 +1,14 @@
 #include "libm.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 roundf(float x)
 {
        union {float f; uint32_t i;} u = {x};
@@ -11,10 +20,10 @@ float roundf(float x)
        if (u.i >> 31)
                x = -x;
        if (e < 0x7f-1) {
-               FORCE_EVAL(x + 0x1p23f);
+               FORCE_EVAL(x + toint);
                return 0*u.f;
        }
-       y = (float)(x + 0x1p23f) - 0x1p23f - x;
+       y = x + toint - toint - x;
        if (y > 0.5f)
                y = y + x - 1;
        else if (y <= -0.5f)

diff --git a/src/math/roundl.c b/src/math/roundl.c
index 8f3f28b348e131fc2bc7f6c3098e398c75f75a67..f4ff6820a982276321f6bfe8198e0975c7016017 100644 (file)
--- a/src/math/roundl.c
+++ b/src/math/roundl.c
@@ -6,11 +6,9 @@ long double roundl(long double x)
        return round(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double roundl(long double x)
 {
        union ldshape u = {x};
@@ -22,10 +20,10 @@ long double roundl(long double x)
        if (u.i.se >> 15)
                x = -x;
        if (e < 0x3fff-1) {
-               FORCE_EVAL(x + TOINT);
+               FORCE_EVAL(x + toint);
                return 0*u.f;
        }
-       y = x + TOINT - TOINT - x;
+       y = x + toint - toint - x;
        if (y > 0.5)
                y = y + x - 1;
        else if (y <= -0.5)

diff --git a/src/math/truncl.c b/src/math/truncl.c
index 3eedb0839713144cf1fbb2f2e05cd7ba6f04aebc..f07b1934097bcbb932e37a8b0710c99738993e39 100644 (file)
--- a/src/math/truncl.c
+++ b/src/math/truncl.c
@@ -6,11 +6,9 @@ long double truncl(long double x)
        return trunc(x);
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-#if LDBL_MANT_DIG == 64
-#define TOINT 0x1p63
-#elif LDBL_MANT_DIG == 113
-#define TOINT 0x1p112
-#endif
+
+static const long double toint = 1/LDBL_EPSILON;
+
 long double truncl(long double x)
 {
        union ldshape u = {x};
@@ -27,7 +25,7 @@ long double truncl(long double x)
        /* y = int(|x|) - |x|, where int(|x|) is an integer neighbor of |x| */
        if (s)
                x = -x;
-       y = x + TOINT - TOINT - x;
+       y = x + toint - toint - x;
        if (y > 0)
                y -= 1;
        x += y;