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 <sys/types.h> 23 #include <machine/endian.h> 24 #include <math.h> 25 26 /* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */ 27 28 #if BYTE_ORDER == BIG_ENDIAN 29 30 typedef union 31 { 32 double value; 33 struct 34 { 35 u_int32_t msw; 36 u_int32_t lsw; 37 } parts; 38 } ieee_double_shape_type; 39 40 #endif 41 42 #if BYTE_ORDER == LITTLE_ENDIAN 43 44 typedef union 45 { 46 double value; 47 struct 48 { 49 u_int32_t lsw; 50 u_int32_t msw; 51 } parts; 52 } ieee_double_shape_type; 53 54 #endif 55 56 /* Get two 32 bit ints from a double. */ 57 58 #define EXTRACT_WORDS(ix0,ix1,d) \ 59 do { \ 60 ieee_double_shape_type ew_u; \ 61 ew_u.value = (d); \ 62 (ix0) = ew_u.parts.msw; \ 63 (ix1) = ew_u.parts.lsw; \ 64 } while (0) 65 66 /* Get the more significant 32 bit int from a double. */ 67 68 #define GET_HIGH_WORD(i,d) \ 69 do { \ 70 ieee_double_shape_type gh_u; \ 71 gh_u.value = (d); \ 72 (i) = gh_u.parts.msw; \ 73 } while (0) 74 75 /* Set a double from two 32 bit ints. */ 76 77 #define INSERT_WORDS(d,ix0,ix1) \ 78 do { \ 79 ieee_double_shape_type iw_u; \ 80 iw_u.parts.msw = (ix0); \ 81 iw_u.parts.lsw = (ix1); \ 82 (d) = iw_u.value; \ 83 } while (0) 84 85 static const double one = 1.0; 86 87 double 88 modf(double x, double *iptr) 89 { 90 int32_t i0,i1,j0; 91 u_int32_t i; 92 EXTRACT_WORDS(i0,i1,x); 93 j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */ 94 if(j0<20) { /* integer part in high x */ 95 if(j0<0) { /* |x|<1 */ 96 INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */ 97 return x; 98 } else { 99 i = (0x000fffff)>>j0; 100 if(((i0&i)|i1)==0) { /* x is integral */ 101 u_int32_t high; 102 *iptr = x; 103 GET_HIGH_WORD(high,x); 104 INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ 105 return x; 106 } else { 107 INSERT_WORDS(*iptr,i0&(~i),0); 108 return x - *iptr; 109 } 110 } 111 } else if (j0>51) { /* no fraction part */ 112 u_int32_t high; 113 if (j0 == 0x400) { /* inf/NaN */ 114 *iptr = x; 115 return 0.0 / x; 116 } 117 *iptr = x*one; 118 GET_HIGH_WORD(high,x); 119 INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ 120 return x; 121 } else { /* fraction part in low x */ 122 i = ((u_int32_t)(0xffffffff))>>(j0-20); 123 if((i1&i)==0) { /* x is integral */ 124 u_int32_t high; 125 *iptr = x; 126 GET_HIGH_WORD(high,x); 127 INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ 128 return x; 129 } else { 130 INSERT_WORDS(*iptr,i0,i1&(~i)); 131 return x - *iptr; 132 } 133 } 134 } 135