ieee754 0.7.1 → 0.7.2
raw patch · 5 files changed
+111/−101 lines, 5 filesPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
API changes (from Hackage documentation)
- Data.AEq: instance (AEq a) => AEq (Maybe a)
- Data.AEq: instance (AEq a) => AEq [a]
+ Data.AEq: instance AEq a => AEq (Maybe a)
+ Data.AEq: instance AEq a => AEq [a]
- Data.AEq: (===) :: (AEq a) => a -> a -> Bool
+ Data.AEq: (===) :: AEq a => a -> a -> Bool
- Data.AEq: (~==) :: (AEq a) => a -> a -> Bool
+ Data.AEq: (~==) :: AEq a => a -> a -> Bool
- Data.AEq: class (Eq a) => AEq a
+ Data.AEq: class Eq a => AEq a
- Numeric.IEEE: bisectIEEE :: (IEEE a) => a -> a -> a
+ Numeric.IEEE: bisectIEEE :: IEEE a => a -> a -> a
- Numeric.IEEE: class (RealFloat a) => IEEE a
+ Numeric.IEEE: class RealFloat a => IEEE a
- Numeric.IEEE: copySign :: (IEEE a) => a -> a -> a
+ Numeric.IEEE: copySign :: IEEE a => a -> a -> a
- Numeric.IEEE: epsilon :: (IEEE a) => a
+ Numeric.IEEE: epsilon :: IEEE a => a
- Numeric.IEEE: identicalIEEE :: (IEEE a) => a -> a -> Bool
+ Numeric.IEEE: identicalIEEE :: IEEE a => a -> a -> Bool
- Numeric.IEEE: infinity :: (IEEE a) => a
+ Numeric.IEEE: infinity :: IEEE a => a
- Numeric.IEEE: maxFinite :: (IEEE a) => a
+ Numeric.IEEE: maxFinite :: IEEE a => a
- Numeric.IEEE: maxNaN :: (RealFloat a) => a -> a -> a
+ Numeric.IEEE: maxNaN :: RealFloat a => a -> a -> a
- Numeric.IEEE: maxNaNPayload :: (IEEE a) => a -> Word64
+ Numeric.IEEE: maxNaNPayload :: IEEE a => a -> Word64
- Numeric.IEEE: maxNum :: (RealFloat a) => a -> a -> a
+ Numeric.IEEE: maxNum :: RealFloat a => a -> a -> a
- Numeric.IEEE: minNaN :: (RealFloat a) => a -> a -> a
+ Numeric.IEEE: minNaN :: RealFloat a => a -> a -> a
- Numeric.IEEE: minNormal :: (IEEE a) => a
+ Numeric.IEEE: minNormal :: IEEE a => a
- Numeric.IEEE: minNum :: (RealFloat a) => a -> a -> a
+ Numeric.IEEE: minNum :: RealFloat a => a -> a -> a
- Numeric.IEEE: nan :: (IEEE a) => a
+ Numeric.IEEE: nan :: IEEE a => a
- Numeric.IEEE: nanPayload :: (IEEE a) => a -> Word64
+ Numeric.IEEE: nanPayload :: IEEE a => a -> Word64
- Numeric.IEEE: nanWithPayload :: (IEEE a) => Word64 -> a
+ Numeric.IEEE: nanWithPayload :: IEEE a => Word64 -> a
- Numeric.IEEE: predIEEE :: (IEEE a) => a -> a
+ Numeric.IEEE: predIEEE :: IEEE a => a -> a
- Numeric.IEEE: sameSignificandBits :: (IEEE a) => a -> a -> Int
+ Numeric.IEEE: sameSignificandBits :: IEEE a => a -> a -> Int
- Numeric.IEEE: succIEEE :: (IEEE a) => a -> a
+ Numeric.IEEE: succIEEE :: IEEE a => a -> a
Files
- NEWS +5/−0
- cbits/double.c +33/−31
- cbits/feqrel_source.c +40/−40
- cbits/float.c +32/−29
- ieee754.cabal +1/−1
NEWS view
@@ -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"
cbits/double.c view
@@ -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; }
cbits/feqrel_source.c view
@@ -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
cbits/float.c view
@@ -3,6 +3,7 @@ #include <math.h> #include <stdint.h> + #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; }
ieee754.cabal view
@@ -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: