1#include "libm.h"
2
3#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4long double modfl(long double x, long double* iptr) {
5    double d;
6    long double r;
7
8    r = modf(x, &d);
9    *iptr = d;
10    return r;
11}
12#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
13
14static const long double toint = 1 / LDBL_EPSILON;
15
16long double modfl(long double x, long double* iptr) {
17    union ldshape u = {x};
18    int e = (u.i.se & 0x7fff) - 0x3fff;
19    int s = u.i.se >> 15;
20    long double absx;
21    long double y;
22
23    /* no fractional part */
24    if (e >= LDBL_MANT_DIG - 1) {
25        *iptr = x;
26        if (isnan(x))
27            return x;
28        return s ? -0.0 : 0.0;
29    }
30
31    /* no integral part*/
32    if (e < 0) {
33        *iptr = s ? -0.0 : 0.0;
34        return x;
35    }
36
37    /* raises spurious inexact */
38    absx = s ? -x : x;
39    y = absx + toint - toint - absx;
40    if (y == 0) {
41        *iptr = x;
42        return s ? -0.0 : 0.0;
43    }
44    if (y > 0)
45        y -= 1;
46    if (s)
47        y = -y;
48    *iptr = x + y;
49    return -y;
50}
51#endif
52