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