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