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[] = "$Id: s_modf.c,v 1.4 1994/03/03 17:04:42 jtc Exp $"; 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 <machine/endian.h> 29 30 #if BYTE_ORDER == LITTLE_ENDIAN 31 #define n0 1 32 #define n1 0 33 #else 34 #define n0 0 35 #define n1 1 36 #endif 37 38 #ifdef __STDC__ 39 static const double one = 1.0; 40 #else 41 static double one = 1.0; 42 #endif 43 44 #ifdef __STDC__ 45 double modf(double x, double *iptr) 46 #else 47 double modf(x, iptr) 48 double x,*iptr; 49 #endif 50 { 51 int i0,i1,j0; 52 unsigned i; 53 i0 = *(n0+(int*)&x); /* high x */ 54 i1 = *(n1+(int*)&x); /* low x */ 55 j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */ 56 if(j0<20) { /* integer part in high x */ 57 if(j0<0) { /* |x|<1 */ 58 *(n0+(int*)iptr) = i0&0x80000000; 59 *(n1+(int*)iptr) = 0; /* *iptr = +-0 */ 60 return x; 61 } else { 62 i = (0x000fffff)>>j0; 63 if(((i0&i)|i1)==0) { /* x is integral */ 64 *iptr = x; 65 *(n0+(int*)&x) &= 0x80000000; 66 *(n1+(int*)&x) = 0; /* return +-0 */ 67 return x; 68 } else { 69 *(n0+(int*)iptr) = i0&(~i); 70 *(n1+(int*)iptr) = 0; 71 return x - *iptr; 72 } 73 } 74 } else if (j0>51) { /* no fraction part */ 75 *iptr = x*one; 76 *(n0+(int*)&x) &= 0x80000000; 77 *(n1+(int*)&x) = 0; /* return +-0 */ 78 return x; 79 } else { /* fraction part in low x */ 80 i = ((unsigned)(0xffffffff))>>(j0-20); 81 if((i1&i)==0) { /* x is integral */ 82 *iptr = x; 83 *(n0+(int*)&x) &= 0x80000000; 84 *(n1+(int*)&x) = 0; /* return +-0 */ 85 return x; 86 } else { 87 *(n0+(int*)iptr) = i0; 88 *(n1+(int*)iptr) = i1&(~i); 89 return x - *iptr; 90 } 91 } 92 } 93