1*f504f610SAugustin Cavalier #include "libm.h" 2*f504f610SAugustin Cavalier 3*f504f610SAugustin Cavalier #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 4*f504f610SAugustin Cavalier #define EPS DBL_EPSILON 5*f504f610SAugustin Cavalier #elif FLT_EVAL_METHOD==2 6*f504f610SAugustin Cavalier #define EPS LDBL_EPSILON 7*f504f610SAugustin Cavalier #endif 8*f504f610SAugustin Cavalier static const double_t toint = 1/EPS; 9*f504f610SAugustin Cavalier ceil(double x)10*f504f610SAugustin Cavalierdouble ceil(double x) 11*f504f610SAugustin Cavalier { 12*f504f610SAugustin Cavalier union {double f; uint64_t i;} u = {x}; 13*f504f610SAugustin Cavalier int e = u.i >> 52 & 0x7ff; 14*f504f610SAugustin Cavalier double_t y; 15*f504f610SAugustin Cavalier 16*f504f610SAugustin Cavalier if (e >= 0x3ff+52 || x == 0) 17*f504f610SAugustin Cavalier return x; 18*f504f610SAugustin Cavalier /* y = int(x) - x, where int(x) is an integer neighbor of x */ 19*f504f610SAugustin Cavalier if (u.i >> 63) 20*f504f610SAugustin Cavalier y = x - toint + toint - x; 21*f504f610SAugustin Cavalier else 22*f504f610SAugustin Cavalier y = x + toint - toint - x; 23*f504f610SAugustin Cavalier /* special case because of non-nearest rounding modes */ 24*f504f610SAugustin Cavalier if (e <= 0x3ff-1) { 25*f504f610SAugustin Cavalier FORCE_EVAL(y); 26*f504f610SAugustin Cavalier return u.i >> 63 ? -0.0 : 1; 27*f504f610SAugustin Cavalier } 28*f504f610SAugustin Cavalier if (y < 0) 29*f504f610SAugustin Cavalier return x + y + 1; 30*f504f610SAugustin Cavalier return x + y; 31*f504f610SAugustin Cavalier } 32