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