1
2/* @(#)s_modf.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14/*
15 * ieee_modf(double x, double *iptr)
16 * return fraction part of x, and return x's integral part in *iptr.
17 * Method:
18 *	Bit twiddling.
19 *
20 * Exception:
21 *	No exception.
22 */
23
24#include "fdlibm.h"
25
26#ifdef __STDC__
27static const double one = 1.0;
28#else
29static double one = 1.0;
30#endif
31
32#ifdef __STDC__
33	double ieee_modf(double x, double *iptr)
34#else
35	double ieee_modf(x, iptr)
36	double x,*iptr;
37#endif
38{
39	int i0,i1,j0;
40	unsigned i;
41	i0 =  __HI(x);		/* high x */
42	i1 =  __LO(x);		/* low  x */
43	j0 = ((i0>>20)&0x7ff)-0x3ff;	/* exponent of x */
44	if(j0<20) {			/* integer part in high x */
45	    if(j0<0) {			/* |x|<1 */
46		__HIp(iptr) = i0&0x80000000;
47		__LOp(iptr) = 0;		/* *iptr = +-0 */
48		return x;
49	    } else {
50		i = (0x000fffff)>>j0;
51		if(((i0&i)|i1)==0) {		/* x is integral */
52		    *iptr = x;
53		    __HI(x) &= 0x80000000;
54		    __LO(x)  = 0;	/* return +-0 */
55		    return x;
56		} else {
57		    __HIp(iptr) = i0&(~i);
58		    __LOp(iptr) = 0;
59		    return x - *iptr;
60		}
61	    }
62	} else if (j0>51) {		/* no fraction part */
63	    *iptr = x*one;
64	    __HI(x) &= 0x80000000;
65	    __LO(x)  = 0;	/* return +-0 */
66	    return x;
67	} else {			/* fraction part in low x */
68	    i = ((unsigned)(0xffffffff))>>(j0-20);
69	    if((i1&i)==0) { 		/* x is integral */
70		*iptr = x;
71		__HI(x) &= 0x80000000;
72		__LO(x)  = 0;	/* return +-0 */
73		return x;
74	    } else {
75		__HIp(iptr) = i0;
76		__LOp(iptr) = i1&(~i);
77		return x - *iptr;
78	    }
79	}
80}
81