1*2813126eSlogin 2*2813126eSlogin 3*2813126eSlogin #include "libm.h" 4*2813126eSlogin 5*2813126eSlogin #if __FLT_EVAL_METHOD__ == 0 || __FLT_EVAL_METHOD__ == 1 6*2813126eSlogin #define EPS __DBL_EPSILON__ 7*2813126eSlogin #elif __FLT_EVAL_METHOD__ == 2 8*2813126eSlogin #define EPS __LDBL_EPSILON__ 9*2813126eSlogin #endif 10*2813126eSlogin static const double toint = 1 / EPS; 11*2813126eSlogin round(double x)12*2813126eSlogindouble round(double x) 13*2813126eSlogin { 14*2813126eSlogin union 15*2813126eSlogin { 16*2813126eSlogin double f; 17*2813126eSlogin uint64_t i; 18*2813126eSlogin } u = {x}; 19*2813126eSlogin 20*2813126eSlogin int e = u.i >> 52 & 0x7ff; 21*2813126eSlogin double y; 22*2813126eSlogin 23*2813126eSlogin if (e >= 0x3ff + 52) 24*2813126eSlogin return x; 25*2813126eSlogin if (u.i >> 63) 26*2813126eSlogin x = -x; 27*2813126eSlogin if (e < 0x3ff - 1) 28*2813126eSlogin { 29*2813126eSlogin /* raise inexact if x!=0 */ 30*2813126eSlogin FORCE_EVAL(x + toint); 31*2813126eSlogin return 0 * u.f; 32*2813126eSlogin } 33*2813126eSlogin y = x + toint - toint - x; 34*2813126eSlogin if (y > 0.5) 35*2813126eSlogin y = y + x - 1; 36*2813126eSlogin else if (y <= -0.5) 37*2813126eSlogin y = y + x + 1; 38*2813126eSlogin else 39*2813126eSlogin y = y + x; 40*2813126eSlogin if (u.i >> 63) 41*2813126eSlogin y = -y; 42*2813126eSlogin return y; 43*2813126eSlogin }