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