diff --git a/Data/AEq.hs b/Data/AEq.hs
new file mode 100644
--- /dev/null
+++ b/Data/AEq.hs
@@ -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
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Numeric/IEEE.hs b/Numeric/IEEE.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/IEEE.hs
@@ -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 #-}
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -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}
diff --git a/ieee.cabal b/ieee.cabal
new file mode 100644
--- /dev/null
+++ b/ieee.cabal
@@ -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
diff --git a/tests/Properties.hs b/tests/Properties.hs
new file mode 100644
--- /dev/null
+++ b/tests/Properties.hs
@@ -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
+