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