ieee (empty) → 0.1
raw patch · 6 files changed
+595/−0 lines, 6 filesdep +basebuild-type:Customsetup-changed
Dependencies added: base
Files
- Data/AEq.hs +133/−0
- LICENSE +30/−0
- Numeric/IEEE.hs +147/−0
- Setup.lhs +8/−0
- ieee.cabal +29/−0
- tests/Properties.hs +248/−0
+ Data/AEq.hs view
@@ -0,0 +1,133 @@+-----------------------------------------------------------------------------+-- |+-- Module : Data.AEq+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+-- A type class for approximate and exact equalilty comparisons and instances +-- for common data types. +module Data.AEq+ where++import Data.Int+import Data.Maybe ( fromMaybe )+import Data.Word+import Data.Complex+import Numeric.IEEE++infix 4 ===, ~==++class Eq a => AEq a where+ -- | A reliable way to test if two values are exactly equal. For floating+ -- point values, this will consider @NaN@ to be (===) to @NaN@.+ (===) :: a -> a -> Bool+ + -- | An approximate equality comparison operator. For @RealFloat@ values, + -- @(~==) x y = (x == y)+ -- || (abs (x - y) < epsilon) + -- || (eqRel delta x y) + -- || (isNaN x && isNaN y)@.+ -- For Complex numbers, the if the real and imaginary parts are not+ -- approximately equal, the polar forms are compared, instead.+ (~==) :: a -> a -> Bool+ ++instance AEq Float where+ (===) x y =+ (x == y) || (isNaN x && isNaN y)+ (~==) x y = + (x == y) || (abs (x - y) < epsilon) || (eqRel delta x y) || (isNaN x && isNaN y)++instance AEq Double where+ (===) x y =+ (x == y) || (isNaN x && isNaN y)+ (~==) x y = + (x == y) || (abs (x - y) < epsilon) || (eqRel delta x y) || (isNaN x && isNaN y)++instance (RealFloat a, AEq a) => AEq (Complex a) where+ (===) (x1 :+ y1) (x2 :+ y2) = ((===) x1 x2) && ((===) y1 y2)+ (~==) z1@(x1 :+ y1) z2@(x2 :+ y2) = + let (r1,c1) = polar z1+ (r2,c2) = polar z2+ c = min c1 c2+ c' = max c1 c2+ in (x1 ~== x2 && y1 ~== y2) || (r1 ~== r2) && ((c1 ~== c2) || (c + 2 * pi ~== c'))++instance AEq Bool where+ (===) = (==)+ (~==) = (==)+ +instance AEq Int where+ (===) = (==)+ (~==) = (==)++instance AEq Int8 where+ (===) = (==)+ (~==) = (==)+ +instance AEq Int16 where+ (===) = (==)+ (~==) = (==)+ +instance AEq Int32 where+ (===) = (==)+ (~==) = (==)++instance AEq Int64 where+ (===) = (==)+ (~==) = (==)++instance AEq Word where+ (===) = (==)+ (~==) = (==)++instance AEq Word8 where+ (===) = (==)+ (~==) = (==)+ +instance AEq Word16 where+ (===) = (==)+ (~==) = (==)+ +instance AEq Word32 where+ (===) = (==)+ (~==) = (==)+ +instance AEq Word64 where+ (===) = (==)+ (~==) = (==)++instance AEq () where+ (===) = (==)+ (~==) = (==)+ +instance (AEq a, AEq b) => AEq (a,b) where+ (===) (a1,b1) (a2,b2) = ((===) a1 a2) && ((===) b1 b2)+ (~==) (a1,b1) (a2,b2) = ((~==) a1 a2) && ((~==) b1 b2)++instance (AEq a, AEq b, AEq c) => AEq (a,b,c) where+ (===) (a1,b1,c1) (a2,b2,c2) = ((===) a1 a2) && ((===) b1 b2) && ((===) c1 c2)+ (~==) (a1,b1,c1) (a2,b2,c2) = ((~==) a1 a2) && ((~==) b1 b2) && ((~==) c1 c2)++instance (AEq a, AEq b, AEq c, AEq d) => AEq (a,b,c,d) where+ (===) (a1,b1,c1,d1) (a2,b2,c2,d2) = ((===) a1 a2) && ((===) b1 b2) && ((===) c1 c2) && ((===) d1 d2)+ (~==) (a1,b1,c1,d1) (a2,b2,c2,d2) = ((~==) a1 a2) && ((~==) b1 b2) && ((~==) c1 c2) && ((~==) d1 d2)++instance (AEq a) => AEq [a] where+ (===) xs ys = and $ zipWith (===) xs ys+ (~==) xs ys = and $ zipWith (~==) xs ys++instance (AEq a) => AEq (Maybe a) where+ (===) x y = fromMaybe True $ do x >>= \x' -> y >>= \y' -> return ((===) x' y')+ (~==) x y = fromMaybe True $ do x >>= \x' -> y >>= \y' -> return ((~==) x' y')+ +instance (AEq a, AEq b) => AEq (Either a b) where+ (===) (Left a1) (Left a2) = (===) a1 a2+ (===) (Right b1) (Right b2) = (===) b1 b2+ (===) _ _ = False++ (~==) (Left a1) (Left a2) = (~==) a1 a2+ (~==) (Right b1) (Right b2) = (~==) b1 b2+ (~==) _ _ = False
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) Patrick Perry <patperry@stanford.edu> 2008++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:++1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.++3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Numeric/IEEE.hs view
@@ -0,0 +1,147 @@+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.IEEE+-- Copyright : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>+-- License : BSD3+-- Maintainer : Patrick Perry <patperry@stanford.edu>+-- Stability : experimental+--+-- Approximate comparison of floating point numbers based on the+-- algorithm in Section 4.2.2 of Knuth's _Seminumerical Algorithms_+-- and NaN-aware minimum and maximum.+-- +-- Relative accuracy within @eps@ is measured using an interval of size @2*r@,+-- where @r = 2^k eps@, and @k@ is the maximum exponent of @x@ and @y@. If +-- @x@ and @y@ lie within this interval, they are considered approximately +-- equal.+-- +-- Note that @x@ and @y@ are compared to relative accuracy, so these functions+-- are not suitable for testing whether a value is approximately zero.+--+-- The implementation is based on the GNU Scientific Library implementation, +-- which is based on the package @fcmp@ by T.C. Belding.+module Numeric.IEEE (+ + -- * NaN-aware minimum and maximum+ maxF,+ minF,+ + -- * Relative comparisons+ delta,+ epsilon,+ epsilon',++ eqRel,+ neqRel,+ ltRel,+ lteRel,+ gtRel,+ gteRel,+ compareRel,+ ) where++-- | A version of 'max' that returns @NaN@ if either argument is @NaN@.+maxF :: RealFloat a => a -> a -> a+maxF a b+ | isNaN a = a+ | b < a = a+ | otherwise = b+{-# INLINE maxF #-}++-- | A version of 'min' that returns @NaN@ if either argument is @NaN@.+minF :: RealFloat a => a -> a -> a+minF a b+ | isNaN a = a+ | b > a = a+ | otherwise = b+{-# INLINE minF #-}+++epsHelp :: RealFloat a => (Int -> Int) -> a+epsHelp = epsHelp' undefined+ where+ epsHelp' :: RealFloat a => a -> (Int -> Int) -> a+ epsHelp' a f =+ let digits = floatDigits a+ in encodeFloat 1 $ f digits++-- | A value suitable for relative comparisons when half of of the +-- digits of precision are important. For @Double@s this value is +-- @7.450580596923828e-9@.+delta :: RealFloat a => a+delta = epsHelp (\digits -> negate $ digits `div` 2 + 1)++-- | The smallest positive floating-point number x such that @1 + x != 1@.+-- Suitable for relative comparisons when all but the least significant digit+-- of precision are important. For @Double@s this value is +-- @2.220446049250313e-16@.+epsilon :: RealFloat a => a+epsilon = epsHelp (\digits -> negate $ digits - 1)++-- | The smallest positive floating-point number x such that @1 - x != 1@.+-- Suitable for relative comparisons when one number is exact and all but the +-- least significant digit of precision in the other number are important. +-- For @Double@s this value is @1.1102230246251565e-16@.+epsilon' :: RealFloat a => a+epsilon' = epsHelp (\digits -> negate $ digits)+++compareRelHelp :: (RealFloat a) => (a -> a -> Bool) -> a -> a -> a -> Bool+compareRelHelp cmp tol x y =+ let e = max (exponent x) (exponent y)+ (epsM,epsE) = decodeFloat tol+ r = encodeFloat epsM (epsE + e)+ diff = x - y+ in+ diff `cmp` r+{-# INLINE compareRelHelp #-}++-- | @eqRel eps x y@. Relative equality comparator.+-- Returns @False@ if either argument is @NaN@.+eqRel :: (RealFloat a) => a -> a -> a -> Bool+eqRel = compareRelHelp (\diff r -> abs diff <= r)+{-# INLINE eqRel #-}++-- | @neqRel eps x y@. Relative inequality comparator.+-- Returns @False@ if either argument is @NaN@.+neqRel :: (RealFloat a) => a -> a -> a -> Bool+neqRel = compareRelHelp (\diff r -> abs diff > r)+{-# INLINE neqRel #-}++-- | @ltRel eps x y@. Relative less-than comparator. +-- Returns @False@ if either argument is @NaN@.+ltRel :: (RealFloat a) => a -> a -> a -> Bool+ltRel = compareRelHelp (\diff r -> diff < -r)+{-# INLINE ltRel #-}++-- | @lteRel eps x y@. Relative less-than-or-equal-to comparator. +-- Returns @False@ if either argument is @NaN@.+lteRel :: (RealFloat a) => a -> a -> a -> Bool+lteRel = compareRelHelp (\diff r -> diff <= r)+{-# INLINE lteRel #-}++-- | @gtRel eps x y@. Relative greater-than comparator.+-- Returns @False@ if either argument is @NaN@.+gtRel :: (RealFloat a) => a -> a -> a -> Bool+gtRel = compareRelHelp (\diff r -> diff > r)+{-# INLINE gtRel #-}++-- | @gteRel eps x y@. Relative greater-than-or-equal-to comparator.+-- Returns @False@ if either argument is @NaN@.+gteRel :: (RealFloat a) => a -> a -> a -> Bool+gteRel = compareRelHelp (\diff r -> diff >= -r)+{-# INLINE gteRel #-}++-- | @compareRel eps x y@ gives an ordering of @x@ and @y@ based on a +-- relative comparison of accuracy @eps@. This will call @error@ if either+-- argument is @NaN@.+compareRel :: RealFloat a => a -> a -> a -> Ordering+compareRel eps x y =+ if ltRel eps x y+ then LT+ else if gtRel eps x y+ then GT+ else if eqRel eps x y+ then EQ+ else error $ "NaN comparison"+{-# INLINE compareRel #-}
+ Setup.lhs view
@@ -0,0 +1,8 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> import System.Cmd+>+> testing _ _ _ _ = system "runhaskell tests/Properties.hs" >> return ()+>+> main = defaultMainWithHooks defaultUserHooks+> {runTests=testing}
+ ieee.cabal view
@@ -0,0 +1,29 @@+name: ieee+version: 0.1+homepage: http://stat.stanford.edu/~patperry/code/ieee+synopsis: Approximate comparisons for IEEE floating point numbers+description:+ Approximate comparison of floating point numbers based on the+ algorithm in Section 4.2.2 of Knuth's _Seminumerical Algorithms_,+ NaN-aware minimum and maximum, and a type class for approximate + comparisons.+ .+category: Math+license: BSD3+license-file: LICENSE+copyright: (c) 2008. Patrick Perry <patperry@stanford.edu>+author: Patrick Perry+maintainer: Patrick Perry <patperry@stanford.edu>+cabal-version: >= 1.2.0+build-type: Custom+tested-with: GHC ==6.8.2++extra-source-files: tests/Properties.hs++library+ exposed-modules: Data.AEq+ Numeric.IEEE++ build-depends: base++ ghc-options: -Wall
+ tests/Properties.hs view
@@ -0,0 +1,248 @@+{-# OPTIONS -fglasgow-exts #-}++import Test.QuickCheck+import Data.Int++import Numeric.IEEE++type D = Double++b = floatRadix (undefined :: D)++nan :: RealFloat a => a+nan = 0 / 0++incSignif :: RealFloat a => Integer -> a -> a+incSignif i x =+ let (m,n) = decodeFloat x+ in encodeFloat (m+i) n+++------------------------- NaN-aware min and max -----------------------------++prop_minF (x :: D) (y :: D) =+ not (isNaN x || isNaN y) ==> minF x y == min x y++prop_minF_nan1 (x :: D) =+ isNaN (minF x nan)+ +prop_minF_nan2 (x :: D) =+ isNaN (minF nan x)+++prop_maxF (x :: D) (y :: D) =+ not (isNaN x || isNaN y) ==> maxF x y == max x y++prop_maxF_nan1 (x :: D) =+ isNaN (maxF x nan)+ +prop_maxF_nan2 (x :: D) =+ isNaN (maxF nan x)+ +----------------------------- NaN comparisons -------------------------------++prop_nan1 (x :: D) (eps :: D) =+ eps >= 0 ==> (not . or) [ eqRel eps x nan+ , neqRel eps x nan+ , ltRel eps x nan+ , lteRel eps x nan+ , gtRel eps x nan+ , gteRel eps x nan+ ]++prop_nan2 (x :: D) (eps :: D) =+ eps >= 0 ==> (not . or) [ eqRel eps nan x+ , neqRel eps nan x+ , ltRel eps nan x+ , lteRel eps nan x+ , gtRel eps nan x+ , gteRel eps nan x+ ]+++--------------------- Comparisons relative to 0 -----------------------------++prop_0_exact (x :: D) =+ not (isNaN x || x == 0) ==> and [ eqRel 0 x x'+ , lteRel 0 x x'+ , gteRel 0 x x'+ ]+ && (not . or) [ neqRel 0 x x'+ , ltRel 0 x x'+ , gtRel 0 x x'+ ]+ && (compareRel epsilon' x x' == EQ)+ where + x' = x++prop_plus_1_exact (x :: D) =+ not (isNaN x || x == 0) ==> and [ neqRel 0 x x'+ , ltRel 0 x x'+ , lteRel 0 x x'+ ]+ && (not . or) [ eqRel 0 x x'+ , gtRel 0 x x'+ , gteRel 0 x x'+ ]+ && (compareRel 0 x x' == LT)+ where + x' = incSignif 1 x++prop_minus_1_exact (x :: D) =+ not (isNaN x || x == 0) ==> and [ neqRel 0 x x'+ , gtRel 0 x x'+ , gteRel 0 x x'+ ]+ && (not . or) [ eqRel 0 x x'+ , ltRel 0 x x'+ , lteRel 0 x x'+ ]+ && (compareRel 0 x x' == GT)+ where + x' = incSignif (-1) x++--------------------- Comparisons relative to epsilon' ----------------------++prop_plus_1' (x :: D) =+ not (isNaN x || x == 0) ==> and [ eqRel epsilon' x x'+ , lteRel epsilon' x x'+ , gteRel epsilon' x x'+ ]+ && (not . or) [ neqRel epsilon' x x'+ , ltRel epsilon' x x'+ , gtRel epsilon' x x'+ ]+ && (compareRel epsilon' x x' == EQ)+ where + x' = incSignif 1 x++prop_minus_1' (x :: D) =+ not (isNaN x || x == 0) ==> and [ eqRel epsilon' x x'+ , lteRel epsilon' x x'+ , gteRel epsilon' x x'+ ]+ && (not . or) [ neqRel epsilon' x x'+ , ltRel epsilon' x x'+ , gtRel epsilon' x x'+ ]+ && (compareRel epsilon' x x' == EQ)+ where + x' = incSignif (-1) x+++prop_plus_b' (x :: D) =+ not (isNaN x || x == 0) ==> and [ neqRel epsilon' x x'+ , ltRel epsilon' x x'+ , lteRel epsilon' x x'+ ]+ && (not . or) [ eqRel epsilon' x x'+ , gtRel epsilon' x x'+ , gteRel epsilon' x x'+ ]+ && (compareRel epsilon' x x' == LT)+ where + x' = incSignif b x++prop_minus_b' (x :: D) =+ not (isNaN x || x == 0) ==> and [ neqRel epsilon' x x'+ , gtRel epsilon' x x'+ , gteRel epsilon' x x'+ ]+ && (not . or) [ eqRel epsilon' x x'+ , ltRel epsilon' x x'+ , lteRel epsilon' x x'+ ]+ && (compareRel epsilon' x x' == GT)+ where + x' = incSignif (-b) x+++--------------------- Comparisons relative to epsilon -----------------------++prop_plus_b (x :: D) =+ not (isNaN x || x == 0) ==> and [ eqRel epsilon x x'+ , lteRel epsilon x x'+ , gteRel epsilon x x'+ ]+ && (not . or) [ neqRel epsilon x x'+ , ltRel epsilon x x'+ , gtRel epsilon x x'+ ]+ && (compareRel epsilon x x' == EQ)+ where + x' = incSignif b x++prop_minus_b (x :: D) =+ not (isNaN x || x == 0) ==> and [ eqRel epsilon x x'+ , lteRel epsilon x x'+ , gteRel epsilon x x'+ ]+ && (not . or) [ neqRel epsilon x x'+ , ltRel epsilon x x'+ , gtRel epsilon x x'+ ]+ && (compareRel epsilon x x' == EQ)+ where + x' = incSignif (-b) x+++prop_plus_b1 (x :: D) =+ not (isNaN x || x == 0) ==> and [ neqRel epsilon x x'+ , ltRel epsilon x x'+ , lteRel epsilon x x'+ ]+ && (not . or) [ eqRel epsilon x x'+ , gtRel epsilon x x'+ , gteRel epsilon x x'+ ]+ && (compareRel epsilon x x' == LT)+ where + x' = incSignif (b+1) x++prop_minus_b1 (x :: D) =+ not (isNaN x || x == 0) ==> and [ neqRel epsilon x x'+ , gtRel epsilon x x'+ , gteRel epsilon x x'+ ]+ && (not . or) [ eqRel epsilon x x'+ , ltRel epsilon x x'+ , lteRel epsilon x x'+ ]+ && (compareRel epsilon x x' == GT)+ where + x' = incSignif (-b-1) x++main = do+ let runT s a = do print s; a++ runT "minF" $ do+ quickCheck prop_minF + quickCheck prop_minF_nan1+ quickCheck prop_minF_nan2+ + runT "maxF" $ do+ quickCheck prop_maxF+ quickCheck prop_maxF_nan1 + quickCheck prop_maxF_nan2++ runT "NaN comparisons" $ do+ quickCheck prop_nan1+ quickCheck prop_nan2+ + runT "comparisons relative to 0" $ do+ quickCheck prop_0_exact+ quickCheck prop_plus_1_exact+ quickCheck prop_minus_1_exact+ + runT "comparisons relative to epsilon'" $ do+ quickCheck prop_plus_1'+ quickCheck prop_minus_1'+ quickCheck prop_plus_b'+ quickCheck prop_minus_b'++ runT "comparisons relative to epsilon" $ do+ quickCheck prop_plus_b+ quickCheck prop_minus_b+ quickCheck prop_plus_b1+ quickCheck prop_minus_b1+