xref: /freebsd/lib/msun/src/s_modf.c (revision 4b9d6057)
1 /*
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  */
11 
12 /*
13  * modf(double x, double *iptr)
14  * return fraction part of x, and return x's integral part in *iptr.
15  * Method:
16  *	Bit twiddling.
17  *
18  * Exception:
19  *	No exception.
20  */
21 
22 #include "math.h"
23 #include "math_private.h"
24 
25 static const double one = 1.0;
26 
27 double
28 modf(double x, double *iptr)
29 {
30 	int32_t i0,i1,j0;
31 	u_int32_t i;
32 	EXTRACT_WORDS(i0,i1,x);
33 	j0 = ((i0>>20)&0x7ff)-0x3ff;	/* exponent of x */
34 	if(j0<20) {			/* integer part in high x */
35 	    if(j0<0) {			/* |x|<1 */
36 	        INSERT_WORDS(*iptr,i0&0x80000000,0);	/* *iptr = +-0 */
37 		return x;
38 	    } else {
39 		i = (0x000fffff)>>j0;
40 		if(((i0&i)|i1)==0) {		/* x is integral */
41 		    u_int32_t high;
42 		    *iptr = x;
43 		    GET_HIGH_WORD(high,x);
44 		    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
45 		    return x;
46 		} else {
47 		    INSERT_WORDS(*iptr,i0&(~i),0);
48 		    return x - *iptr;
49 		}
50 	    }
51 	} else if (j0>51) {		/* no fraction part */
52 	    u_int32_t high;
53 	    if (j0 == 0x400) {		/* inf/NaN */
54 		*iptr = x;
55 		return 0.0 / x;
56 	    }
57 	    *iptr = x*one;
58 	    GET_HIGH_WORD(high,x);
59 	    INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
60 	    return x;
61 	} else {			/* fraction part in low x */
62 	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
63 	    if((i1&i)==0) { 		/* x is integral */
64 	        u_int32_t high;
65 		*iptr = x;
66 		GET_HIGH_WORD(high,x);
67 		INSERT_WORDS(x,high&0x80000000,0);	/* return +-0 */
68 		return x;
69 	    } else {
70 	        INSERT_WORDS(*iptr,i0,i1&(~i));
71 		return x - *iptr;
72 	    }
73 	}
74 }
75