Home | History | Annotate | Download | only in fdlibm
      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__
     27 static const double one = 1.0;
     28 #else
     29 static 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