module Numeric.MathFunctions.Comparison
(
relativeError
, eqRelErr
, addUlps
, ulpDistance
, ulpDelta
, within
) where
import Control.Monad.ST (runST)
import Data.Primitive.ByteArray (newByteArray, readByteArray, writeByteArray)
import Data.Word (Word64)
import Data.Int (Int64)
relativeError :: Double -> Double -> Double
relativeError :: Double -> Double -> Double
relativeError Double
a Double
b
| Double
a forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
&& Double
b forall a. Eq a => a -> a -> Bool
== Double
0 = Double
0
| Bool
otherwise = forall a. Num a => a -> a
abs (Double
a forall a. Num a => a -> a -> a
- Double
b) forall a. Fractional a => a -> a -> a
/ forall a. Ord a => a -> a -> a
max (forall a. Num a => a -> a
abs Double
a) (forall a. Num a => a -> a
abs Double
b)
eqRelErr :: Double
-> Double
-> Double
-> Bool
eqRelErr :: Double -> Double -> Double -> Bool
eqRelErr Double
eps Double
a Double
b = Double -> Double -> Double
relativeError Double
a Double
b forall a. Ord a => a -> a -> Bool
< Double
eps
addUlps :: Int -> Double -> Double
addUlps :: Int -> Double -> Double
addUlps Int
n Double
a = forall a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
MutableByteArray s
buf <- forall (m :: * -> *).
PrimMonad m =>
Int -> m (MutableByteArray (PrimState m))
newByteArray Int
8
Word64
ai0 <- forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> a -> m ()
writeByteArray MutableByteArray s
buf Int
0 Double
a forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> m a
readByteArray MutableByteArray s
buf Int
0
let big :: Word64
big = Word64
0x8000000000000000
order :: Word64 -> Int64
order :: Word64 -> Int64
order Word64
i | Word64
i forall a. Ord a => a -> a -> Bool
< Word64
big = forall a b. (Integral a, Num b) => a -> b
fromIntegral Word64
i
| Bool
otherwise = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a. Bounded a => a
maxBound forall a. Num a => a -> a -> a
- (Word64
i forall a. Num a => a -> a -> a
- Word64
big)
unorder :: Int64 -> Word64
unorder :: Int64 -> Word64
unorder Int64
i | Int64
i forall a. Ord a => a -> a -> Bool
>= Int64
0 = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int64
i
| Bool
otherwise = Word64
big forall a. Num a => a -> a -> a
+ (forall a. Bounded a => a
maxBound forall a. Num a => a -> a -> a
- (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int64
i))
let ai0' :: Word64
ai0' = Int64 -> Word64
unorder forall a b. (a -> b) -> a -> b
$ Word64 -> Int64
order Word64
ai0 forall a. Num a => a -> a -> a
+ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> a -> m ()
writeByteArray MutableByteArray s
buf Int
0 Word64
ai0' forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> m a
readByteArray MutableByteArray s
buf Int
0
ulpDistance :: Double
-> Double
-> Word64
ulpDistance :: Double -> Double -> Word64
ulpDistance Double
a Double
b = forall a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
MutableByteArray s
buf <- forall (m :: * -> *).
PrimMonad m =>
Int -> m (MutableByteArray (PrimState m))
newByteArray Int
8
Word64
ai0 <- forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> a -> m ()
writeByteArray MutableByteArray s
buf Int
0 Double
a forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> m a
readByteArray MutableByteArray s
buf Int
0
Word64
bi0 <- forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> a -> m ()
writeByteArray MutableByteArray s
buf Int
0 Double
b forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> m a
readByteArray MutableByteArray s
buf Int
0
let big :: Word64
big = Word64
0x8000000000000000
order :: Word64 -> Word64
order Word64
i | Word64
i forall a. Ord a => a -> a -> Bool
< Word64
big = Word64
i forall a. Num a => a -> a -> a
+ Word64
big
| Bool
otherwise = forall a. Bounded a => a
maxBound forall a. Num a => a -> a -> a
- Word64
i
ai :: Word64
ai = Word64 -> Word64
order Word64
ai0
bi :: Word64
bi = Word64 -> Word64
order Word64
bi0
d :: Word64
d | Word64
ai forall a. Ord a => a -> a -> Bool
> Word64
bi = Word64
ai forall a. Num a => a -> a -> a
- Word64
bi
| Bool
otherwise = Word64
bi forall a. Num a => a -> a -> a
- Word64
ai
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$! Word64
d
ulpDelta :: Double
-> Double
-> Int64
ulpDelta :: Double -> Double -> Int64
ulpDelta Double
a Double
b = forall a. (forall s. ST s a) -> a
runST forall a b. (a -> b) -> a -> b
$ do
MutableByteArray s
buf <- forall (m :: * -> *).
PrimMonad m =>
Int -> m (MutableByteArray (PrimState m))
newByteArray Int
8
Word64
ai0 <- forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> a -> m ()
writeByteArray MutableByteArray s
buf Int
0 Double
a forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> m a
readByteArray MutableByteArray s
buf Int
0
Word64
bi0 <- forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> a -> m ()
writeByteArray MutableByteArray s
buf Int
0 Double
b forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> forall a (m :: * -> *).
(Prim a, PrimMonad m) =>
MutableByteArray (PrimState m) -> Int -> m a
readByteArray MutableByteArray s
buf Int
0
let big :: Word64
big = Word64
0x8000000000000000 :: Word64
order :: Word64 -> Word64
order Word64
i | Word64
i forall a. Ord a => a -> a -> Bool
< Word64
big = Word64
i forall a. Num a => a -> a -> a
+ Word64
big
| Bool
otherwise = forall a. Bounded a => a
maxBound forall a. Num a => a -> a -> a
- Word64
i
ai :: Word64
ai = Word64 -> Word64
order Word64
ai0
bi :: Word64
bi = Word64 -> Word64
order Word64
bi0
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$! forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ Word64
bi forall a. Num a => a -> a -> a
- Word64
ai
within :: Int
-> Double -> Double -> Bool
within :: Int -> Double -> Double -> Bool
within Int
ulps Double
a Double
b
| Int
ulps forall a. Ord a => a -> a -> Bool
< Int
0 = Bool
False
| Bool
otherwise = Double -> Double -> Word64
ulpDistance Double
a Double
b forall a. Ord a => a -> a -> Bool
<= forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ulps