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