12116Sjkh/* @(#)s_modf.c 5.1 93/09/24 */ 22116Sjkh/* 32116Sjkh * ==================================================== 42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 52116Sjkh * 62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 72116Sjkh * Permission to use, copy, modify, and distribute this 88870Srgrimes * software is freely granted, provided that this notice 92116Sjkh * is preserved. 102116Sjkh * ==================================================== 112116Sjkh */ 122116Sjkh 13226606Sdas#include <sys/cdefs.h> 14226606Sdas__FBSDID("$FreeBSD$"); 152116Sjkh 162116Sjkh/* 178870Srgrimes * modf(double x, double *iptr) 182116Sjkh * return fraction part of x, and return x's integral part in *iptr. 192116Sjkh * Method: 202116Sjkh * Bit twiddling. 212116Sjkh * 222116Sjkh * Exception: 232116Sjkh * No exception. 242116Sjkh */ 252116Sjkh 26226606Sdas#include <sys/types.h> 27226606Sdas#include <machine/endian.h> 28226606Sdas#include <math.h> 292116Sjkh 30226606Sdas/* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */ 31226606Sdas 32226606Sdas#if BYTE_ORDER == BIG_ENDIAN 33226606Sdas 34226606Sdastypedef union 35226606Sdas{ 36226606Sdas double value; 37226606Sdas struct 38226606Sdas { 39226606Sdas u_int32_t msw; 40226606Sdas u_int32_t lsw; 41226606Sdas } parts; 42226606Sdas} ieee_double_shape_type; 43226606Sdas 44226606Sdas#endif 45226606Sdas 46226606Sdas#if BYTE_ORDER == LITTLE_ENDIAN 47226606Sdas 48226606Sdastypedef union 49226606Sdas{ 50226606Sdas double value; 51226606Sdas struct 52226606Sdas { 53226606Sdas u_int32_t lsw; 54226606Sdas u_int32_t msw; 55226606Sdas } parts; 56226606Sdas} ieee_double_shape_type; 57226606Sdas 58226606Sdas#endif 59226606Sdas 60226606Sdas/* Get two 32 bit ints from a double. */ 61226606Sdas 62226606Sdas#define EXTRACT_WORDS(ix0,ix1,d) \ 63226606Sdasdo { \ 64226606Sdas ieee_double_shape_type ew_u; \ 65226606Sdas ew_u.value = (d); \ 66226606Sdas (ix0) = ew_u.parts.msw; \ 67226606Sdas (ix1) = ew_u.parts.lsw; \ 68226606Sdas} while (0) 69226606Sdas 70226606Sdas/* Get the more significant 32 bit int from a double. */ 71226606Sdas 72226606Sdas#define GET_HIGH_WORD(i,d) \ 73226606Sdasdo { \ 74226606Sdas ieee_double_shape_type gh_u; \ 75226606Sdas gh_u.value = (d); \ 76226606Sdas (i) = gh_u.parts.msw; \ 77226606Sdas} while (0) 78226606Sdas 79226606Sdas/* Set a double from two 32 bit ints. */ 80226606Sdas 81226606Sdas#define INSERT_WORDS(d,ix0,ix1) \ 82226606Sdasdo { \ 83226606Sdas ieee_double_shape_type iw_u; \ 84226606Sdas iw_u.parts.msw = (ix0); \ 85226606Sdas iw_u.parts.lsw = (ix1); \ 86226606Sdas (d) = iw_u.value; \ 87226606Sdas} while (0) 88226606Sdas 892116Sjkhstatic const double one = 1.0; 902116Sjkh 9197413Salfreddouble 9297413Salfredmodf(double x, double *iptr) 932116Sjkh{ 942116Sjkh int32_t i0,i1,j0; 952116Sjkh u_int32_t i; 962116Sjkh EXTRACT_WORDS(i0,i1,x); 972116Sjkh j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */ 982116Sjkh if(j0<20) { /* integer part in high x */ 992116Sjkh if(j0<0) { /* |x|<1 */ 1002116Sjkh INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */ 1012116Sjkh return x; 1022116Sjkh } else { 1032116Sjkh i = (0x000fffff)>>j0; 1042116Sjkh if(((i0&i)|i1)==0) { /* x is integral */ 1052116Sjkh u_int32_t high; 1062116Sjkh *iptr = x; 1072116Sjkh GET_HIGH_WORD(high,x); 1082116Sjkh INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ 1092116Sjkh return x; 1102116Sjkh } else { 1112116Sjkh INSERT_WORDS(*iptr,i0&(~i),0); 1122116Sjkh return x - *iptr; 1132116Sjkh } 1142116Sjkh } 1152116Sjkh } else if (j0>51) { /* no fraction part */ 1162116Sjkh u_int32_t high; 117165838Sdas if (j0 == 0x400) { /* inf/NaN */ 118165838Sdas *iptr = x; 119165838Sdas return 0.0 / x; 120165838Sdas } 1212116Sjkh *iptr = x*one; 1222116Sjkh GET_HIGH_WORD(high,x); 1232116Sjkh INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ 1242116Sjkh return x; 1252116Sjkh } else { /* fraction part in low x */ 1262116Sjkh i = ((u_int32_t)(0xffffffff))>>(j0-20); 1272116Sjkh if((i1&i)==0) { /* x is integral */ 1282116Sjkh u_int32_t high; 1292116Sjkh *iptr = x; 1302116Sjkh GET_HIGH_WORD(high,x); 1312116Sjkh INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */ 1322116Sjkh return x; 1332116Sjkh } else { 1342116Sjkh INSERT_WORDS(*iptr,i0,i1&(~i)); 1352116Sjkh return x - *iptr; 1362116Sjkh } 1372116Sjkh } 1382116Sjkh} 139