#include <math.h>
#include <float.h>
#include <errno.h>

/* modf suitable for IEEE double-precision */

#define	MASK	0x7ffL
#define SIGN	0x80000000
#define	SHIFT	20
#define	BIAS	1022L

#ifdef OLD_WAY
typedef	union
{
	double	x;
	struct
	{
		long	hi;		/* big-endian */
		long	lo;
	};
} Cheat;
#endif

double
modf(double d, double *ip)
{
	FPdbleword x;
	int e;

	if(-1 < d && d < 1) {
		*ip = 0;
		return d;
	}
	x.x = d;
	x.hi &= ~SIGN;
	e = (x.hi >> SHIFT) & MASK;
	if(e == MASK || e == 0){
		errno = EDOM;
		*ip = (d > 0)? HUGE_VAL : -HUGE_VAL;
		return 0;
	}
	e -= BIAS;
	if(e <= SHIFT+1) {
		x.hi &= ~(0x1fffffL >> e);
		x.lo = 0;
	} else
	if(e <= SHIFT+33)
		x.lo &= ~(0x7fffffffL >> (e-SHIFT-2));
	if(d > 0){
		*ip = x.x;
		return d - x.x;
	}else{
		*ip = -x.x;
		return d + x.x;
	}
}