-- https://tel.archives-ouvertes.fr/tel-03116750/document
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RecordWildCards #-}
module Data.Double.Shaman
  ( Shaman
  , significativeBits
  , significativeDigits
  , SDouble
  ) where

import Numeric.MathFunctions.Comparison ( addUlps)
import Numeric.MathFunctions.Constants
import GHC.TypeLits
import Data.Double.Approximate

-- | Double-precision floating point numbers that throw exceptions if
--   the accumulated errors grow large enough to cause unstable branching.
--
--   If @SDouble n@ works without throwing any exceptions, it'll be safe to
--   use @DoubleRelAbs n 0@ instead for a sizable performance boost.
--
-- >>> sin pi == (0 :: SDouble 0)
-- *** Exception: Insufficient precision.
-- ...
--
-- @SDouble 0@ failed so @DoubleRelAbs 0 0@ will lead to an unstable branch. In
-- other words, it'll return @False@ when it should have returned @True@:
--
-- >>> sin pi == (0 :: DoubleRelAbs 0 0)
-- False
--
-- Comparing to within 1 ULP stabalizes the branch:
--
-- >>> sin pi == (0 :: SDouble 1)
-- True
--
-- >>> sin pi == (0 :: DoubleRelAbs 1 0)
-- True
--
newtype SDouble (n::Nat) = SDouble Shaman
  deriving (Integer -> SDouble n
SDouble n -> SDouble n
SDouble n -> SDouble n -> SDouble n
(SDouble n -> SDouble n -> SDouble n)
-> (SDouble n -> SDouble n -> SDouble n)
-> (SDouble n -> SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (Integer -> SDouble n)
-> Num (SDouble n)
forall a.
(a -> a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (Integer -> a)
-> Num a
forall (n :: Nat). Integer -> SDouble n
forall (n :: Nat). SDouble n -> SDouble n
forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
fromInteger :: Integer -> SDouble n
$cfromInteger :: forall (n :: Nat). Integer -> SDouble n
signum :: SDouble n -> SDouble n
$csignum :: forall (n :: Nat). SDouble n -> SDouble n
abs :: SDouble n -> SDouble n
$cabs :: forall (n :: Nat). SDouble n -> SDouble n
negate :: SDouble n -> SDouble n
$cnegate :: forall (n :: Nat). SDouble n -> SDouble n
* :: SDouble n -> SDouble n -> SDouble n
$c* :: forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
- :: SDouble n -> SDouble n -> SDouble n
$c- :: forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
+ :: SDouble n -> SDouble n -> SDouble n
$c+ :: forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
Num, Num (SDouble n)
Num (SDouble n)
-> (SDouble n -> SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (Rational -> SDouble n)
-> Fractional (SDouble n)
Rational -> SDouble n
SDouble n -> SDouble n
SDouble n -> SDouble n -> SDouble n
forall a.
Num a
-> (a -> a -> a) -> (a -> a) -> (Rational -> a) -> Fractional a
forall (n :: Nat). Num (SDouble n)
forall (n :: Nat). Rational -> SDouble n
forall (n :: Nat). SDouble n -> SDouble n
forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
fromRational :: Rational -> SDouble n
$cfromRational :: forall (n :: Nat). Rational -> SDouble n
recip :: SDouble n -> SDouble n
$crecip :: forall (n :: Nat). SDouble n -> SDouble n
/ :: SDouble n -> SDouble n -> SDouble n
$c/ :: forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
$cp1Fractional :: forall (n :: Nat). Num (SDouble n)
Fractional, Fractional (SDouble n)
SDouble n
Fractional (SDouble n)
-> SDouble n
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n -> SDouble n)
-> (SDouble n -> SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> (SDouble n -> SDouble n)
-> Floating (SDouble n)
SDouble n -> SDouble n
SDouble n -> SDouble n -> SDouble n
forall a.
Fractional a
-> a
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a -> a)
-> (a -> a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> (a -> a)
-> Floating a
forall (n :: Nat). Fractional (SDouble n)
forall (n :: Nat). SDouble n
forall (n :: Nat). SDouble n -> SDouble n
forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
log1mexp :: SDouble n -> SDouble n
$clog1mexp :: forall (n :: Nat). SDouble n -> SDouble n
log1pexp :: SDouble n -> SDouble n
$clog1pexp :: forall (n :: Nat). SDouble n -> SDouble n
expm1 :: SDouble n -> SDouble n
$cexpm1 :: forall (n :: Nat). SDouble n -> SDouble n
log1p :: SDouble n -> SDouble n
$clog1p :: forall (n :: Nat). SDouble n -> SDouble n
atanh :: SDouble n -> SDouble n
$catanh :: forall (n :: Nat). SDouble n -> SDouble n
acosh :: SDouble n -> SDouble n
$cacosh :: forall (n :: Nat). SDouble n -> SDouble n
asinh :: SDouble n -> SDouble n
$casinh :: forall (n :: Nat). SDouble n -> SDouble n
tanh :: SDouble n -> SDouble n
$ctanh :: forall (n :: Nat). SDouble n -> SDouble n
cosh :: SDouble n -> SDouble n
$ccosh :: forall (n :: Nat). SDouble n -> SDouble n
sinh :: SDouble n -> SDouble n
$csinh :: forall (n :: Nat). SDouble n -> SDouble n
atan :: SDouble n -> SDouble n
$catan :: forall (n :: Nat). SDouble n -> SDouble n
acos :: SDouble n -> SDouble n
$cacos :: forall (n :: Nat). SDouble n -> SDouble n
asin :: SDouble n -> SDouble n
$casin :: forall (n :: Nat). SDouble n -> SDouble n
tan :: SDouble n -> SDouble n
$ctan :: forall (n :: Nat). SDouble n -> SDouble n
cos :: SDouble n -> SDouble n
$ccos :: forall (n :: Nat). SDouble n -> SDouble n
sin :: SDouble n -> SDouble n
$csin :: forall (n :: Nat). SDouble n -> SDouble n
logBase :: SDouble n -> SDouble n -> SDouble n
$clogBase :: forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
** :: SDouble n -> SDouble n -> SDouble n
$c** :: forall (n :: Nat). SDouble n -> SDouble n -> SDouble n
sqrt :: SDouble n -> SDouble n
$csqrt :: forall (n :: Nat). SDouble n -> SDouble n
log :: SDouble n -> SDouble n
$clog :: forall (n :: Nat). SDouble n -> SDouble n
exp :: SDouble n -> SDouble n
$cexp :: forall (n :: Nat). SDouble n -> SDouble n
pi :: SDouble n
$cpi :: forall (n :: Nat). SDouble n
$cp1Floating :: forall (n :: Nat). Fractional (SDouble n)
Floating, Num (SDouble n)
Ord (SDouble n)
Num (SDouble n)
-> Ord (SDouble n) -> (SDouble n -> Rational) -> Real (SDouble n)
SDouble n -> Rational
forall a. Num a -> Ord a -> (a -> Rational) -> Real a
forall (n :: Nat). KnownNat n => Num (SDouble n)
forall (n :: Nat). KnownNat n => Ord (SDouble n)
forall (n :: Nat). KnownNat n => SDouble n -> Rational
toRational :: SDouble n -> Rational
$ctoRational :: forall (n :: Nat). KnownNat n => SDouble n -> Rational
$cp2Real :: forall (n :: Nat). KnownNat n => Ord (SDouble n)
$cp1Real :: forall (n :: Nat). KnownNat n => Num (SDouble n)
Real, Fractional (SDouble n)
Real (SDouble n)
Real (SDouble n)
-> Fractional (SDouble n)
-> (forall b. Integral b => SDouble n -> (b, SDouble n))
-> (forall b. Integral b => SDouble n -> b)
-> (forall b. Integral b => SDouble n -> b)
-> (forall b. Integral b => SDouble n -> b)
-> (forall b. Integral b => SDouble n -> b)
-> RealFrac (SDouble n)
SDouble n -> b
SDouble n -> b
SDouble n -> b
SDouble n -> b
SDouble n -> (b, SDouble n)
forall b. Integral b => SDouble n -> b
forall b. Integral b => SDouble n -> (b, SDouble n)
forall a.
Real a
-> Fractional a
-> (forall b. Integral b => a -> (b, a))
-> (forall b. Integral b => a -> b)
-> (forall b. Integral b => a -> b)
-> (forall b. Integral b => a -> b)
-> (forall b. Integral b => a -> b)
-> RealFrac a
forall (n :: Nat). KnownNat n => Fractional (SDouble n)
forall (n :: Nat). KnownNat n => Real (SDouble n)
forall (n :: Nat) b. (KnownNat n, Integral b) => SDouble n -> b
forall (n :: Nat) b.
(KnownNat n, Integral b) =>
SDouble n -> (b, SDouble n)
floor :: SDouble n -> b
$cfloor :: forall (n :: Nat) b. (KnownNat n, Integral b) => SDouble n -> b
ceiling :: SDouble n -> b
$cceiling :: forall (n :: Nat) b. (KnownNat n, Integral b) => SDouble n -> b
round :: SDouble n -> b
$cround :: forall (n :: Nat) b. (KnownNat n, Integral b) => SDouble n -> b
truncate :: SDouble n -> b
$ctruncate :: forall (n :: Nat) b. (KnownNat n, Integral b) => SDouble n -> b
properFraction :: SDouble n -> (b, SDouble n)
$cproperFraction :: forall (n :: Nat) b.
(KnownNat n, Integral b) =>
SDouble n -> (b, SDouble n)
$cp2RealFrac :: forall (n :: Nat). KnownNat n => Fractional (SDouble n)
$cp1RealFrac :: forall (n :: Nat). KnownNat n => Real (SDouble n)
RealFrac, Floating (SDouble n)
RealFrac (SDouble n)
RealFrac (SDouble n)
-> Floating (SDouble n)
-> (SDouble n -> Integer)
-> (SDouble n -> Int)
-> (SDouble n -> (Int, Int))
-> (SDouble n -> (Integer, Int))
-> (Integer -> Int -> SDouble n)
-> (SDouble n -> Int)
-> (SDouble n -> SDouble n)
-> (Int -> SDouble n -> SDouble n)
-> (SDouble n -> Bool)
-> (SDouble n -> Bool)
-> (SDouble n -> Bool)
-> (SDouble n -> Bool)
-> (SDouble n -> Bool)
-> (SDouble n -> SDouble n -> SDouble n)
-> RealFloat (SDouble n)
Int -> SDouble n -> SDouble n
Integer -> Int -> SDouble n
SDouble n -> Bool
SDouble n -> Int
SDouble n -> Integer
SDouble n -> (Int, Int)
SDouble n -> (Integer, Int)
SDouble n -> SDouble n
SDouble n -> SDouble n -> SDouble n
forall a.
RealFrac a
-> Floating a
-> (a -> Integer)
-> (a -> Int)
-> (a -> (Int, Int))
-> (a -> (Integer, Int))
-> (Integer -> Int -> a)
-> (a -> Int)
-> (a -> a)
-> (Int -> a -> a)
-> (a -> Bool)
-> (a -> Bool)
-> (a -> Bool)
-> (a -> Bool)
-> (a -> Bool)
-> (a -> a -> a)
-> RealFloat a
forall (n :: Nat). KnownNat n => Floating (SDouble n)
forall (n :: Nat). KnownNat n => RealFrac (SDouble n)
forall (n :: Nat). KnownNat n => Int -> SDouble n -> SDouble n
forall (n :: Nat). KnownNat n => Integer -> Int -> SDouble n
forall (n :: Nat). KnownNat n => SDouble n -> Bool
forall (n :: Nat). KnownNat n => SDouble n -> Int
forall (n :: Nat). KnownNat n => SDouble n -> Integer
forall (n :: Nat). KnownNat n => SDouble n -> (Int, Int)
forall (n :: Nat). KnownNat n => SDouble n -> (Integer, Int)
forall (n :: Nat). KnownNat n => SDouble n -> SDouble n
forall (n :: Nat).
KnownNat n =>
SDouble n -> SDouble n -> SDouble n
atan2 :: SDouble n -> SDouble n -> SDouble n
$catan2 :: forall (n :: Nat).
KnownNat n =>
SDouble n -> SDouble n -> SDouble n
isIEEE :: SDouble n -> Bool
$cisIEEE :: forall (n :: Nat). KnownNat n => SDouble n -> Bool
isNegativeZero :: SDouble n -> Bool
$cisNegativeZero :: forall (n :: Nat). KnownNat n => SDouble n -> Bool
isDenormalized :: SDouble n -> Bool
$cisDenormalized :: forall (n :: Nat). KnownNat n => SDouble n -> Bool
isInfinite :: SDouble n -> Bool
$cisInfinite :: forall (n :: Nat). KnownNat n => SDouble n -> Bool
isNaN :: SDouble n -> Bool
$cisNaN :: forall (n :: Nat). KnownNat n => SDouble n -> Bool
scaleFloat :: Int -> SDouble n -> SDouble n
$cscaleFloat :: forall (n :: Nat). KnownNat n => Int -> SDouble n -> SDouble n
significand :: SDouble n -> SDouble n
$csignificand :: forall (n :: Nat). KnownNat n => SDouble n -> SDouble n
exponent :: SDouble n -> Int
$cexponent :: forall (n :: Nat). KnownNat n => SDouble n -> Int
encodeFloat :: Integer -> Int -> SDouble n
$cencodeFloat :: forall (n :: Nat). KnownNat n => Integer -> Int -> SDouble n
decodeFloat :: SDouble n -> (Integer, Int)
$cdecodeFloat :: forall (n :: Nat). KnownNat n => SDouble n -> (Integer, Int)
floatRange :: SDouble n -> (Int, Int)
$cfloatRange :: forall (n :: Nat). KnownNat n => SDouble n -> (Int, Int)
floatDigits :: SDouble n -> Int
$cfloatDigits :: forall (n :: Nat). KnownNat n => SDouble n -> Int
floatRadix :: SDouble n -> Integer
$cfloatRadix :: forall (n :: Nat). KnownNat n => SDouble n -> Integer
$cp2RealFloat :: forall (n :: Nat). KnownNat n => Floating (SDouble n)
$cp1RealFloat :: forall (n :: Nat). KnownNat n => RealFrac (SDouble n)
RealFloat)

instance Show (SDouble n) where
  showsPrec :: Int -> SDouble n -> ShowS
showsPrec Int
d (SDouble Shaman
s) = Int -> Shaman -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
d Shaman
s

instance Read (SDouble n) where
  readsPrec :: Int -> ReadS (SDouble n)
readsPrec Int
d String
inp = [ (Shaman -> SDouble n
forall (n :: Nat). Shaman -> SDouble n
SDouble Shaman
v, String
r) | (Shaman
v, String
r) <- Int -> ReadS Shaman
forall a. Read a => Int -> ReadS a
readsPrec Int
d String
inp ]


-- shamanUlpDistance :: Shaman -> Word64
-- shamanUlpDistance (Shaman n e) = ulpDistance n (n+e)

compareSDouble :: KnownNat n => SDouble n -> SDouble n -> Ordering
compareSDouble :: SDouble n -> SDouble n -> Ordering
compareSDouble a :: SDouble n
a@(SDouble Shaman
a') b :: SDouble n
b@(SDouble Shaman
b')
   | SDouble n
a SDouble n -> SDouble n -> Bool
forall a. Eq a => a -> a -> Bool
== SDouble n
b    = Ordering
EQ
   | Bool
otherwise = Shaman -> Double
shamanValue Shaman
a' Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Shaman -> Double
shamanValue Shaman
b'

toDoubleRelAbs :: SDouble n -> DoubleRelAbs n 0
toDoubleRelAbs :: SDouble n -> DoubleRelAbs n 0
toDoubleRelAbs (SDouble (Shaman Double
v Double
_e)) = Double -> DoubleRelAbs n 0
forall (abs :: Nat) (rel :: Nat). Double -> DoubleRelAbs abs rel
DoubleRelAbs Double
v

-- If (a :: DoubleRelAbs n 0) == (b :: DoubleRelAbs n 0)
instance KnownNat n => Eq (SDouble n) where
  a :: SDouble n
a@(SDouble Shaman
a') == :: SDouble n -> SDouble n -> Bool
== b :: SDouble n
b@(SDouble Shaman
b')
    | Bool
blessedResult Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== Bool
approxResult = Bool
blessedResult
    | Bool
otherwise                     = String -> Bool
forall a. HasCallStack => String -> a
error String
"Insufficient precision."
    where
      blessedResult :: Bool
blessedResult = Shaman
a' Shaman -> Shaman -> Bool
forall a. Eq a => a -> a -> Bool
== Shaman
b'
      approxResult :: Bool
approxResult = SDouble n -> DoubleRelAbs n 0
forall (n :: Nat). SDouble n -> DoubleRelAbs n 0
toDoubleRelAbs SDouble n
a DoubleRelAbs n 0 -> DoubleRelAbs n 0 -> Bool
forall a. Eq a => a -> a -> Bool
== SDouble n -> DoubleRelAbs n 0
forall (n :: Nat). SDouble n -> DoubleRelAbs n 0
toDoubleRelAbs SDouble n
b

instance KnownNat n => Ord (SDouble n) where
  compare :: SDouble n -> SDouble n -> Ordering
compare = SDouble n -> SDouble n -> Ordering
forall (n :: Nat). KnownNat n => SDouble n -> SDouble n -> Ordering
compareSDouble





-- | Double-precision floating point numbers with error-bounds.
--
-- Some digits can be represented exactly and have essentially an infinitely number of significant digits:
--
-- >>> significativeDigits 1
-- Infinity
--
-- Some fractional numbers can also be represented exactly:
--
-- >>> significativeDigits 0.5
-- Infinity
--
-- Other numbers are merely approximations:
--
-- >>> significativeDigits 0.1
-- 16.255619765854984
--
-- Pi is an irrational number so we can't represent it with infinite precision:
--
-- >>> significativeDigits pi
-- 15.849679651557175
--
-- @sin pi@ should theoretically be zero but we cannot do better than saying it is near zero:
--
-- >>> sin pi
-- 1.2246467991473532e-16
--
-- The error margins are greater than value itself so we have no significant digits:
--
-- >>> significativeDigits (sin pi)
-- 0.0
--
-- Since 'near zero' is not zero, the following fails when using Doubles:
--
-- >>> sin pi == (0 :: Double)
-- False
--
-- Equality testing for Shaman numbers tests whether the two intervals
-- overlap:
--
-- >>> sin pi == (0 :: Shaman)
-- True
data Shaman = Shaman
  { Shaman -> Double
shamanValue :: {-# UNPACK #-} !Double
  , Shaman -> Double
shamanError :: {-# UNPACK #-} !Double
  }

instance Show Shaman where
  showsPrec :: Int -> Shaman -> ShowS
showsPrec Int
d Shaman{Double
shamanError :: Double
shamanValue :: Double
shamanError :: Shaman -> Double
shamanValue :: Shaman -> Double
..} = Bool -> ShowS -> ShowS
showParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10) (ShowS -> ShowS) -> ShowS -> ShowS
forall a b. (a -> b) -> a -> b
$
    Double -> ShowS
forall a. Show a => a -> ShowS
shows Double
shamanValue ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
'±' ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> ShowS
forall a. Show a => a -> ShowS
shows Double
shamanError

instance Read Shaman where
  readsPrec :: Int -> ReadS Shaman
readsPrec Int
d = Bool -> ReadS Shaman -> ReadS Shaman
forall a. Bool -> ReadS a -> ReadS a
readParen (Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
app_prec) (ReadS Shaman -> ReadS Shaman) -> ReadS Shaman -> ReadS Shaman
forall a b. (a -> b) -> a -> b
$ \String
r ->
      [ (Shaman :: Double -> Double -> Shaman
Shaman{Double
shamanError :: Double
shamanValue :: Double
shamanError :: Double
shamanValue :: Double
..}, String
t')
      | (Double
shamanValue, Char
'±':String
t) <- ReadS Double
forall a. Read a => ReadS a
reads String
r
      , (Double
shamanError, String
t') <- ReadS Double
forall a. Read a => ReadS a
reads String
t]
    where app_prec :: Int
app_prec = Int
10

-- | Number of significant bits (base 2).
significativeBits :: Shaman -> Double
significativeBits :: Shaman -> Double
significativeBits Shaman
v = Shaman -> Double
significativeValue Shaman
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
log Double
2

-- | Number of significant digits (base 10).
significativeDigits :: Shaman -> Double
significativeDigits :: Shaman -> Double
significativeDigits Shaman
v = Shaman -> Double
significativeValue Shaman
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
log Double
10

significativeValue :: Shaman -> Double
significativeValue :: Shaman -> Double
significativeValue (Shaman Double
v Double
e)
  | Double
e Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = Double
m_pos_inf
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
e   = Double
0
  | Double
v Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double
forall a. Floating a => a -> a
log (Double -> Double
forall a. Num a => a -> a
abs Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
  | Bool
otherwise =
    let relError :: Double
relError = Double -> Double
forall a. Num a => a -> a
abs (Double
e Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
v) in
    if Double
relError Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1
      then Double
0
      else Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
log Double
relError

instance Num Shaman where
  Shaman Double
x Double
dx + :: Shaman -> Shaman -> Shaman
+ Shaman Double
y Double
dy =
    case Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
y of
      !Double
xy -> Double -> Double -> Shaman
Shaman Double
xy (Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
dyDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double -> Double -> Double -> Double
twoSum Double
x Double
y Double
xy)

  Shaman Double
x Double
dx * :: Shaman -> Shaman -> Shaman
* Shaman Double
y Double
dy =
    case Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y of
      !Double
z -> Double -> Double -> Shaman
Shaman Double
z (Double
dx Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dy Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double -> Double -> Double
fma Double
x Double
y (-Double
z))
  
  abs :: Shaman -> Shaman
abs (Shaman Double
x Double
dx) = Double -> Double -> Shaman
Shaman (Double -> Double
forall a. Num a => a -> a
abs Double
x) Double
dx

  signum :: Shaman -> Shaman
signum (Shaman Double
x Double
dx) = Double -> Double -> Shaman
Shaman (Double -> Double
forall a. Num a => a -> a
signum Double
x) Double
dx

  negate :: Shaman -> Shaman
negate (Shaman Double
x Double
dx) = Double -> Double -> Shaman
Shaman (Double -> Double
forall a. Num a => a -> a
negate Double
x) Double
dx

  fromInteger :: Integer -> Shaman
fromInteger Integer
i =
    let d :: Double
d = Integer -> Double
forall a. Num a => Integer -> a
fromInteger Integer
i
    in Double -> Double -> Shaman
Shaman Double
d (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Integer -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Integer -> Double) -> Integer -> Double
forall a b. (a -> b) -> a -> b
$ Integer
i Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Double -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Double
d)

instance Fractional Shaman where
  Shaman Double
x Double
dx / :: Shaman -> Shaman -> Shaman
/ Shaman Double
y Double
dy =
    let z :: Double
z = Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
y
        numerator :: Double
numerator = Double
dx Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double
fma Double
y Double
z (-Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dy
        denominator :: Double
denominator = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dy
    in Double -> Double -> Shaman
Shaman Double
z (Double
numerator Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
denominator)
  fromRational :: Rational -> Shaman
fromRational Rational
r =
    let d :: Double
d = Rational -> Double
forall a. Fractional a => Rational -> a
fromRational Rational
r
    in Double -> Double -> Shaman
Shaman Double
d (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Rational -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Rational -> Double) -> Rational -> Double
forall a b. (a -> b) -> a -> b
$ Rational
r Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a b. (Real a, Fractional b) => a -> b
realToFrac Double
d)

instance Floating Shaman where
  pi :: Shaman
pi = Double -> Double -> Shaman
Shaman Double
forall a. Floating a => a
pi (Int -> Double -> Double
addUlps Int
1 Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. Floating a => a
pi)
  sqrt :: Shaman -> Shaman
sqrt (Shaman Double
x Double
dx)
    | Double
z Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0, Double
dx Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 =
      Double -> Double -> Shaman
Shaman Double
z Double
0
    | Double
z Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 =
      Double -> Double -> Shaman
Shaman Double
z (Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double
forall a. Num a => a -> a
abs Double
dx)) -- FIXME: Use long double.
    | Bool
otherwise =
      Double -> Double -> Shaman
Shaman Double
z ((Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
dx) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
z))
    where
      r :: Double
r = Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> Double
fma Double
z Double
z (Double -> Double
forall a. Num a => a -> a
negate Double
x)
      z :: Double
z = Double -> Double
forall a. Floating a => a -> a
sqrt Double
x
  log :: Shaman -> Shaman
log = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_log Double -> Double
forall a. Floating a => a -> a
log
  exp :: Shaman -> Shaman
exp = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_exp Double -> Double
forall a. Floating a => a -> a
exp
  ** :: Shaman -> Shaman -> Shaman
(**) = (Double -> Double -> Double -> Double -> Double -> Double)
-> (Double -> Double -> Double) -> Shaman -> Shaman -> Shaman
foreignBinOp Double -> Double -> Double -> Double -> Double -> Double
erf_pow Double -> Double -> Double
forall a. Floating a => a -> a -> a
(**)
  sin :: Shaman -> Shaman
sin = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_sin Double -> Double
forall a. Floating a => a -> a
sin
  cos :: Shaman -> Shaman
cos = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_cos Double -> Double
forall a. Floating a => a -> a
cos
  tan :: Shaman -> Shaman
tan = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_tan Double -> Double
forall a. Floating a => a -> a
tan
  asin :: Shaman -> Shaman
asin = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_asin Double -> Double
forall a. Floating a => a -> a
asin
  acos :: Shaman -> Shaman
acos = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_acos Double -> Double
forall a. Floating a => a -> a
acos
  atan :: Shaman -> Shaman
atan = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_atan Double -> Double
forall a. Floating a => a -> a
atan
  sinh :: Shaman -> Shaman
sinh = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_sinh Double -> Double
forall a. Floating a => a -> a
sinh
  cosh :: Shaman -> Shaman
cosh = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_cosh Double -> Double
forall a. Floating a => a -> a
cosh
  tanh :: Shaman -> Shaman
tanh = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_tanh Double -> Double
forall a. Floating a => a -> a
tanh
  asinh :: Shaman -> Shaman
asinh = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_asinh Double -> Double
forall a. Floating a => a -> a
asinh
  acosh :: Shaman -> Shaman
acosh = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_acosh Double -> Double
forall a. Floating a => a -> a
acosh
  atanh :: Shaman -> Shaman
atanh = (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
erf_atanh Double -> Double
forall a. Floating a => a -> a
atanh

instance Real Shaman where
  toRational :: Shaman -> Rational
toRational (Shaman Double
v Double
_) = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
v

instance RealFrac Shaman where
  properFraction :: Shaman -> (b, Shaman)
properFraction (Shaman Double
v Double
e) =
    case Double -> (b, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction Double
v of
      (b
i, Double
v') -> (b
i, Double -> Double -> Shaman
Shaman Double
v' Double
e)
  truncate :: Shaman -> b
truncate (Shaman Double
v Double
_e) = Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
truncate Double
v
  round :: Shaman -> b
round (Shaman Double
v Double
_e) = Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round Double
v
  ceiling :: Shaman -> b
ceiling (Shaman Double
v Double
e) = Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double
vDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
e)
  floor :: Shaman -> b
floor (Shaman Double
v Double
e) = Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
vDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
e)

instance RealFloat Shaman where
  floatRadix :: Shaman -> Integer
floatRadix = Double -> Integer
forall a. RealFloat a => a -> Integer
floatRadix (Double -> Integer) -> (Shaman -> Double) -> Shaman -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  floatDigits :: Shaman -> Int
floatDigits = Double -> Int
forall a. RealFloat a => a -> Int
floatDigits (Double -> Int) -> (Shaman -> Double) -> Shaman -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  floatRange :: Shaman -> (Int, Int)
floatRange = Double -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange (Double -> (Int, Int))
-> (Shaman -> Double) -> Shaman -> (Int, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  decodeFloat :: Shaman -> (Integer, Int)
decodeFloat = Double -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat (Double -> (Integer, Int))
-> (Shaman -> Double) -> Shaman -> (Integer, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  encodeFloat :: Integer -> Int -> Shaman
encodeFloat Integer
m Int
e = Double -> Double -> Shaman
Shaman (Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
m Int
e) Double
0 -- FIXME: 1 ULP error bound?
  exponent :: Shaman -> Int
exponent = Double -> Int
forall a. RealFloat a => a -> Int
exponent (Double -> Int) -> (Shaman -> Double) -> Shaman -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  significand :: Shaman -> Shaman
significand Shaman
s = Double -> Double -> Shaman
Shaman (Double -> Double
forall a. RealFloat a => a -> a
significand (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Shaman -> Double
shamanValue Shaman
s) Double
0
  scaleFloat :: Int -> Shaman -> Shaman
scaleFloat Int
p (Shaman Double
v Double
e) = Double -> Double -> Shaman
Shaman (Int -> Double -> Double
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
p Double
v) Double
e
  isNaN :: Shaman -> Bool
isNaN = Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN (Double -> Bool) -> (Shaman -> Double) -> Shaman -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  isInfinite :: Shaman -> Bool
isInfinite = Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite (Double -> Bool) -> (Shaman -> Double) -> Shaman -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  isDenormalized :: Shaman -> Bool
isDenormalized = Double -> Bool
forall a. RealFloat a => a -> Bool
isDenormalized (Double -> Bool) -> (Shaman -> Double) -> Shaman -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  isNegativeZero :: Shaman -> Bool
isNegativeZero = Double -> Bool
forall a. RealFloat a => a -> Bool
isNegativeZero (Double -> Bool) -> (Shaman -> Double) -> Shaman -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  isIEEE :: Shaman -> Bool
isIEEE = Double -> Bool
forall a. RealFloat a => a -> Bool
isIEEE (Double -> Bool) -> (Shaman -> Double) -> Shaman -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Shaman -> Double
shamanValue
  atan2 :: Shaman -> Shaman -> Shaman
atan2 = (Double -> Double -> Double -> Double -> Double -> Double)
-> (Double -> Double -> Double) -> Shaman -> Shaman -> Shaman
foreignBinOp Double -> Double -> Double -> Double -> Double -> Double
erf_atan2 Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2


instance Eq Shaman where
  Shaman
x == :: Shaman -> Shaman -> Bool
== Shaman
y = Shaman -> Shaman -> Ordering
compareShaman Shaman
x Shaman
y Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
EQ

instance Ord Shaman where
  compare :: Shaman -> Shaman -> Ordering
compare = Shaman -> Shaman -> Ordering
compareShaman

-- Returns EQ iff the two numbers overlap.
compareShaman :: Shaman -> Shaman -> Ordering
compareShaman :: Shaman -> Shaman -> Ordering
compareShaman (Shaman Double
x Double
dx) (Shaman Double
y Double
dy)
  | Double -> Double
forall a. Num a => a -> a
abs (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double -> Double
forall a. Num a => a -> a
abs Double
dxDouble -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs Double
dy = Ordering
EQ
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
y             = Ordering
LT
  | Bool
otherwise         = Ordering
GT

twoSum :: Double -> Double -> Double -> Double
twoSum :: Double -> Double -> Double -> Double
twoSum Double
x Double
y Double
xy = Double -> Double
forall a. Num a => a -> a
abs (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x') Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Num a => a -> a
abs (Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y')
  where
    x' :: Double
x' = Double
xyDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y
    y' :: Double
y' = Double
xyDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x

foreignOp :: (Double -> Double -> Double -> Double) -> (Double -> Double) -> Shaman -> Shaman
foreignOp :: (Double -> Double -> Double -> Double)
-> (Double -> Double) -> Shaman -> Shaman
foreignOp Double -> Double -> Double -> Double
mkDelta Double -> Double
fn (Shaman Double
x Double
dx) =
  let !z :: Double
z = Double -> Double
fn Double
x
  in Double -> Double -> Shaman
Shaman Double
z (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> Double
mkDelta Double
x Double
dx Double
z)

foreignBinOp ::
    (Double -> Double -> Double -> Double -> Double -> Double) ->
    (Double -> Double -> Double) ->
    Shaman -> Shaman -> Shaman
foreignBinOp :: (Double -> Double -> Double -> Double -> Double -> Double)
-> (Double -> Double -> Double) -> Shaman -> Shaman -> Shaman
foreignBinOp Double -> Double -> Double -> Double -> Double -> Double
mkDelta Double -> Double -> Double
fn (Shaman Double
x Double
dx) (Shaman Double
y Double
dy) =
  let !z :: Double
z = Double -> Double -> Double
fn Double
x Double
y
  in Double -> Double -> Shaman
Shaman Double
z (Double -> Double -> Double -> Double -> Double -> Double
mkDelta Double
x Double
dx Double
y Double
dy Double
z)

foreign import ccall unsafe "fma" fma :: Double -> Double -> Double -> Double

foreign import ccall unsafe "erf_log" erf_log :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_exp" erf_exp :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_pow" erf_pow :: Double -> Double -> Double -> Double -> Double -> Double

foreign import ccall unsafe "erf_sin" erf_sin :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_cos" erf_cos :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_tan" erf_tan :: Double -> Double -> Double -> Double

foreign import ccall unsafe "erf_asin" erf_asin :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_acos" erf_acos :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_atan" erf_atan :: Double -> Double -> Double -> Double

foreign import ccall unsafe "erf_sinh" erf_sinh :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_cosh" erf_cosh :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_tanh" erf_tanh :: Double -> Double -> Double -> Double

foreign import ccall unsafe "erf_asinh" erf_asinh :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_acosh" erf_acosh :: Double -> Double -> Double -> Double
foreign import ccall unsafe "erf_atanh" erf_atanh :: Double -> Double -> Double -> Double

foreign import ccall unsafe "erf_atan2" erf_atan2 :: Double -> Double -> Double -> Double -> Double -> Double