# numeric-tools / Numeric / ApproxEq.hs

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81``` ```-- | -- Module : Numeric.ApproxEq -- Copyright : (c) 2011 Aleksey Khudyakov, Bryan O'Sullivan -- License : BSD3 -- -- Maintainer : Aleksey Khudyakov -- Stability : experimental -- Portability : portable -- -- Different implementations of approximate equality for floating -- point values. There are multiple ways to implement approximate -- equality. They have different semantics and it's up to programmer -- to choose right one. module Numeric.ApproxEq ( eqRelative , eqRelCompl , eqAbsolute , within ) where import Control.Monad.ST (runST) import Data.Complex (Complex(..)) import Data.Primitive.ByteArray (newByteArray, readByteArray, writeByteArray) import Data.Int (Int64) -- | Relative difference between two numbers are less than predefined -- value. For example 1 is approximately equal to 1.0001 with 1e-4 -- precision. Same is true for 10000 and 10001. -- -- This method of camparison doesn't work for numbers which are -- approximately 0. 'eqAbsolute' should be used instead. eqRelative :: (Fractional a, Ord a) => a -- ^ Relative precision -> a -> a -> Bool eqRelative eps a b = abs (a - b) <= abs (eps * a) {-# INLINE eqRelative #-} -- | Relative equality for comlex numbers. eqRelCompl :: (RealFloat a, Ord a) => a -- ^ Relative precision -> Complex a -> Complex a -> Bool eqRelCompl eps a b = d <= eps * r where (d :+ _) = abs (a - b) (r :+ _) = abs a {-# INLINE eqRelCompl #-} -- | Difference between values is less than specified precision. eqAbsolute :: (Num a, Ord a) => a -- ^ Absolute precision -> a -> a -> Bool eqAbsolute eps a b = abs (a - b) <= eps {-# INLINE eqAbsolute #-} -- | Compare two 'Double' values for approximate equality, using -- Dawson's method. -- -- The required accuracy is specified in ULPs (units of least -- precision). If the two numbers differ by the given number of ULPs -- or less, this function returns @True@. -- -- Algorithm is based on Bruce Dawson's \"Comparing floating point -- numbers\": -- within :: Int -- ^ Number of ULPs of accuracy desired. -> Double -> Double -> Bool within ulps a b = runST \$ do buf <- newByteArray 8 ai0 <- writeByteArray buf 0 a >> readByteArray buf 0 bi0 <- writeByteArray buf 0 b >> readByteArray buf 0 let big = 0x8000000000000000 :: Int64 ai | ai0 < 0 = big - ai0 | otherwise = ai0 bi | bi0 < 0 = big - bi0 | otherwise = bi0 return \$! abs (ai - bi) <= fromIntegral ulps ```