{-# LANGUAGE GADTs #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ScopedTypeVariables #-}
-- | See: https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
module Data.Double.Approximate
( SafeDouble
, DoubleRelAbs(..)
) where
import Control.DeepSeq
import Data.Proxy
import GHC.TypeLits
import Numeric.MathFunctions.Comparison
import Numeric.MathFunctions.Constants
import System.Random
import Text.Read
-- | Relatively safe double floating-point type with a relative error
-- margin of 10 <https://en.wikipedia.org/wiki/Unit_in_the_last_place ULPs>
-- and an absolute margin around zero of
-- @10*<https://en.wikipedia.org/wiki/Machine_epsilon epsilon>@.
--
-- Warning: All numbers within
-- @10*<https://en.wikipedia.org/wiki/Machine_epsilon epsilon>@ of zero will be considered zero.
--
-- >>> m_epsilon * 10
-- 2.220446049250313e-15
--
-- >>> realToFrac (m_epsilon * 10) == (0::SafeDouble)
-- False
--
-- >>> realToFrac (m_epsilon * 9) == (0::SafeDouble)
-- True
--
-- >>> 1e-20 == (5e-20 :: Double)
-- False
-- >>> 1e-20 == (5e-20 :: SafeDouble)
-- True
--
-- 'pi' and 'sin' are approximations:
--
-- >>> sin pi
-- 1.2246467991473532e-16
--
-- >>> sin pi == (0 :: Double)
-- False
--
-- >>> sin pi == (0 :: SafeDouble)
-- True
--
type SafeDouble = DoubleRelAbs 10 10
-- | Custom double floating-point type with a relative error margin of
-- @rel@ number of
-- <https://en.wikipedia.org/wiki/Unit_in_the_last_place ULPs> and an
-- absolute error margin of @abs@ times
-- <https://en.wikipedia.org/wiki/Machine_epsilon epsilon>.
--
-- The relative error margin is the primary tool for good numerical
-- robustness and can relatively safely be set to a high number such
-- as 100. The absolute error margin is a last ditch attempt at fixing
-- broken algorithms and dramatically limits the resolution around zero.
-- If possible, use a low absolute error margin.
newtype DoubleRelAbs (abs :: Nat) (rel :: Nat) = DoubleRelAbs Double
deriving (Num, Enum, Floating, Fractional, Real, RealFloat, RealFrac, Random, NFData)
instance (KnownNat abs, KnownNat rel) => Eq (DoubleRelAbs abs rel) where
DoubleRelAbs d1 == DoubleRelAbs d2 =
within (fromIntegral (natVal @rel Proxy)) d1 d2 ||
(abs d1 < m_epsilon * fromIntegral (natVal @abs Proxy) &&
abs d2 < m_epsilon * fromIntegral (natVal @abs Proxy))
instance (KnownNat abs, KnownNat rel) => Ord (DoubleRelAbs abs rel) where
lhs@(DoubleRelAbs d1) `compare` rhs@(DoubleRelAbs d2)
| lhs == rhs = EQ
| otherwise = d1 `compare` d2
instance Show (DoubleRelAbs abs rel) where
showsPrec i (DoubleRelAbs d) = showsPrec i d
instance Read (DoubleRelAbs abs rel) where
readPrec = DoubleRelAbs <$> readPrec