1#include "libm.h"
2
3double modf(double x, double* iptr) {
4    union {
5        double f;
6        uint64_t i;
7    } u = {x};
8    uint64_t mask;
9    int e = (int)(u.i >> 52 & 0x7ff) - 0x3ff;
10
11    /* no fractional part */
12    if (e >= 52) {
13        *iptr = x;
14        if (e == 0x400 && u.i << 12 != 0) /* nan */
15            return x;
16        u.i &= 1ULL << 63;
17        return u.f;
18    }
19
20    /* no integral part*/
21    if (e < 0) {
22        u.i &= 1ULL << 63;
23        *iptr = u.f;
24        return x;
25    }
26
27    mask = -1ULL >> 12 >> e;
28    if ((u.i & mask) == 0) {
29        *iptr = x;
30        u.i &= 1ULL << 63;
31        return u.f;
32    }
33    u.i &= ~mask;
34    *iptr = u.f;
35    return x - u.f;
36}
37