```-- |
-- Module    : Numeric.ApproxEq
-- Copyright : (c) 2011 Aleksey Khudyakov, Bryan O'Sullivan
--
-- Maintainer  : Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- 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 Data.Complex             (Complex(..))
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\":
-- <http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm>
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

```