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