diff --git a/NEWS b/NEWS --- a/NEWS +++ b/NEWS @@ -1,3 +1,8 @@ +Changes in 0.7.2: + +* Fix aliasing warnings in C code + + Changes in 0.7.1: * Rename package from "ieee" to "ieee754" diff --git a/cbits/double.c b/cbits/double.c --- a/cbits/double.c +++ b/cbits/double.c @@ -14,35 +14,40 @@ #define FEQREL feqrel #include "feqrel_source.c" +union double_t { + double d; + uint64_t w; +}; + int identical (double x, double y) { - uint64_t *ux = (uint64_t *)(&x); - uint64_t *uy = (uint64_t *)(&y); - return *ux == *uy; + union double_t ux = { x }; + union double_t uy = { y }; + return ux.w == uy.w; } /* ported from tango/math/IEEE.d */ double nextup (double x) { - uint64_t *ps = (uint64_t *)&x; + union double_t ps = { x }; - if ((*ps & 0x7FF0000000000000ULL) == 0x7FF0000000000000ULL) { + if ((ps.w & 0x7FF0000000000000ULL) == 0x7FF0000000000000ULL) { /* First, deal with NANs and infinity */ if (x == -INFINITY) return -REAL_MAX; return x; // +INF and NAN are unchanged. } - if (*ps & 0x8000000000000000ULL) { /* Negative number */ - if (*ps == 0x8000000000000000ULL) { /* it was negative zero */ - *ps = 0x0000000000000001ULL; /* change to smallest subnormal */ - return x; + if (ps.w & 0x8000000000000000ULL) { /* Negative number */ + if (ps.w == 0x8000000000000000ULL) { /* it was negative zero */ + ps.w = 0x0000000000000001ULL; /* change to smallest subnormal */ + return ps.d; } - --*ps; + --ps.w; } else { /* Positive number */ - ++*ps; + ++ps.w; } - return x; + return ps.d; } /* ported from tango/math/IEEE.d */ @@ -58,41 +63,38 @@ { if (!((x>=0 && y>=0) || (x<=0 && y<=0))) return NAN; - double u; - - uint64_t *ul = (uint64_t *)&u; - uint64_t *xl = (uint64_t *)&x; - uint64_t *yl = (uint64_t *)&y; - uint64_t m = ( ((*xl) & 0x7FFFFFFFFFFFFFFFULL) - + ((*yl) & 0x7FFFFFFFFFFFFFFFULL) ) >> 1; - m |= ((*xl) & 0x8000000000000000ULL); - *ul = m; + union double_t ul; + union double_t xl = { x }; + union double_t yl = { y }; + uint64_t m = ( (xl.w & 0x7FFFFFFFFFFFFFFFULL) + + (yl.w & 0x7FFFFFFFFFFFFFFFULL) ) >> 1; + m |= (xl.w & 0x8000000000000000ULL); + ul.w = m; - return u; + return ul.d; } double mknan (uint64_t payload) { - double x = NAN; - uint64_t *ux = (uint64_t *)(&x); + union double_t x = { NAN }; /* get sign, exponent, and quiet bit from NAN */ - *ux &= 0xFFF8000000000000ULL; + x.w &= 0xFFF8000000000000ULL; /* ignore sign, exponent, and quiet bit in payload */ payload &= 0x0007FFFFFFFFFFFFULL; - *ux |= payload; + x.w |= payload; - return x; + return x.d; } uint64_t getnan (double x) { - uint64_t payload = *(uint64_t *)(&x); + union double_t payload = { x }; - /* clear ignore sign, exponent, and quiet bit */ - payload &= 0x0007FFFFFFFFFFFFULL; - return payload; + /* clear sign, exponent, and quiet bit */ + payload.w &= 0x0007FFFFFFFFFFFFULL; + return payload.w; } diff --git a/cbits/feqrel_source.c b/cbits/feqrel_source.c --- a/cbits/feqrel_source.c +++ b/cbits/feqrel_source.c @@ -40,68 +40,68 @@ int FEQREL (REAL x, REAL y) { - /* Public Domain. Original Author: Don Clugston, 18 Aug 2005. - * Ported to C by Patrick Perry, 26 Feb 2010. - */ - if (x == y) return REAL_MANT_DIG; /* ensure diff!= 0, cope with INF. */ + /* Public Domain. Original Author: Don Clugston, 18 Aug 2005. + * Ported to C by Patrick Perry, 26 Feb 2010. + */ + if (x == y) return REAL_MANT_DIG; /* ensure diff!= 0, cope with INF. */ - REAL diff = REAL_ABS(x - y); + REAL diff = REAL_ABS(x - y); - uint16_t *pa = (uint16_t *)(&x); - uint16_t *pb = (uint16_t *)(&y); - uint16_t *pd = (uint16_t *)(&diff); + union { REAL r; uint16_t w[sizeof(REAL)/2]; } pa = { x }; + union { REAL r; uint16_t w[sizeof(REAL)/2]; } pb = { y }; + union { REAL r; uint16_t w[sizeof(REAL)/2]; } pd = { diff }; - /* The difference in abs(exponent) between x or y and abs(x-y) - * is equal to the number of significand bits of x which are - * equal to y. If negative, x and y have different exponents. - * If positive, x and y are equal to 'bitsdiff' bits. - * AND with 0x7FFF to form the absolute value. - * To avoid out-by-1 errors, we subtract 1 so it rounds down - * if the exponents were different. This means 'bitsdiff' is - * always 1 lower than we want, except that if bitsdiff==0, - * they could have 0 or 1 bits in common. - */ + /* The difference in abs(exponent) between x or y and abs(x-y) + * is equal to the number of significand bits of x which are + * equal to y. If negative, x and y have different exponents. + * If positive, x and y are equal to 'bitsdiff' bits. + * AND with 0x7FFF to form the absolute value. + * To avoid out-by-1 errors, we subtract 1 so it rounds down + * if the exponents were different. This means 'bitsdiff' is + * always 1 lower than we want, except that if bitsdiff==0, + * they could have 0 or 1 bits in common. + */ #if REAL_MANT_DIG == 53 /* double */ - int bitsdiff = (( ((pa[REAL_EXPPOS_INT16] & REAL_EXPMASK) - + (pb[REAL_EXPPOS_INT16] & REAL_EXPMASK) + int bitsdiff = (( ((pa.w[REAL_EXPPOS_INT16] & REAL_EXPMASK) + + (pb.w[REAL_EXPPOS_INT16] & REAL_EXPMASK) - ((uint16_t) 0x8000 - REAL_EXPMASK)) >> 1) - - (pd[REAL_EXPPOS_INT16] & REAL_EXPMASK)) >> 4; + - (pd.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)) >> 4; #elif REAL_MANT_DIG == 24 /* float */ - int bitsdiff = (( ((pa[REAL_EXPPOS_INT16] & REAL_EXPMASK) - + (pb[REAL_EXPPOS_INT16] & REAL_EXPMASK) + int bitsdiff = (( ((pa.w[REAL_EXPPOS_INT16] & REAL_EXPMASK) + + (pb.w[REAL_EXPPOS_INT16] & REAL_EXPMASK) - ((uint16_t) 0x8000 - REAL_EXPMASK)) >> 1) - - (pd[REAL_EXPPOS_INT16] & REAL_EXPMASK)) >> 7; + - (pd.w[REAL_EXPPOS_INT16] & REAL_EXPMASK)) >> 7; #else # error unsuported floating-point mantissa size #endif - if ((pd[REAL_EXPPOS_INT16] & REAL_EXPMASK) == 0) { - /* Difference is denormal - * For denormals, we need to add the number of zeros that - * lie at the start of diff's significand. - * We do this by multiplying by 2^REAL_MANT_DIG - */ - diff *= REAL_RECIP_EPSILON; + if ((pd.w[REAL_EXPPOS_INT16] & REAL_EXPMASK) == 0) { + /* Difference is denormal + * For denormals, we need to add the number of zeros that + * lie at the start of diff's significand. + * We do this by multiplying by 2^REAL_MANT_DIG + */ + pd.r *= REAL_RECIP_EPSILON; #if REAL_MANT_DIG == 53 /* double */ return (bitsdiff + REAL_MANT_DIG - - (pd[REAL_EXPPOS_INT16] >> 4)); + - (pd.w[REAL_EXPPOS_INT16] >> 4)); #elif REAL_MANT_DIG == 24 /* float */ return (bitsdiff + REAL_MANT_DIG - - (pd[REAL_EXPPOS_INT16] >> 7)); + - (pd.w[REAL_EXPPOS_INT16] >> 7)); #else # error unsuported floating-point mantissa size #endif - } + } - if (bitsdiff > 0) - return bitsdiff + 1; /* add the 1 we subtracted before */ + if (bitsdiff > 0) + return bitsdiff + 1; /* add the 1 we subtracted before */ - /* Avoid out-by-1 errors when factor is almost 2. */ - return (bitsdiff == 0 - && !((pa[REAL_EXPPOS_INT16] - ^ pb[REAL_EXPPOS_INT16]) & REAL_EXPMASK)) ? 1 : 0; + /* Avoid out-by-1 errors when factor is almost 2. */ + return (bitsdiff == 0 + && !((pa.w[REAL_EXPPOS_INT16] + ^ pb.w[REAL_EXPPOS_INT16]) & REAL_EXPMASK)) ? 1 : 0; } #undef REAL_RECIP_EPSILON diff --git a/cbits/float.c b/cbits/float.c --- a/cbits/float.c +++ b/cbits/float.c @@ -3,6 +3,7 @@ #include #include + #define REAL float #define REAL_ABS fabsf @@ -14,35 +15,40 @@ #define FEQREL feqrelf #include "feqrel_source.c" +union float_t { + float f; + uint32_t w; +}; + int identicalf (float x, float y) { - uint32_t *ux = (uint32_t *)(&x); - uint32_t *uy = (uint32_t *)(&y); - return *ux == *uy; + union float_t ux = { x }; + union float_t uy = { y }; + return ux.w == uy.w; } /* ported from tango/math/IEEE.d */ float nextupf (float x) { - uint32_t *ps = (uint32_t *)&x; + union float_t ps = { x }; - if ((*ps & 0x7F800000) == 0x7F800000) { + if ((ps.w & 0x7F800000) == 0x7F800000) { /* First, deal with NANs and infinity */ if (x == -INFINITY) return -REAL_MAX; return x; /* +INF and NAN are unchanged. */ } - if (*ps & 0x80000000) { /* Negative number */ - if (*ps == 0x80000000) { /* it was negative zero */ - *ps = 0x00000001; /* change to smallest subnormal */ - return x; + if (ps.w & 0x80000000) { /* Negative number */ + if (ps.w == 0x80000000) { /* it was negative zero */ + ps.w = 0x00000001; /* change to smallest subnormal */ + return ps.f; } - --*ps; + --ps.w; } else { /* Positive number */ - ++*ps; + ++ps.w; } - return x; + return ps.f; } /* ported from tango/math/IEEE.d */ @@ -58,40 +64,37 @@ { if (!((x>=0 && y>=0) || (x<=0 && y<=0))) return NAN; - float u; - - uint32_t *ul = (uint32_t *)&u; - uint32_t *xl = (uint32_t *)&x; - uint32_t *yl = (uint32_t *)&y; - uint32_t m = (((*xl) & 0x7FFFFFFF) + ((*yl) & 0x7FFFFFFF)) >> 1; - m |= ((*xl) & 0x80000000); - *ul = m; + union float_t ul; + union float_t xl = { x }; + union float_t yl = { y }; + uint32_t m = ((xl.w & 0x7FFFFFFF) + (yl.w & 0x7FFFFFFF)) >> 1; + m |= (xl.w & 0x80000000); + ul.w = m; - return u; + return ul.f; } float mknanf (uint32_t payload) { - float x = NAN; - uint32_t *ux = (uint32_t *)(&x); + union float_t ux = { NAN }; /* get sign, exponent, and quiet bit from NAN */ - *ux &= 0xFFC00000; + ux.w &= 0xFFC00000; /* ignore sign, exponent, and quiet bit in payload */ payload &= 0x003FFFFF; - *ux |= payload; + ux.w |= payload; - return x; + return ux.f; } uint32_t getnanf (float x) { - uint32_t payload = *(uint32_t *)(&x); + union float_t payload = { x }; /* clear sign, exponent, and quiet bit */ - payload &= 0x003FFFFF; - return payload; + payload.w &= 0x003FFFFF; + return payload.w; } diff --git a/ieee754.cabal b/ieee754.cabal --- a/ieee754.cabal +++ b/ieee754.cabal @@ -1,5 +1,5 @@ name: ieee754 -version: 0.7.1 +version: 0.7.2 homepage: http://github.com/patperry/hs-ieee754 synopsis: Utilities for dealing with IEEE floating point numbers description: