-- SPDX-License-Identifier: MPL-2.0

{- | Numeric representations and deterministic IEEE utilities. Standard integer,
rational, binary32/64 and complex types are re-exported. Conversions never
silently promote precision: 'BigFloat' operations explicitly select the larger
operand precision, decimal/fixed operations the larger scale. The host IEEE
rounding environment is not modified. Explicit rounding returns software flags.
-}
module HMath.Numeric (
    Int8,
    Int16,
    Int32,
    Int64,
    Word8,
    Word16,
    Word32,
    Word64,
    Float32,
    Float64,
    Complex (..),
    Rational,
    RoundingMode (..),
    NumericException (..),
    Rounded (..),
    roundRational,
    IEEEClass (..),
    classify,
    NumericLimits (..),
    numericLimits,
    integerLimits,
    float16Limits,
    bfloat16Limits,
    nextAfter,
    frexp,
    ldexp,
    ulpDistance,
    withinULPs,
    safeIntegralCast,
    saturatingCast,
    realToIntegral,
    narrowFloat,
    Float16,
    float16,
    fromFloat16,
    float16Bits,
    float16FromBits,
    BFloat16,
    bfloat16,
    fromBFloat16,
    bfloat16Bits,
    bfloat16FromBits,
    BigFloat,
    bigFloat,
    bigPrecision,
    bigRational,
    bigToDouble,
    bigAdd,
    bigSubtract,
    bigMultiply,
    bigDivide,
    bigSqrt,
    bigExp,
    bigLog,
    extended,
    mixedAdd,
    mixedMultiply,
    Fixed,
    fixed,
    fixedScale,
    fixedCoefficient,
    fixedRational,
    fixedAdd,
    fixedMultiply,
    fixedDivide,
    Decimal,
    decimal,
    decimalRational,
    decimalAdd,
    decimalMultiply,
    decimalDivide,
    Interval,
    interval,
    intervalBounds,
    intervalAdd,
    intervalSubtract,
    intervalMultiply,
    intervalDivide,
    intervalContains,
    Dual (..),
    derivative,
    Quaternion (..),
    quaternionConjugate,
    quaternionNorm,
    quaternionInverse,
    rotateVector,
    DualQuaternion (..),
    dualQuaternionMultiply,
    dualQuaternionInverse,
    rigidTransform,
    transformPoint,
    kahanSum,
    neumaierSum,
    pairwiseSum,
    logsumexp,
    sumErrorBound,
) where

import Data.Bits (shiftL, shiftR, testBit, (.&.), (.|.))
import Data.Complex (Complex (..))
import Data.Int (Int16, Int32, Int64, Int8)
import qualified Data.List as L
import Data.Ratio (denominator, numerator, (%))
import Data.Word (Word16, Word32, Word64, Word8)
import GHC.Float (castDoubleToWord64)
import HMath.Error (MathError (..))

-- | Alias for Haskell Float (IEEE binary32 on supported targets).
type Float32 = Float

-- | Alias for Haskell Double (IEEE binary64 on supported targets).
type Float64 = Double

{- | Software rounding direction: nearest-even,toward zero,toward positive infinity,or
toward negative infinity. Host CPU rounding mode is unchanged.
-}
data RoundingMode = NearestEven | TowardZero | TowardPositive | TowardNegative
    deriving (RoundingMode -> RoundingMode -> Bool
(RoundingMode -> RoundingMode -> Bool)
-> (RoundingMode -> RoundingMode -> Bool) -> Eq RoundingMode
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RoundingMode -> RoundingMode -> Bool
== :: RoundingMode -> RoundingMode -> Bool
$c/= :: RoundingMode -> RoundingMode -> Bool
/= :: RoundingMode -> RoundingMode -> Bool
Eq, Int -> RoundingMode -> ShowS
[RoundingMode] -> ShowS
RoundingMode -> String
(Int -> RoundingMode -> ShowS)
-> (RoundingMode -> String)
-> ([RoundingMode] -> ShowS)
-> Show RoundingMode
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RoundingMode -> ShowS
showsPrec :: Int -> RoundingMode -> ShowS
$cshow :: RoundingMode -> String
show :: RoundingMode -> String
$cshowList :: [RoundingMode] -> ShowS
showList :: [RoundingMode] -> ShowS
Show)

{- | Software operation flags: invalid operation,division by zero,overflow,underflow,or
inexact conversion.
-}
data NumericException = InvalidOperation | DivisionByZero | NumericOverflow | NumericUnderflow | Inexact
    deriving (NumericException -> NumericException -> Bool
(NumericException -> NumericException -> Bool)
-> (NumericException -> NumericException -> Bool)
-> Eq NumericException
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: NumericException -> NumericException -> Bool
== :: NumericException -> NumericException -> Bool
$c/= :: NumericException -> NumericException -> Bool
/= :: NumericException -> NumericException -> Bool
Eq, Int -> NumericException -> ShowS
[NumericException] -> ShowS
NumericException -> String
(Int -> NumericException -> ShowS)
-> (NumericException -> String)
-> ([NumericException] -> ShowS)
-> Show NumericException
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> NumericException -> ShowS
showsPrec :: Int -> NumericException -> ShowS
$cshow :: NumericException -> String
show :: NumericException -> String
$cshowList :: [NumericException] -> ShowS
showList :: [NumericException] -> ShowS
Show)

{- | Rounded numeric result paired with software exception flags. Inspect flags when
intentional precision loss matters.
-}
data Rounded a = Rounded
    { forall a. Rounded a -> a
roundedValue :: a
    -- ^ Value produced by the explicitly selected software rounding operation.
    , forall a. Rounded a -> [NumericException]
numericalExceptions :: [NumericException]
    -- ^ Software rounding flags; these do not modify the host IEEE exception environment.
    }
    deriving (Rounded a -> Rounded a -> Bool
(Rounded a -> Rounded a -> Bool)
-> (Rounded a -> Rounded a -> Bool) -> Eq (Rounded a)
forall a. Eq a => Rounded a -> Rounded a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Rounded a -> Rounded a -> Bool
== :: Rounded a -> Rounded a -> Bool
$c/= :: forall a. Eq a => Rounded a -> Rounded a -> Bool
/= :: Rounded a -> Rounded a -> Bool
Eq, Int -> Rounded a -> ShowS
[Rounded a] -> ShowS
Rounded a -> String
(Int -> Rounded a -> ShowS)
-> (Rounded a -> String)
-> ([Rounded a] -> ShowS)
-> Show (Rounded a)
forall a. Show a => Int -> Rounded a -> ShowS
forall a. Show a => [Rounded a] -> ShowS
forall a. Show a => Rounded a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Rounded a -> ShowS
showsPrec :: Int -> Rounded a -> ShowS
$cshow :: forall a. Show a => Rounded a -> String
show :: Rounded a -> String
$cshowList :: forall a. Show a => [Rounded a] -> ShowS
showList :: [Rounded a] -> ShowS
Show)

{- | Given software rounding mode then exact Rational, return Integer rounded by
nearest-even, toward zero, toward positive infinity, or toward negative infinity. No
host rounding flags are changed.
-}
roundRational :: RoundingMode -> Rational -> Integer
roundRational :: RoundingMode -> Rational -> Integer
roundRational RoundingMode
NearestEven = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
round
roundRational RoundingMode
TowardZero = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
truncate
roundRational RoundingMode
TowardPositive = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling
roundRational RoundingMode
TowardNegative = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor

{- | IEEE classification distinguishing NaN,signed infinity,signed zero,subnormal,and normal
finite values.
-}
data IEEEClass
    = NotANumber
    | PositiveInfinity
    | NegativeInfinity
    | PositiveZero
    | NegativeZero
    | Subnormal
    | Normal
    deriving (IEEEClass -> IEEEClass -> Bool
(IEEEClass -> IEEEClass -> Bool)
-> (IEEEClass -> IEEEClass -> Bool) -> Eq IEEEClass
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IEEEClass -> IEEEClass -> Bool
== :: IEEEClass -> IEEEClass -> Bool
$c/= :: IEEEClass -> IEEEClass -> Bool
/= :: IEEEClass -> IEEEClass -> Bool
Eq, Int -> IEEEClass -> ShowS
[IEEEClass] -> ShowS
IEEEClass -> String
(Int -> IEEEClass -> ShowS)
-> (IEEEClass -> String)
-> ([IEEEClass] -> ShowS)
-> Show IEEEClass
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IEEEClass -> ShowS
showsPrec :: Int -> IEEEClass -> ShowS
$cshow :: IEEEClass -> String
show :: IEEEClass -> String
$cshowList :: [IEEEClass] -> ShowS
showList :: [IEEEClass] -> ShowS
Show)

{- | Classify a RealFloat value as NaN, signed infinity, signed zero, subnormal, or normal.
Signed zero is distinguished with isNegativeZero.
-}
classify :: (RealFloat a) => a -> IEEEClass
classify :: forall a. RealFloat a => a -> IEEEClass
classify a
x
    | a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x = IEEEClass
NotANumber
    | a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x = if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 then IEEEClass
PositiveInfinity else IEEEClass
NegativeInfinity
    | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = if a -> Bool
forall a. RealFloat a => a -> Bool
isNegativeZero a
x then IEEEClass
NegativeZero else IEEEClass
PositiveZero
    | a -> Bool
forall a. RealFloat a => a -> Bool
isDenormalized a
x = IEEEClass
Subnormal
    | Bool
otherwise = IEEEClass
Normal

{- | Format limits and precision metadata derived from RealFloat or the explicit 16-bit
format definition.
-}
data NumericLimits a = NumericLimits
    { forall a. NumericLimits a -> a
minimumNormal :: a
    -- ^ Smallest positive normal value of the represented floating format.
    , forall a. NumericLimits a -> a
minimumSubnormal :: a
    -- ^ Smallest positive subnormal value of the represented floating format.
    , forall a. NumericLimits a -> a
maximumFinite :: a
    -- ^ Largest finite positive value of the represented floating format.
    , forall a. NumericLimits a -> a
machineEpsilon :: a
    -- ^ Gap between one and its next larger representable value.
    , forall a. NumericLimits a -> Int
significandDigits :: Int
    -- ^ Number of radix digits in the floating significand.
    , forall a. NumericLimits a -> Int
minimumExponent :: Int
    -- ^ Lower exponent bound returned by RealFloat.floatRange.
    , forall a. NumericLimits a -> Int
maximumExponent :: Int
    -- ^ Upper exponent bound returned by RealFloat.floatRange.
    }
    deriving (NumericLimits a -> NumericLimits a -> Bool
(NumericLimits a -> NumericLimits a -> Bool)
-> (NumericLimits a -> NumericLimits a -> Bool)
-> Eq (NumericLimits a)
forall a. Eq a => NumericLimits a -> NumericLimits a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => NumericLimits a -> NumericLimits a -> Bool
== :: NumericLimits a -> NumericLimits a -> Bool
$c/= :: forall a. Eq a => NumericLimits a -> NumericLimits a -> Bool
/= :: NumericLimits a -> NumericLimits a -> Bool
Eq, Int -> NumericLimits a -> ShowS
[NumericLimits a] -> ShowS
NumericLimits a -> String
(Int -> NumericLimits a -> ShowS)
-> (NumericLimits a -> String)
-> ([NumericLimits a] -> ShowS)
-> Show (NumericLimits a)
forall a. Show a => Int -> NumericLimits a -> ShowS
forall a. Show a => [NumericLimits a] -> ShowS
forall a. Show a => NumericLimits a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> NumericLimits a -> ShowS
showsPrec :: Int -> NumericLimits a -> ShowS
$cshow :: forall a. Show a => NumericLimits a -> String
show :: NumericLimits a -> String
$cshowList :: forall a. Show a => [NumericLimits a] -> ShowS
showList :: [NumericLimits a] -> ShowS
Show)

-- | Derived from 'RealFloat' radix/digits/range, not hardcoded for Double.
numericLimits :: (RealFloat a) => a -> NumericLimits a
numericLimits :: forall a. RealFloat a => a -> NumericLimits a
numericLimits a
witness =
    a -> a -> a -> a -> Int -> Int -> Int -> NumericLimits a
forall a. a -> a -> a -> a -> Int -> Int -> Int -> NumericLimits a
NumericLimits
        (Integer -> Int -> a
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
1 (Int
lo Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
        (Integer -> Int -> a
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
1 (Int
lo Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p))
        (Integer -> Int -> a
forall a. RealFloat a => Integer -> Int -> a
encodeFloat (Integer
b Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) (Int
hi Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p))
        (Integer -> Int -> a
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
1 (Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p))
        Int
p
        Int
lo
        Int
hi
  where
    b :: Integer
b = Integer -> Integer
forall a. Integral a => a -> Integer
toInteger (a -> Integer
forall a. RealFloat a => a -> Integer
floatRadix a
witness); p :: Int
p = a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
witness; (Int
lo, Int
hi) = a -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange a
witness

-- | Limits returned as exact binary64 values for the two 16-bit storage formats.
float16Limits, bfloat16Limits :: NumericLimits Double
float16Limits :: NumericLimits Double
float16Limits = Double
-> Double
-> Double
-> Double
-> Int
-> Int
-> Int
-> NumericLimits Double
forall a. a -> a -> a -> a -> Int -> Int -> Int -> NumericLimits a
NumericLimits (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
14)) (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
24)) Double
65504 (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
10)) Int
11 (-Int
13) Int
16
bfloat16Limits :: NumericLimits Double
bfloat16Limits = Double
-> Double
-> Double
-> Double
-> Int
-> Int
-> Int
-> NumericLimits Double
forall a. a -> a -> a -> a -> Int -> Int -> Int -> NumericLimits a
NumericLimits (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
126)) (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
133)) (Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
255 Int
120) (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
7)) Int
8 (-Int
125) Int
128

{- | Given a bounded integral type witness, return its exact (minimum,maximum) as Integers.
The witness value is used only to select the type.
-}
integerLimits :: (Integral a, Bounded a) => a -> (Integer, Integer)
integerLimits :: forall a. (Integral a, Bounded a) => a -> (Integer, Integer)
integerLimits a
x = (a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
forall a. Bounded a => a
minBound a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
x), a -> Integer
forall a. Integral a => a -> Integer
toInteger (a
forall a. Bounded a => a
maxBound a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
x))

{- | Adjacent representable value toward target, including gradual underflow.
Equal signed zeros return the target; NaN propagates. Supports binary formats.
-}
nextAfter :: (RealFloat a) => a -> a -> a
nextAfter :: forall a. RealFloat a => a -> a -> a
nextAfter a
x a
target
    | a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
target = a
0 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0
    | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
target = a
target
    | a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x = if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 then NumericLimits a -> a
forall a. NumericLimits a -> a
maximumFinite NumericLimits a
lim else a -> a
forall a. Num a => a -> a
negate (NumericLimits a -> a
forall a. NumericLimits a -> a
maximumFinite NumericLimits a
lim)
    | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = if a
target a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 then NumericLimits a -> a
forall a. NumericLimits a -> a
minimumSubnormal NumericLimits a
lim else a -> a
forall a. Num a => a -> a
negate (NumericLimits a -> a
forall a. NumericLimits a -> a
minimumSubnormal NumericLimits a
lim)
    | Bool
otherwise = a -> a
forall a. Num a => a -> a
signum a
x a -> a -> a
forall a. Num a => a -> a -> a
* Rational -> a
forall a. Fractional a => Rational -> a
fromRational (a -> Rational
forall a. Real a => a -> Rational
toRational a
ax Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ (if Bool
increase then Rational
unit else Rational -> Rational
forall a. Num a => a -> a
negate Rational
unit))
  where
    lim :: NumericLimits a
lim = a -> NumericLimits a
forall a. RealFloat a => a -> NumericLimits a
numericLimits a
x
    ax :: a
ax = a -> a
forall a. Num a => a -> a
abs a
x
    increase :: Bool
increase = (a
target a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
x) Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== (a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0)
    (Integer
m, Int
e) = a -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat a
ax
    p :: Int
p = a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
x
    b :: Integer
b = Integer -> Integer
forall a. Integral a => a -> Integer
toInteger (a -> Integer
forall a. RealFloat a => a -> Integer
floatRadix a
x)
    quantum :: Int
quantum = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (NumericLimits a -> Int
forall a. NumericLimits a -> Int
minimumExponent NumericLimits a
lim Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p) (if Bool -> Bool
not Bool
increase Bool -> Bool -> Bool
&& Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
b Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) then Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 else Int
e)
    unit :: Rational
unit = if Int
quantum Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
b Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
quantum) else Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
b Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int -> Int
forall a. Num a => a -> a
negate Int
quantum)

-- | x=m*2^e, with abs m in [1/2,1) for nonzero finite binary inputs.
frexp :: (RealFloat a) => a -> (a, Int)
frexp :: forall a. RealFloat a => a -> (a, Int)
frexp a
x
    | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x = (a
x, Int
0)
    | Bool
otherwise = let (Integer
_, Int
e) = a -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat a
x; k :: Int
k = Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
x in (Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (-Int
k) a
x, Int
k)

{- | Given real value then Int exponent, scale by the numeric radix to that power using
scaleFloat (binary types use powers of two). IEEE overflow/underflow follows the host
type; no Either check is added.
-}
ldexp :: (RealFloat a) => a -> Int -> a
ldexp :: forall a. RealFloat a => a -> Int -> a
ldexp a
x Int
n = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
n a
x

-- | IEEE binary64 distance. Signed zeros compare equal; NaN is rejected.
ulpDistance :: Double -> Double -> Either MathError Integer
ulpDistance :: Double -> Double -> Either MathError Integer
ulpDistance Double
a Double
b
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
a Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
b = MathError -> Either MathError Integer
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"ULP distance is undefined for NaN")
    | Bool
otherwise = Integer -> Either MathError Integer
forall a b. b -> Either a b
Right (Integer -> Integer
forall a. Num a => a -> a
abs (Double -> Integer
ordered Double
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Double -> Integer
ordered Double
b))
  where
    ordered :: Double -> Integer
ordered Double
x =
        let bits :: Word64
bits = Double -> Word64
castDoubleToWord64 Double
x; mag :: Integer
mag = Word64 -> Integer
forall a. Integral a => a -> Integer
toInteger (Word64
bits Word64 -> Word64 -> Word64
forall a. Bits a => a -> a -> a
.&. Word64
0x7fffffffffffffff)
         in if Word64 -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Word64
bits Int
63 then Integer -> Integer
forall a. Num a => a -> a
negate Integer
mag else Integer
mag

{- | Given nonnegative allowed ULP distance then two Doubles, test their binary64
representable-value distance. Signed zeros compare equal. Negative tolerance or NaN
gives False.
-}
withinULPs :: Integer -> Double -> Double -> Bool
withinULPs :: Integer -> Double -> Double -> Bool
withinULPs Integer
n Double
a Double
b = Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
0 Bool -> Bool -> Bool
&& (MathError -> Bool)
-> (Integer -> Bool) -> Either MathError Integer -> Bool
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Bool -> MathError -> Bool
forall a b. a -> b -> a
const Bool
False) (Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
n) (Double -> Double -> Either MathError Integer
ulpDistance Double
a Double
b)

{- | Convert an integral value to a bounded integral destination only if in range. Return
Left DomainError rather than wrap or truncate out-of-range input.
-}
safeIntegralCast :: (Integral a, Integral b, Bounded b) => a -> Either MathError b
safeIntegralCast :: forall a b.
(Integral a, Integral b, Bounded b) =>
a -> Either MathError b
safeIntegralCast a
x = Either MathError b
r
  where
    y :: b
y = a -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
x
    r :: Either MathError b
r =
        if a -> Integer
forall a. Integral a => a -> Integer
toInteger a
x Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< b -> Integer
forall a. Integral a => a -> Integer
toInteger (b
forall a. Bounded a => a
minBound b -> b -> b
forall a. a -> a -> a
`asTypeOf` b
y) Bool -> Bool -> Bool
|| a -> Integer
forall a. Integral a => a -> Integer
toInteger a
x Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> b -> Integer
forall a. Integral a => a -> Integer
toInteger (b
forall a. Bounded a => a
maxBound b -> b -> b
forall a. a -> a -> a
`asTypeOf` b
y)
            then MathError -> Either MathError b
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"integral cast out of range")
            else b -> Either MathError b
forall a b. b -> Either a b
Right b
y

{- | Convert an integral value to a bounded integral destination, clamping to
minBound/maxBound when out of range. No error result is produced.
-}
saturatingCast :: (Integral a, Integral b, Bounded b) => a -> b
saturatingCast :: forall a b. (Integral a, Integral b, Bounded b) => a -> b
saturatingCast a
x = b
y
  where
    y :: b
y =
        Integer -> b
forall a. Num a => Integer -> a
fromInteger
            (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max (b -> Integer
forall a. Integral a => a -> Integer
toInteger (b
forall a. Bounded a => a
minBound b -> b -> b
forall a. a -> a -> a
`asTypeOf` b
y)) (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min (b -> Integer
forall a. Integral a => a -> Integer
toInteger (b
forall a. Bounded a => a
maxBound b -> b -> b
forall a. a -> a -> a
`asTypeOf` b
y)) (a -> Integer
forall a. Integral a => a -> Integer
toInteger a
x)))

{- | Given software rounding mode then finite real value, round exactly through Rational
then range-check the bounded integral destination. Non-finite or out-of-range input
returns Left DomainError.
-}
realToIntegral :: (RealFloat a, Integral b, Bounded b) => RoundingMode -> a -> Either MathError b
realToIntegral :: forall a b.
(RealFloat a, Integral b, Bounded b) =>
RoundingMode -> a -> Either MathError b
realToIntegral RoundingMode
mode a
x
    | a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x = MathError -> Either MathError b
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"non-finite integral cast")
    | Bool
otherwise = Integer -> Either MathError b
forall a b.
(Integral a, Integral b, Bounded b) =>
a -> Either MathError b
safeIntegralCast (RoundingMode -> Rational -> Integer
roundRational RoundingMode
mode (a -> Rational
forall a. Real a => a -> Rational
toRational a
x))

{- | A narrowing conversion must be finite and exact; choose explicit quantized
Float16/BigFloat rounding APIs when loss of precision is intentional.
-}
narrowFloat :: (RealFloat a, RealFloat b) => a -> Either MathError b
narrowFloat :: forall a b. (RealFloat a, RealFloat b) => a -> Either MathError b
narrowFloat a
x
    | a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x = MathError -> Either MathError b
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"non-finite float cast")
    | b -> Bool
forall a. RealFloat a => a -> Bool
isInfinite b
y = MathError -> Either MathError b
forall a b. a -> Either a b
Left MathError
NonFiniteResult
    | b -> Rational
forall a. Real a => a -> Rational
toRational b
y Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= a -> Rational
forall a. Real a => a -> Rational
toRational a
x = MathError -> Either MathError b
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"inexact float cast")
    | Bool
otherwise = b -> Either MathError b
forall a b. b -> Either a b
Right b
y
  where
    y :: b
y = if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
isNegativeZero a
x then -b
0 else Rational -> b
forall a. Fractional a => Rational -> a
fromRational (a -> Rational
forall a. Real a => a -> Rational
toRational a
x)

{- | IEEE binary16 and bfloat16 storage. Conversion rounds directly from the
exact binary64 value (no intermediate binary32 rounding), preserving signed
zero and infinity. NaNs are canonical quiet NaNs. All four rounding modes
report inexact, overflow and tininess-after-rounding underflow flags.
-}
newtype Float16 = Float16 Word16 deriving (Float16 -> Float16 -> Bool
(Float16 -> Float16 -> Bool)
-> (Float16 -> Float16 -> Bool) -> Eq Float16
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Float16 -> Float16 -> Bool
== :: Float16 -> Float16 -> Bool
$c/= :: Float16 -> Float16 -> Bool
/= :: Float16 -> Float16 -> Bool
Eq, Int -> Float16 -> ShowS
[Float16] -> ShowS
Float16 -> String
(Int -> Float16 -> ShowS)
-> (Float16 -> String) -> ([Float16] -> ShowS) -> Show Float16
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Float16 -> ShowS
showsPrec :: Int -> Float16 -> ShowS
$cshow :: Float16 -> String
show :: Float16 -> String
$cshowList :: [Float16] -> ShowS
showList :: [Float16] -> ShowS
Show)

{- | Opaque bfloat16 storage with 8 exponent and 7 explicit fraction bits. Use bfloat16 for
rounded conversion or bfloat16FromBits for exact storage construction.
-}
newtype BFloat16 = BFloat16 Word16 deriving (BFloat16 -> BFloat16 -> Bool
(BFloat16 -> BFloat16 -> Bool)
-> (BFloat16 -> BFloat16 -> Bool) -> Eq BFloat16
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BFloat16 -> BFloat16 -> Bool
== :: BFloat16 -> BFloat16 -> Bool
$c/= :: BFloat16 -> BFloat16 -> Bool
/= :: BFloat16 -> BFloat16 -> Bool
Eq, Int -> BFloat16 -> ShowS
[BFloat16] -> ShowS
BFloat16 -> String
(Int -> BFloat16 -> ShowS)
-> (BFloat16 -> String) -> ([BFloat16] -> ShowS) -> Show BFloat16
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BFloat16 -> ShowS
showsPrec :: Int -> BFloat16 -> ShowS
$cshow :: BFloat16 -> String
show :: BFloat16 -> String
$cshowList :: [BFloat16] -> ShowS
showList :: [BFloat16] -> ShowS
Show)

{- | Return raw IEEE binary16 storage bits, including sign, exponent and fraction, as
Word16.
-}
float16Bits :: Float16 -> Word16
float16Bits :: Float16 -> Word16
float16Bits (Float16 Word16
w) = Word16
w

{- | Construct binary16 storage from raw Word16 without validation or rounding; all bit
patterns, including NaNs, are representable.
-}
float16FromBits :: Word16 -> Float16
float16FromBits :: Word16 -> Float16
float16FromBits = Word16 -> Float16
Float16

-- | Return raw bfloat16 storage bits as Word16.
bfloat16Bits :: BFloat16 -> Word16
bfloat16Bits :: BFloat16 -> Word16
bfloat16Bits (BFloat16 Word16
w) = Word16
w

{- | Construct bfloat16 storage from raw Word16 without conversion; all bit patterns are
retained.
-}
bfloat16FromBits :: Word16 -> BFloat16
bfloat16FromBits :: Word16 -> BFloat16
bfloat16FromBits = Word16 -> BFloat16
BFloat16

{- | Given rounding mode then Double, round directly to IEEE binary16 and return Rounded
storage plus software flags. Preserve signed zero/infinity; canonicalize NaN. Flags
report inexact, overflow, and tininess-after-rounding underflow.
-}
float16 :: RoundingMode -> Double -> Rounded Float16
float16 :: RoundingMode -> Double -> Rounded Float16
float16 RoundingMode
mode Double
x = let Rounded Word16
w [NumericException]
flags = Int -> Int -> RoundingMode -> Double -> Rounded Word16
encodeSmall Int
5 Int
10 RoundingMode
mode Double
x in Float16 -> [NumericException] -> Rounded Float16
forall a. a -> [NumericException] -> Rounded a
Rounded (Word16 -> Float16
Float16 Word16
w) [NumericException]
flags

{- | Given rounding mode then Double, round directly to bfloat16 without an intermediate
binary32 rounding. Return Rounded storage and software flags, preserving signed
zero/infinity and canonicalizing NaN.
-}
bfloat16 :: RoundingMode -> Double -> Rounded BFloat16
bfloat16 :: RoundingMode -> Double -> Rounded BFloat16
bfloat16 RoundingMode
mode Double
x = let Rounded Word16
w [NumericException]
flags = Int -> Int -> RoundingMode -> Double -> Rounded Word16
encodeSmall Int
8 Int
7 RoundingMode
mode Double
x in BFloat16 -> [NumericException] -> Rounded BFloat16
forall a. a -> [NumericException] -> Rounded a
Rounded (Word16 -> BFloat16
BFloat16 Word16
w) [NumericException]
flags

{- | Decode binary16 storage to Double, preserving finite values, signed zero, infinity and
NaN classification. Every finite binary16 value is exactly representable as Double.
-}
fromFloat16 :: Float16 -> Double
fromFloat16 :: Float16 -> Double
fromFloat16 (Float16 Word16
w) = Int -> Int -> Word16 -> Double
decodeSmall Int
5 Int
10 Word16
w

{- | Decode bfloat16 storage to Double, preserving finite values, signed zero and non-finite
classification.
-}
fromBFloat16 :: BFloat16 -> Double
fromBFloat16 :: BFloat16 -> Double
fromBFloat16 (BFloat16 Word16
w) = Int -> Int -> Word16 -> Double
decodeSmall Int
8 Int
7 Word16
w

powerTwo :: Int -> Rational
powerTwo :: Int -> Rational
powerTwo Int
e
    | Int
e Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
e)
    | Bool
otherwise = Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int -> Int
forall a. Num a => a -> a
negate Int
e)
integerBits :: Integer -> Int
integerBits :: Integer -> Int
integerBits Integer
n = Integer -> Int -> Int
forall {t} {t}. (Integral t, Num t) => t -> t -> t
go (Integer -> Integer
forall a. Num a => a -> a
abs Integer
n) Int
0 where go :: t -> t -> t
go t
0 t
k = t
k; go t
v t
k = t -> t -> t
go (t
v t -> t -> t
forall a. Integral a => a -> a -> a
`quot` t
2) (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
floorLog2 :: Rational -> Int
floorLog2 :: Rational -> Int
floorLog2 Rational
x =
    let candidate :: Int
candidate = Integer -> Int
integerBits (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
x) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Integer -> Int
integerBits (Rational -> Integer
forall a. Ratio a -> a
denominator Rational
x)
     in if Rational
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Rational
powerTwo Int
candidate then Int
candidate Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 else Int
candidate
encodeSmall :: Int -> Int -> RoundingMode -> Double -> Rounded Word16
encodeSmall :: Int -> Int -> RoundingMode -> Double -> Rounded Word16
encodeSmall Int
eb Int
fb RoundingMode
mode Double
x
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x =
        Word16 -> [NumericException] -> Rounded Word16
forall a. a -> [NumericException] -> Rounded a
Rounded
            (Word16
expMask Word16 -> Word16 -> Word16
forall a. Bits a => a -> a -> a
.|. Int -> Word16
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
fb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) :: Int))
            [NumericException
InvalidOperation | Bool -> Bool
not (Word64 -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit (Double -> Word64
castDoubleToWord64 Double
x) Int
51)]
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x = Word16 -> [NumericException] -> Rounded Word16
forall a. a -> [NumericException] -> Rounded a
Rounded (Word16
sgn Word16 -> Word16 -> Word16
forall a. Bits a => a -> a -> a
.|. Word16
expMask) []
    | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Word16 -> [NumericException] -> Rounded Word16
forall a. a -> [NumericException] -> Rounded a
Rounded Word16
sgn []
    | Bool
otherwise =
        let q :: Rational
q = Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x
            e :: Int
e = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
emin (Rational -> Int
floorLog2 (Rational -> Rational
forall a. Num a => a -> a
abs Rational
q))
            quantum :: Rational
quantum = Int -> Rational
powerTwo (Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
fb)
            mant :: Integer
mant = RoundingMode -> Rational -> Integer
roundRational RoundingMode
mode (Rational
q Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
quantum)
            rounded :: Rational
rounded = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
mant Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
quantum
            inexact :: Bool
inexact = Rational
rounded Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
q
            maxValue :: Rational
maxValue = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
fb Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Int -> Rational
powerTwo (Int
emax Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
fb)
            overflow :: Bool
overflow = Rational -> Rational
forall a. Num a => a -> a
abs Rational
rounded Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
maxValue
            toInf :: Bool
toInf = RoundingMode
mode RoundingMode -> RoundingMode -> Bool
forall a. Eq a => a -> a -> Bool
== RoundingMode
NearestEven Bool -> Bool -> Bool
|| (RoundingMode
mode RoundingMode -> RoundingMode -> Bool
forall a. Eq a => a -> a -> Bool
== RoundingMode
TowardPositive Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) Bool -> Bool -> Bool
|| (RoundingMode
mode RoundingMode -> RoundingMode -> Bool
forall a. Eq a => a -> a -> Bool
== RoundingMode
TowardNegative Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0)
            magnitude :: Integer
magnitude = Integer -> Integer
forall a. Num a => a -> a
abs Integer
mant
            (Int
expField, Integer
frac) =
                if Integer
magnitude Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
fb
                    then (Int
0, Integer
magnitude)
                    else if Integer
magnitude Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
fb Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) then (Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
bias, Integer
0) else (Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
bias, Integer
magnitude Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
fb)
            bits :: Word16
bits =
                if Bool
overflow
                    then (if Bool
toInf then Word16
expMask else Word16
expMask Word16 -> Word16 -> Word16
forall a. Num a => a -> a -> a
- Word16
1)
                    else Integer -> Word16
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
expField Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftL` Int
fb) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
frac)
            flags :: [NumericException]
flags =
                [NumericException
NumericOverflow | Bool
overflow]
                    [NumericException] -> [NumericException] -> [NumericException]
forall a. [a] -> [a] -> [a]
++ [NumericException
NumericUnderflow | Bool
inexact Bool -> Bool -> Bool
&& Rational -> Rational
forall a. Num a => a -> a
abs Rational
rounded Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Rational
powerTwo Int
emin]
                    [NumericException] -> [NumericException] -> [NumericException]
forall a. [a] -> [a] -> [a]
++ [NumericException
Inexact | Bool
inexact Bool -> Bool -> Bool
|| Bool
overflow]
         in Word16 -> [NumericException] -> Rounded Word16
forall a. a -> [NumericException] -> Rounded a
Rounded (Word16
sgn Word16 -> Word16 -> Word16
forall a. Bits a => a -> a -> a
.|. Word16
bits) [NumericException]
flags
  where
    bias :: Int
bias = Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
eb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
    emin :: Int
emin = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
bias
    emax :: Int
emax = Int
bias
    expMask :: Word16
expMask = Int -> Word16
forall a b. (Integral a, Num b) => a -> b
fromIntegral (((Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
eb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) :: Int) Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftL` Int
fb)
    sgn :: Word16
sgn = if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNegativeZero Double
x then Word16
0x8000 else Word16
0

decodeSmall :: Int -> Int -> Word16 -> Double
decodeSmall :: Int -> Int -> Word16 -> Double
decodeSmall Int
eb Int
fb Word16
w
    | Int
e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
eb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 = if Int
f Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
s Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 else Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
    | Int
e Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 =
        if Int
f Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then (if Word16 -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Word16
w Int
15 then -Double
0 else Double
0) else Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
f) (Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
bias Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
fb)
    | Bool
otherwise = Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat (Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
fb Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
f)) (Int
e Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
bias Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
fb)
  where
    e :: Int
e = Word16 -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Word16
w Word16 -> Int -> Word16
forall a. Bits a => a -> Int -> a
`shiftR` Int
fb) Word16 -> Word16 -> Word16
forall a. Bits a => a -> a -> a
.&. Int -> Word16
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
eb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) :: Int))
    f :: Int
f = Word16 -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word16
w Word16 -> Word16 -> Word16
forall a. Bits a => a -> a -> a
.&. Int -> Word16
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
fb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) :: Int)) :: Int
    bias :: Int
bias = Int
2 Int -> Int -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
eb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
    s :: Double
s = if Word16 -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Word16
w Int
15 then -Double
1 else Double
1

{- | An unbounded-exponent binary floating value backed by an exact Rational,
rounded to a runtime number of significant bits (2..1,000,000). Unlike Double,
this is genuinely configurable precision. No NaN/infinity or hidden CPU flags.
Algebraic operations round once, nearest-even; explicit construction can select
any 'RoundingMode'. Transcendentals use guard bits and rational series.
-}
data BigFloat = BigFloat Int Rational deriving (BigFloat -> BigFloat -> Bool
(BigFloat -> BigFloat -> Bool)
-> (BigFloat -> BigFloat -> Bool) -> Eq BigFloat
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BigFloat -> BigFloat -> Bool
== :: BigFloat -> BigFloat -> Bool
$c/= :: BigFloat -> BigFloat -> Bool
/= :: BigFloat -> BigFloat -> Bool
Eq, Int -> BigFloat -> ShowS
[BigFloat] -> ShowS
BigFloat -> String
(Int -> BigFloat -> ShowS)
-> (BigFloat -> String) -> ([BigFloat] -> ShowS) -> Show BigFloat
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BigFloat -> ShowS
showsPrec :: Int -> BigFloat -> ShowS
$cshow :: BigFloat -> String
show :: BigFloat -> String
$cshowList :: [BigFloat] -> ShowS
showList :: [BigFloat] -> ShowS
Show)

-- | Return the runtime number of significant binary bits retained by a BigFloat.
bigPrecision :: BigFloat -> Int
bigPrecision :: BigFloat -> Int
bigPrecision (BigFloat Int
p Rational
_) = Int
p

{- | Return the exact Rational value of the already-rounded BigFloat, without another
conversion or rounding.
-}
bigRational :: BigFloat -> Rational
bigRational :: BigFloat -> Rational
bigRational (BigFloat Int
_ Rational
x) = Rational
x

{- | Given precision in [2,1000000] bits, rounding mode, then Rational, construct
configurable-precision BigFloat plus Inexact flag if rounded. Invalid precision returns
Left InvalidArgument. Exponents are unbounded; no IEEE NaN/infinity is stored.
-}
bigFloat :: Int -> RoundingMode -> Rational -> Either MathError (Rounded BigFloat)
bigFloat :: Int
-> RoundingMode -> Rational -> Either MathError (Rounded BigFloat)
bigFloat Int
p RoundingMode
mode Rational
x
    | Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 Bool -> Bool -> Bool
|| Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
1000000 = MathError -> Either MathError (Rounded BigFloat)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"precision must be 2..1000000 bits")
    | Rational
x Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 = Rounded BigFloat -> Either MathError (Rounded BigFloat)
forall a b. b -> Either a b
Right (BigFloat -> [NumericException] -> Rounded BigFloat
forall a. a -> [NumericException] -> Rounded a
Rounded (Int -> Rational -> BigFloat
BigFloat Int
p Rational
0) [])
    | Bool
otherwise =
        let quantum :: Rational
quantum = Int -> Rational
powerTwo (Rational -> Int
floorLog2 (Rational -> Rational
forall a. Num a => a -> a
abs Rational
x) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
            y :: Rational
y = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (RoundingMode -> Rational -> Integer
roundRational RoundingMode
mode (Rational
x Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
quantum)) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
quantum
         in Rounded BigFloat -> Either MathError (Rounded BigFloat)
forall a b. b -> Either a b
Right (BigFloat -> [NumericException] -> Rounded BigFloat
forall a. a -> [NumericException] -> Rounded a
Rounded (Int -> Rational -> BigFloat
BigFloat Int
p Rational
y) [NumericException
Inexact | Rational
x Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
y])

-- | Software extended precision: 64 significant bits, unbounded exponent.
extended :: Rational -> Either MathError BigFloat
extended :: Rational -> Either MathError BigFloat
extended Rational
x = Rounded BigFloat -> BigFloat
forall a. Rounded a -> a
roundedValue (Rounded BigFloat -> BigFloat)
-> Either MathError (Rounded BigFloat) -> Either MathError BigFloat
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int
-> RoundingMode -> Rational -> Either MathError (Rounded BigFloat)
bigFloat Int
64 RoundingMode
NearestEven Rational
x

{- | Round BigFloat to host Double through fromRational. This may overflow/underflow or lose
precision; it returns no software flags or Either validation.
-}
bigToDouble :: BigFloat -> Double
bigToDouble :: BigFloat -> Double
bigToDouble = Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Rational -> Double)
-> (BigFloat -> Rational) -> BigFloat -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. BigFloat -> Rational
bigRational

bigResult :: Int -> Rational -> Either MathError BigFloat
bigResult :: Int -> Rational -> Either MathError BigFloat
bigResult Int
p Rational
x = Rounded BigFloat -> BigFloat
forall a. Rounded a -> a
roundedValue (Rounded BigFloat -> BigFloat)
-> Either MathError (Rounded BigFloat) -> Either MathError BigFloat
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int
-> RoundingMode -> Rational -> Either MathError (Rounded BigFloat)
bigFloat Int
p RoundingMode
NearestEven Rational
x
bigAdd, bigSubtract, bigMultiply, bigDivide :: BigFloat -> BigFloat -> Either MathError BigFloat

{- | Add two BigFloats exactly then round once to the larger operand precision,
nearest-even. Return the new BigFloat in Either.
-}
bigAdd :: BigFloat -> BigFloat -> Either MathError BigFloat
bigAdd BigFloat
a BigFloat
b = Int -> Rational -> Either MathError BigFloat
bigResult (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (BigFloat -> Int
bigPrecision BigFloat
a) (BigFloat -> Int
bigPrecision BigFloat
b)) (BigFloat -> Rational
bigRational BigFloat
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ BigFloat -> Rational
bigRational BigFloat
b)

{- | Subtract the second BigFloat from the first, rounding once nearest-even at the larger
operand precision.
-}
bigSubtract :: BigFloat -> BigFloat -> Either MathError BigFloat
bigSubtract BigFloat
a BigFloat
b = Int -> Rational -> Either MathError BigFloat
bigResult (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (BigFloat -> Int
bigPrecision BigFloat
a) (BigFloat -> Int
bigPrecision BigFloat
b)) (BigFloat -> Rational
bigRational BigFloat
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- BigFloat -> Rational
bigRational BigFloat
b)

{- | Multiply two BigFloats exactly then round once nearest-even to the larger operand
precision.
-}
bigMultiply :: BigFloat -> BigFloat -> Either MathError BigFloat
bigMultiply BigFloat
a BigFloat
b = Int -> Rational -> Either MathError BigFloat
bigResult (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (BigFloat -> Int
bigPrecision BigFloat
a) (BigFloat -> Int
bigPrecision BigFloat
b)) (BigFloat -> Rational
bigRational BigFloat
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* BigFloat -> Rational
bigRational BigFloat
b)

{- | Divide two BigFloats and round once nearest-even to the larger precision. Zero divisor
returns Left DomainError.
-}
bigDivide :: BigFloat -> BigFloat -> Either MathError BigFloat
bigDivide BigFloat
_ (BigFloat Int
_ Rational
0) = MathError -> Either MathError BigFloat
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BigFloat division by zero")
bigDivide BigFloat
a BigFloat
b = Int -> Rational -> Either MathError BigFloat
bigResult (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (BigFloat -> Int
bigPrecision BigFloat
a) (BigFloat -> Int
bigPrecision BigFloat
b)) (BigFloat -> Rational
bigRational BigFloat
a Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ BigFloat -> Rational
bigRational BigFloat
b)

mixedAdd, mixedMultiply :: BigFloat -> Double -> Either MathError BigFloat

{- | Given BigFloat then finite Double, add via exact Rational intermediates at max(BigFloat
precision,53) bits, nearest-even. Non-finite Double fails.
-}
mixedAdd :: BigFloat -> Double -> Either MathError BigFloat
mixedAdd BigFloat
a Double
x
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x = MathError -> Either MathError BigFloat
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"non-finite mixed operand")
    | Bool
otherwise = Int -> Rational -> Either MathError BigFloat
bigResult (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
53 (BigFloat -> Int
bigPrecision BigFloat
a)) (BigFloat -> Rational
bigRational BigFloat
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x)

{- | Given BigFloat then finite Double, multiply at max(BigFloat precision,53) bits using
exact intermediates and nearest-even. Non-finite Double fails.
-}
mixedMultiply :: BigFloat -> Double -> Either MathError BigFloat
mixedMultiply BigFloat
a Double
x
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x = MathError -> Either MathError BigFloat
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"non-finite mixed operand")
    | Bool
otherwise = Int -> Rational -> Either MathError BigFloat
bigResult (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
53 (BigFloat -> Int
bigPrecision BigFloat
a)) (BigFloat -> Rational
bigRational BigFloat
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x)

integerSqrt :: Integer -> Integer
integerSqrt :: Integer -> Integer
integerSqrt Integer
n
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
2 = Integer
n
    | Bool
otherwise = Integer -> Integer
go (Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ ((Integer -> Int
integerBits Integer
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2))
  where
    go :: Integer -> Integer
go Integer
x = let y :: Integer
y = (Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
x) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
2 in if Integer
y Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
x then Integer
x else Integer -> Integer
go Integer
y

{- | Compute nonnegative BigFloat square root using integer/rational guard precision then
round to the input precision. Negative input returns Left DomainError; zero succeeds.
-}
bigSqrt :: BigFloat -> Either MathError BigFloat
bigSqrt :: BigFloat -> Either MathError BigFloat
bigSqrt (BigFloat Int
p Rational
x)
    | Rational
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
0 = MathError -> Either MathError BigFloat
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"negative BigFloat square root")
    | Rational
x Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0 = Int -> Rational -> Either MathError BigFloat
bigResult Int
p Rational
0
    | Bool
otherwise =
        let e :: Int
e = Rational -> Int
floorLog2 Rational
x Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
8
            q :: Rational
q = Int -> Rational
powerTwo Int
e
            n :: Integer
n = Integer -> Integer
integerSqrt (Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Rational
x Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ (Rational
q Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
q)))
         in Int -> Rational -> Either MathError BigFloat
bigResult Int
p (Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
n Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
q)

{- | Absolute tail controlled rational series; final rounding is accurate to
the requested precision up to one final ulp, not a correctly-rounded promise.
-}
bigExp :: BigFloat -> Either MathError BigFloat
bigExp :: BigFloat -> Either MathError BigFloat
bigExp (BigFloat Int
p Rational
x)
    | Rational -> Rational
forall a. Num a => a -> a
abs Rational
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
100000 =
        MathError -> Either MathError BigFloat
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"BigFloat exponential argument exceeds resource limit 100000")
    | Rational
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
0 = do y <- BigFloat -> Either MathError BigFloat
bigExp (Int -> Rational -> BigFloat
BigFloat (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
16) (-Rational
x)); bigResult p (1 / bigRational y)
    | Bool
otherwise =
        let reductions :: Int
reductions = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Rational -> Int
floorLog2 (Rational -> Rational -> Rational
forall a. Ord a => a -> a -> a
max Rational
1 Rational
x) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
            y :: Rational
y = Rational
x Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Int -> Rational
powerTwo Int
reductions
            eps :: Rational
eps = Int -> Rational
powerTwo (Int -> Int
forall a. Num a => a -> a
negate (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
reductions Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
24))
            series :: t -> Rational -> Rational -> Rational
series t
k Rational
term Rational
acc =
                let term' :: Rational
term' = Rational
term Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ t -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral t
k; acc' :: Rational
acc' = Rational
acc Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
term'
                 in if Rational -> Rational
forall a. Num a => a -> a
abs Rational
term' Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
eps then Rational
acc' else t -> Rational -> Rational -> Rational
series (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Rational
term' Rational
acc'
            initial :: Rational
initial = Integer -> Rational -> Rational -> Rational
forall {t}. Integral t => t -> Rational -> Rational -> Rational
series (Integer
1 :: Integer) Rational
1 Rational
1
            sq :: t -> Rational -> Either MathError Rational
sq t
0 Rational
v = Rational -> Either MathError Rational
forall a b. b -> Either a b
Right Rational
v
            sq t
k Rational
v = do rounded <- Int -> Rational -> Either MathError BigFloat
bigResult (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
reductions Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
24) (Rational
v Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
v); sq (k - 1) (bigRational rounded)
         in Int -> Rational -> Either MathError Rational
forall {t}.
(Eq t, Num t) =>
t -> Rational -> Either MathError Rational
sq Int
reductions Rational
initial Either MathError Rational
-> (Rational -> Either MathError BigFloat)
-> Either MathError BigFloat
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Int -> Rational -> Either MathError BigFloat
bigResult Int
p

{- | Compute natural logarithm at the BigFloat precision using rational series and guard
bits. Require positive input; invalid domain/resource limits return Left MathError.
Final rounding has the module's documented accuracy limits.
-}
bigLog :: BigFloat -> Either MathError BigFloat
bigLog :: BigFloat -> Either MathError BigFloat
bigLog (BigFloat Int
p Rational
x)
    | Rational
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
0 = MathError -> Either MathError BigFloat
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BigFloat logarithm needs positive input")
    | Bool
otherwise =
        let e :: Int
e = Rational -> Int
floorLog2 Rational
x
            mant :: Rational
mant = Rational
x Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Int -> Rational
powerTwo Int
e
            eps :: Rational
eps = Int -> Rational
powerTwo (Int -> Int
forall a. Num a => a -> a
negate (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerBits (Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Int
forall a. Num a => a -> a
abs Int
e) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
24))
            atanhSeries :: Rational -> Rational
atanhSeries Rational
y = Integer -> Rational -> Rational -> Rational
forall {t}. Integral t => t -> Rational -> Rational -> Rational
go (Integer
1 :: Integer) Rational
y Rational
y
              where
                go :: t -> Rational -> Rational -> Rational
go t
k Rational
term Rational
acc =
                    let term' :: Rational
term' = Rational
term Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
y; add :: Rational
add = Rational
term' Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ t -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral (t
2 t -> t -> t
forall a. Num a => a -> a -> a
* t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
                     in if Rational -> Rational
forall a. Num a => a -> a
abs Rational
add Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
eps then Rational
2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* (Rational
acc Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
add) else t -> Rational -> Rational -> Rational
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
+ t
1) Rational
term' (Rational
acc Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
add)
            l2 :: Rational
l2 = Rational -> Rational
atanhSeries (Rational
1 Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
3)
         in Int -> Rational -> Either MathError BigFloat
bigResult Int
p (Rational -> Rational
atanhSeries ((Rational
mant Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
1) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ (Rational
mant Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Rational
1)) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Int -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
e Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
l2)

{- | Runtime-scale base-ten fixed point: coefficient * 10^(-scale). All
constructors and arithmetic use explicit rounding and report inexactness.
-}
data Fixed = Fixed Int Integer deriving (Fixed -> Fixed -> Bool
(Fixed -> Fixed -> Bool) -> (Fixed -> Fixed -> Bool) -> Eq Fixed
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Fixed -> Fixed -> Bool
== :: Fixed -> Fixed -> Bool
$c/= :: Fixed -> Fixed -> Bool
/= :: Fixed -> Fixed -> Bool
Eq, Int -> Fixed -> ShowS
[Fixed] -> ShowS
Fixed -> String
(Int -> Fixed -> ShowS)
-> (Fixed -> String) -> ([Fixed] -> ShowS) -> Show Fixed
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Fixed -> ShowS
showsPrec :: Int -> Fixed -> ShowS
$cshow :: Fixed -> String
show :: Fixed -> String
$cshowList :: [Fixed] -> ShowS
showList :: [Fixed] -> ShowS
Show)

{- | Return nonnegative decimal scale s of a Fixed value; its exact value is coefficient /
10^s.
-}
fixedScale :: Fixed -> Int
fixedScale :: Fixed -> Int
fixedScale (Fixed Int
s Integer
_) = Int
s

{- | Return exact stored Integer coefficient of a Fixed value, before division by
10^fixedScale.
-}
fixedCoefficient :: Fixed -> Integer
fixedCoefficient :: Fixed -> Integer
fixedCoefficient (Fixed Int
_ Integer
c) = Integer
c

-- | Return exact coefficient/10^scale represented by a Fixed value, without rounding.
fixedRational :: Fixed -> Rational
fixedRational :: Fixed -> Rational
fixedRational (Fixed Int
s Integer
n) = Integer
n Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Integer
10 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
s)

{- | Given decimal scale in [0,10000], rounding mode, then Rational, return Rounded Fixed
and Inexact when required. Invalid scale fails. This scale counts decimal fractional
places, not significant digits.
-}
fixed :: Int -> RoundingMode -> Rational -> Either MathError (Rounded Fixed)
fixed :: Int -> RoundingMode -> Rational -> Either MathError (Rounded Fixed)
fixed Int
s RoundingMode
mode Rational
x
    | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
10000 = MathError -> Either MathError (Rounded Fixed)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"fixed scale must be 0..10000")
    | Bool
otherwise =
        let c :: Integer
c = RoundingMode -> Rational -> Integer
roundRational RoundingMode
mode (Rational
x Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer
10 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
s)); v :: Fixed
v = Int -> Integer -> Fixed
Fixed Int
s Integer
c
         in Rounded Fixed -> Either MathError (Rounded Fixed)
forall a b. b -> Either a b
Right (Fixed -> [NumericException] -> Rounded Fixed
forall a. a -> [NumericException] -> Rounded a
Rounded Fixed
v [NumericException
Inexact | Fixed -> Rational
fixedRational Fixed
v Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
/= Rational
x])

{- | Add two Fixed values using the larger decimal scale, nearest-even, returning Rounded
result and software flags.
-}
fixedAdd :: Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedAdd :: Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedAdd Fixed
a Fixed
b = Int -> RoundingMode -> Rational -> Either MathError (Rounded Fixed)
fixed (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Fixed -> Int
fixedScale Fixed
a) (Fixed -> Int
fixedScale Fixed
b)) RoundingMode
NearestEven (Fixed -> Rational
fixedRational Fixed
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Fixed -> Rational
fixedRational Fixed
b)

{- | Given rounding mode then two Fixed values, multiply and round to their larger scale,
returning Rounded result and flags.
-}
fixedMultiply :: RoundingMode -> Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedMultiply :: RoundingMode -> Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedMultiply RoundingMode
mode Fixed
a Fixed
b = Int -> RoundingMode -> Rational -> Either MathError (Rounded Fixed)
fixed (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Fixed -> Int
fixedScale Fixed
a) (Fixed -> Int
fixedScale Fixed
b)) RoundingMode
mode (Fixed -> Rational
fixedRational Fixed
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Fixed -> Rational
fixedRational Fixed
b)

{- | Given rounding mode then Fixed numerator and denominator, divide and round to their
larger scale. Zero divisor returns Left DomainError; successful rounding flags are
retained.
-}
fixedDivide :: RoundingMode -> Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedDivide :: RoundingMode -> Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedDivide RoundingMode
mode Fixed
a Fixed
b
    | Fixed -> Integer
fixedCoefficient Fixed
b Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = MathError -> Either MathError (Rounded Fixed)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"fixed division by zero")
    | Bool
otherwise = Int -> RoundingMode -> Rational -> Either MathError (Rounded Fixed)
fixed (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Fixed -> Int
fixedScale Fixed
a) (Fixed -> Int
fixedScale Fixed
b)) RoundingMode
mode (Fixed -> Rational
fixedRational Fixed
a Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Fixed -> Rational
fixedRational Fixed
b)

{- | Opaque decimal fixed-point value with fractional-place scale, built through decimal
with explicit rounding. This is not an implicit significant-digit decimal context.
-}
newtype Decimal = Decimal Fixed deriving (Decimal -> Decimal -> Bool
(Decimal -> Decimal -> Bool)
-> (Decimal -> Decimal -> Bool) -> Eq Decimal
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Decimal -> Decimal -> Bool
== :: Decimal -> Decimal -> Bool
$c/= :: Decimal -> Decimal -> Bool
/= :: Decimal -> Decimal -> Bool
Eq, Int -> Decimal -> ShowS
[Decimal] -> ShowS
Decimal -> String
(Int -> Decimal -> ShowS)
-> (Decimal -> String) -> ([Decimal] -> ShowS) -> Show Decimal
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Decimal -> ShowS
showsPrec :: Int -> Decimal -> ShowS
$cshow :: Decimal -> String
show :: Decimal -> String
$cshowList :: [Decimal] -> ShowS
showList :: [Decimal] -> ShowS
Show)

{- | Given decimal fractional-place scale, rounding mode, then Rational, construct Rounded
Decimal through fixed. Scale must lie in [0,10000]; no significant-digit interpretation
is implied.
-}
decimal :: Int -> RoundingMode -> Rational -> Either MathError (Rounded Decimal)
decimal :: Int
-> RoundingMode -> Rational -> Either MathError (Rounded Decimal)
decimal Int
p RoundingMode
mode Rational
x = do Rounded v flags <- Int -> RoundingMode -> Rational -> Either MathError (Rounded Fixed)
fixed Int
p RoundingMode
mode Rational
x; Right (Rounded (Decimal v) flags)

-- | Return exact Rational value stored in a Decimal without further rounding.
decimalRational :: Decimal -> Rational
decimalRational :: Decimal -> Rational
decimalRational (Decimal Fixed
f) = Fixed -> Rational
fixedRational Fixed
f

{- | Add two Decimals at their larger decimal scale, nearest-even, returning Rounded result
and software flags.
-}
decimalAdd :: Decimal -> Decimal -> Either MathError (Rounded Decimal)
decimalAdd :: Decimal -> Decimal -> Either MathError (Rounded Decimal)
decimalAdd (Decimal Fixed
a) (Decimal Fixed
b) = do Rounded v flags <- Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedAdd Fixed
a Fixed
b; Right (Rounded (Decimal v) flags)

decimalMultiply
    , decimalDivide ::
        RoundingMode -> Decimal -> Decimal -> Either MathError (Rounded Decimal)

{- | Given rounding mode then two Decimals, multiply at the larger fractional-place scale
and return Rounded result with flags.
-}
decimalMultiply :: RoundingMode
-> Decimal -> Decimal -> Either MathError (Rounded Decimal)
decimalMultiply RoundingMode
mode (Decimal Fixed
a) (Decimal Fixed
b) = do Rounded v flags <- RoundingMode -> Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedMultiply RoundingMode
mode Fixed
a Fixed
b; Right (Rounded (Decimal v) flags)

{- | Given rounding mode then two Decimals, divide at their larger scale and return Rounded
result. Zero divisor returns Left DomainError.
-}
decimalDivide :: RoundingMode
-> Decimal -> Decimal -> Either MathError (Rounded Decimal)
decimalDivide RoundingMode
mode (Decimal Fixed
a) (Decimal Fixed
b) = do Rounded v flags <- RoundingMode -> Fixed -> Fixed -> Either MathError (Rounded Fixed)
fixedDivide RoundingMode
mode Fixed
a Fixed
b; Right (Rounded (Decimal v) flags)

{- | Finite closed binary64 intervals. Arithmetic endpoints are rounded outward
using exact rationals, so results enclose exact real arithmetic. Unbounded
output is rejected; division rejects denominators containing zero.
-}
data Interval = Interval Double Double deriving (Interval -> Interval -> Bool
(Interval -> Interval -> Bool)
-> (Interval -> Interval -> Bool) -> Eq Interval
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Interval -> Interval -> Bool
== :: Interval -> Interval -> Bool
$c/= :: Interval -> Interval -> Bool
/= :: Interval -> Interval -> Bool
Eq, Int -> Interval -> ShowS
[Interval] -> ShowS
Interval -> String
(Int -> Interval -> ShowS)
-> (Interval -> String) -> ([Interval] -> ShowS) -> Show Interval
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Interval -> ShowS
showsPrec :: Int -> Interval -> ShowS
$cshow :: Interval -> String
show :: Interval -> String
$cshowList :: [Interval] -> ShowS
showList :: [Interval] -> ShowS
Show)

{- | Given finite lower and upper bounds in increasing order, construct a closed binary64
Interval. Reversed or non-finite bounds fail. Equal endpoints are valid.
-}
interval :: Double -> Double -> Either MathError Interval
interval :: Double -> Double -> Either MathError Interval
interval Double
a Double
b
    | (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\Double
x -> Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x) [Double
a, Double
b] Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
b = MathError -> Either MathError Interval
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"invalid finite interval")
    | Bool
otherwise = Interval -> Either MathError Interval
forall a b. b -> Either a b
Right (Double -> Double -> Interval
Interval Double
a Double
b)

-- | Return (lower,upper) binary64 endpoints of a validated closed interval.
intervalBounds :: Interval -> (Double, Double)
intervalBounds :: Interval -> (Double, Double)
intervalBounds (Interval Double
a Double
b) = (Double
a, Double
b)

{- | Given interval then Double, test closed-bound membership lower\<=x\<=upper. NaN
comparisons give False.
-}
intervalContains :: Interval -> Double -> Bool
intervalContains :: Interval -> Double -> Bool
intervalContains (Interval Double
a Double
b) Double
x = Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
b

outward :: Rational -> Rational -> Either MathError Interval
outward :: Rational -> Rational -> Either MathError Interval
outward Rational
lo Rational
hi = Double -> Double -> Either MathError Interval
interval (Rational -> Double
forall {a}. RealFloat a => Rational -> a
down Rational
lo) (Rational -> Double
forall {a}. RealFloat a => Rational -> a
up Rational
hi)
  where
    down :: Rational -> a
down Rational
q =
        let x :: a
x = Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
q
         in if a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x then a
x else if a -> Rational
forall a. Real a => a -> Rational
toRational a
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
q then a -> a -> a
forall a. RealFloat a => a -> a -> a
nextAfter a
x (-a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0) else a
x
    up :: Rational -> a
up Rational
q =
        let x :: a
x = Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
q
         in if a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x then a
x else if a -> Rational
forall a. Real a => a -> Rational
toRational a
x Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
q then a -> a -> a
forall a. RealFloat a => a -> a -> a
nextAfter a
x (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0) else a
x
intervalAdd
    , intervalSubtract
    , intervalMultiply
    , intervalDivide ::
        Interval -> Interval -> Either MathError Interval

{- | Add finite intervals with exact Rational endpoint arithmetic rounded outward, enclosing
the exact sum. Results requiring unbounded endpoints fail.
-}
intervalAdd :: Interval -> Interval -> Either MathError Interval
intervalAdd (Interval Double
a Double
b) (Interval Double
c Double
d) = Rational -> Rational -> Either MathError Interval
outward (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Double -> Rational
forall a. Real a => a -> Rational
toRational Double
c) (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Double -> Rational
forall a. Real a => a -> Rational
toRational Double
d)

{- | Subtract intervals with endpoints [a-d,b-c], rounding outward. Non-finite/unbounded
resulting endpoints return Left MathError.
-}
intervalSubtract :: Interval -> Interval -> Either MathError Interval
intervalSubtract (Interval Double
a Double
b) (Interval Double
c Double
d) = Rational -> Rational -> Either MathError Interval
outward (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
d) (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Double -> Rational
forall a. Real a => a -> Rational
toRational Double
c)

{- | Multiply intervals using the minimum/maximum of all four endpoint products and outward
rounding. Reject unbounded output.
-}
intervalMultiply :: Interval -> Interval -> Either MathError Interval
intervalMultiply (Interval Double
a Double
b) (Interval Double
c Double
d) =
    let xs :: [Rational]
xs = [Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Double -> Rational
forall a. Real a => a -> Rational
toRational Double
y | Double
x <- [Double
a, Double
b], Double
y <- [Double
c, Double
d]]
     in Rational -> Rational -> Either MathError Interval
outward ([Rational] -> Rational
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Rational]
xs) ([Rational] -> Rational
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Rational]
xs)

{- | Divide intervals using endpoint quotients with outward rounding. A divisor interval
containing zero returns Left DomainError; unbounded output also fails.
-}
intervalDivide :: Interval -> Interval -> Either MathError Interval
intervalDivide (Interval Double
a Double
b) (Interval Double
c Double
d)
    | Double
c Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
&& Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 = MathError -> Either MathError Interval
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"interval divisor contains zero")
    | Bool
otherwise =
        let xs :: [Rational]
xs = [Double -> Rational
forall a. Real a => a -> Rational
toRational Double
x Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Double -> Rational
forall a. Real a => a -> Rational
toRational Double
y | Double
x <- [Double
a, Double
b], Double
y <- [Double
c, Double
d]]
         in Rational -> Rational -> Either MathError Interval
outward ([Rational] -> Rational
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Rational]
xs) ([Rational] -> Rational
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Rational]
xs)

{- | Forward-mode automatic differentiation, f(a+b epsilon), epsilon^2=0.
Elementary functions follow host IEEE domains; checked algorithms can inspect
both components. abs at zero uses derivative zero (a chosen subgradient).
-}
data Dual a = Dual
    { forall a. Dual a -> a
primal :: a
    -- ^ Primal component of a forward-mode Dual number.
    , forall a. Dual a -> a
tangent :: a
    -- ^ Directional derivative component of a forward-mode Dual number.
    }
    deriving (Dual a -> Dual a -> Bool
(Dual a -> Dual a -> Bool)
-> (Dual a -> Dual a -> Bool) -> Eq (Dual a)
forall a. Eq a => Dual a -> Dual a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Dual a -> Dual a -> Bool
== :: Dual a -> Dual a -> Bool
$c/= :: forall a. Eq a => Dual a -> Dual a -> Bool
/= :: Dual a -> Dual a -> Bool
Eq, Int -> Dual a -> ShowS
[Dual a] -> ShowS
Dual a -> String
(Int -> Dual a -> ShowS)
-> (Dual a -> String) -> ([Dual a] -> ShowS) -> Show (Dual a)
forall a. Show a => Int -> Dual a -> ShowS
forall a. Show a => [Dual a] -> ShowS
forall a. Show a => Dual a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Dual a -> ShowS
showsPrec :: Int -> Dual a -> ShowS
$cshow :: forall a. Show a => Dual a -> String
show :: Dual a -> String
$cshowList :: forall a. Show a => [Dual a] -> ShowS
showList :: [Dual a] -> ShowS
Show)

instance (Num a) => Num (Dual a) where
    Dual a
a a
b + :: Dual a -> Dual a -> Dual a
+ Dual a
c a
d = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
c) (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
d)
    Dual a
a a
b - :: Dual a -> Dual a -> Dual a
- Dual a
c a
d = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
c) (a
b a -> a -> a
forall a. Num a => a -> a -> a
- a
d)
    Dual a
a a
b * :: Dual a -> Dual a -> Dual a
* Dual a
c a
d = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
c) (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
d a -> a -> a
forall a. Num a => a -> a -> a
+ a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
c)
    negate :: Dual a -> Dual a
negate (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Num a => a -> a
negate a
a) (a -> a
forall a. Num a => a -> a
negate a
b)
    abs :: Dual a -> Dual a
abs (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Num a => a -> a
abs a
a) (a -> a
forall a. Num a => a -> a
signum a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b)
    signum :: Dual a -> Dual a
signum (Dual a
a a
_) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Num a => a -> a
signum a
a) a
0
    fromInteger :: Integer -> Dual a
fromInteger Integer
n = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
n) a
0
instance (Fractional a) => Fractional (Dual a) where
    recip :: Dual a -> Dual a
recip (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Fractional a => a -> a
recip a
a) (-a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a))
    fromRational :: Rational -> Dual a
fromRational Rational
r = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
r) a
0
instance (Floating a) => Floating (Dual a) where
    pi :: Dual a
pi = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual a
forall a. Floating a => a
pi a
0
    exp :: Dual a -> Dual a
exp (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
exp a
a) (a
b a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
exp a
a)
    log :: Dual a -> Dual a
log (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
log a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
a)
    sqrt :: Dual a -> Dual a
sqrt (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
sqrt a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sqrt a
a))
    sin :: Dual a -> Dual a
sin (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
sin a
a) (a
b a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos a
a)
    cos :: Dual a -> Dual a
cos (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
cos a
a) (-a
b a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sin a
a)
    tan :: Dual a -> Dual a
tan (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
tan a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a -> a
forall a. Floating a => a -> a
cos a
a a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos a
a))
    asin :: Dual a -> Dual a
asin (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
asin a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a))
    acos :: Dual a -> Dual a
acos (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
acos a
a) (-a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a))
    atan :: Dual a -> Dual a
atan (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
atan a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a))
    sinh :: Dual a -> Dual a
sinh (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
sinh a
a) (a
b a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cosh a
a)
    cosh :: Dual a -> Dual a
cosh (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
cosh a
a) (a
b a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sinh a
a)
    tanh :: Dual a -> Dual a
tanh (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
tanh a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a -> a
forall a. Floating a => a -> a
cosh a
a a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cosh a
a))
    asinh :: Dual a -> Dual a
asinh (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
asinh a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a))
    acosh :: Dual a -> Dual a
acosh (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
acosh a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
1))
    atanh :: Dual a -> Dual a
atanh (Dual a
a a
b) = a -> a -> Dual a
forall a. a -> a -> Dual a
Dual (a -> a
forall a. Floating a => a -> a
atanh a
a) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
a))

{- | Given a function on Dual numbers then primal input x, seed Dual x 1 and extract the
output tangent. Arithmetic/domains follow Dual's Num/Floating instances; no Either
validation is added.
-}
derivative :: (Num a) => (Dual a -> Dual a) -> a -> a
derivative :: forall a. Num a => (Dual a -> Dual a) -> a -> a
derivative Dual a -> Dual a
f a
x = Dual a -> a
forall a. Dual a -> a
tangent (Dual a -> Dual a
f (a -> a -> Dual a
forall a. a -> a -> Dual a
Dual a
x a
1))

-- | Hamilton quaternion (real, i, j, k); multiplication is noncommutative.
data Quaternion a = Quaternion a a a a deriving (Quaternion a -> Quaternion a -> Bool
(Quaternion a -> Quaternion a -> Bool)
-> (Quaternion a -> Quaternion a -> Bool) -> Eq (Quaternion a)
forall a. Eq a => Quaternion a -> Quaternion a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Quaternion a -> Quaternion a -> Bool
== :: Quaternion a -> Quaternion a -> Bool
$c/= :: forall a. Eq a => Quaternion a -> Quaternion a -> Bool
/= :: Quaternion a -> Quaternion a -> Bool
Eq, Int -> Quaternion a -> ShowS
[Quaternion a] -> ShowS
Quaternion a -> String
(Int -> Quaternion a -> ShowS)
-> (Quaternion a -> String)
-> ([Quaternion a] -> ShowS)
-> Show (Quaternion a)
forall a. Show a => Int -> Quaternion a -> ShowS
forall a. Show a => [Quaternion a] -> ShowS
forall a. Show a => Quaternion a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Quaternion a -> ShowS
showsPrec :: Int -> Quaternion a -> ShowS
$cshow :: forall a. Show a => Quaternion a -> String
show :: Quaternion a -> String
$cshowList :: forall a. Show a => [Quaternion a] -> ShowS
showList :: [Quaternion a] -> ShowS
Show)

instance (RealFloat a) => Num (Quaternion a) where
    Quaternion a
a a
b a
c a
d + :: Quaternion a -> Quaternion a -> Quaternion a
+ Quaternion a
e a
f a
g a
h = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
e) (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
f) (a
c a -> a -> a
forall a. Num a => a -> a -> a
+ a
g) (a
d a -> a -> a
forall a. Num a => a -> a -> a
+ a
h)
    Quaternion a
a a
b a
c a
d - :: Quaternion a -> Quaternion a -> Quaternion a
- Quaternion a
e a
f a
g a
h = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
e) (a
b a -> a -> a
forall a. Num a => a -> a -> a
- a
f) (a
c a -> a -> a
forall a. Num a => a -> a -> a
- a
g) (a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
h)
    Quaternion a
a a
b a
c a
d * :: Quaternion a -> Quaternion a -> Quaternion a
* Quaternion a
e a
f a
g a
h =
        a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion
            (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
e a -> a -> a
forall a. Num a => a -> a -> a
- a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
f a -> a -> a
forall a. Num a => a -> a -> a
- a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
g a -> a -> a
forall a. Num a => a -> a -> a
- a
d a -> a -> a
forall a. Num a => a -> a -> a
* a
h)
            (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
f a -> a -> a
forall a. Num a => a -> a -> a
+ a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
e a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
h a -> a -> a
forall a. Num a => a -> a -> a
- a
d a -> a -> a
forall a. Num a => a -> a -> a
* a
g)
            (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
g a -> a -> a
forall a. Num a => a -> a -> a
- a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
h a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
e a -> a -> a
forall a. Num a => a -> a -> a
+ a
d a -> a -> a
forall a. Num a => a -> a -> a
* a
f)
            (a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
h a -> a -> a
forall a. Num a => a -> a -> a
+ a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
g a -> a -> a
forall a. Num a => a -> a -> a
- a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
f a -> a -> a
forall a. Num a => a -> a -> a
+ a
d a -> a -> a
forall a. Num a => a -> a -> a
* a
e)
    negate :: Quaternion a -> Quaternion a
negate (Quaternion a
a a
b a
c a
d) = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (-a
a) (-a
b) (-a
c) (-a
d)
    fromInteger :: Integer -> Quaternion a
fromInteger Integer
n = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
n) a
0 a
0 a
0
    abs :: Quaternion a -> Quaternion a
abs Quaternion a
q = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (Quaternion a -> a
forall a. RealFloat a => Quaternion a -> a
quaternionNorm Quaternion a
q) a
0 a
0 a
0
    signum :: Quaternion a -> Quaternion a
signum Quaternion a
q = let n :: a
n = Quaternion a -> a
forall a. RealFloat a => Quaternion a -> a
quaternionNorm Quaternion a
q in if a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion a
0 a
0 a
0 a
0 else a -> Quaternion a -> Quaternion a
forall a. Num a => a -> Quaternion a -> Quaternion a
qScale (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
n) Quaternion a
q

{- | For Quaternion real i j k, negate the three imaginary components and preserve the real
part.
-}
quaternionConjugate :: (Num a) => Quaternion a -> Quaternion a
quaternionConjugate :: forall a. Num a => Quaternion a -> Quaternion a
quaternionConjugate (Quaternion a
a a
b a
c a
d) = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion a
a (-a
b) (-a
c) (-a
d)

{- | Return scaled Euclidean norm of four quaternion components, reducing avoidable
intermediate overflow. Empty/zero-component quaternion has norm zero.
-}
quaternionNorm :: (RealFloat a) => Quaternion a -> a
quaternionNorm :: forall a. RealFloat a => Quaternion a -> a
quaternionNorm (Quaternion a
a a
b a
c a
d)
    | a
scale a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a
0
    | Bool
otherwise = a
scale a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sqrt ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [let v :: a
v = a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale in a
v a -> a -> a
forall a. Num a => a -> a -> a
* a
v | a
x <- [a
a, a
b, a
c, a
d]])
  where
    scale :: a
scale = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs [a
a, a
b, a
c, a
d])

qScale :: (Num a) => a -> Quaternion a -> Quaternion a
qScale :: forall a. Num a => a -> Quaternion a -> Quaternion a
qScale a
s (Quaternion a
a a
b a
c a
d) = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
a) (a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
b) (a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
c) (a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
d)
qDivide :: (Fractional a) => Quaternion a -> a -> Quaternion a
qDivide :: forall a. Fractional a => Quaternion a -> a -> Quaternion a
qDivide (Quaternion a
a a
b a
c a
d) a
s = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion (a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) (a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) (a
c a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) (a
d a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s)

{- | Invert a finite nonzero quaternion using scaled norm arithmetic. Zero/non-finite input
returns DomainError; non-finite inverse returns NonFiniteResult.
-}
quaternionInverse :: (RealFloat a) => Quaternion a -> Either MathError (Quaternion a)
quaternionInverse :: forall a.
RealFloat a =>
Quaternion a -> Either MathError (Quaternion a)
quaternionInverse q :: Quaternion a
q@(Quaternion a
a a
b a
c a
d)
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\a
x -> a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x) [a
a, a
b, a
c, a
d] Bool -> Bool -> Bool
|| [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs [a
a, a
b, a
c, a
d]) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 =
        MathError -> Either MathError (Quaternion a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"invalid quaternion inverse")
    | Bool
otherwise =
        let scale :: a
scale = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs [a
a, a
b, a
c, a
d])
            qs :: Quaternion a
qs = Quaternion a -> a -> Quaternion a
forall a. Fractional a => Quaternion a -> a -> Quaternion a
qDivide Quaternion a
q a
scale
            Quaternion a
e a
f a
g a
h = Quaternion a
qs
            denom :: a
denom = a
e a -> a -> a
forall a. Num a => a -> a -> a
* a
e a -> a -> a
forall a. Num a => a -> a -> a
+ a
f a -> a -> a
forall a. Num a => a -> a -> a
* a
f a -> a -> a
forall a. Num a => a -> a -> a
+ a
g a -> a -> a
forall a. Num a => a -> a -> a
* a
g a -> a -> a
forall a. Num a => a -> a -> a
+ a
h a -> a -> a
forall a. Num a => a -> a -> a
* a
h
            out :: Quaternion a
out = Quaternion a -> a -> Quaternion a
forall a. Fractional a => Quaternion a -> a -> Quaternion a
qDivide (Quaternion a -> a -> Quaternion a
forall a. Fractional a => Quaternion a -> a -> Quaternion a
qDivide (Quaternion a -> Quaternion a
forall a. Num a => Quaternion a -> Quaternion a
quaternionConjugate Quaternion a
qs) a
denom) a
scale
         in if Quaternion a -> Bool
forall a. RealFloat a => Quaternion a -> Bool
qFinite Quaternion a
out then Quaternion a -> Either MathError (Quaternion a)
forall a b. b -> Either a b
Right Quaternion a
out else MathError -> Either MathError (Quaternion a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult

qFinite :: (RealFloat a) => Quaternion a -> Bool
qFinite :: forall a. RealFloat a => Quaternion a -> Bool
qFinite (Quaternion a
a a
b a
c a
d) = (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\a
x -> Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x)) [a
a, a
b, a
c, a
d]

{- | Given nonzero finite quaternion q then (x,y,z), return imaginary part of
q*(0,x,y,z)*inverse(q). Normalization of q is implicit in the inverse; invalid
input/result returns Left MathError.
-}
rotateVector :: (RealFloat a) => Quaternion a -> (a, a, a) -> Either MathError (a, a, a)
rotateVector :: forall a.
RealFloat a =>
Quaternion a -> (a, a, a) -> Either MathError (a, a, a)
rotateVector Quaternion a
q (a
x, a
y, a
z) = do
    qi <- Quaternion a -> Either MathError (Quaternion a)
forall a.
RealFloat a =>
Quaternion a -> Either MathError (Quaternion a)
quaternionInverse Quaternion a
q
    let Quaternion _ a b c = q * Quaternion 0 x y z * qi
    if all (\a
v -> Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
v Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
v)) [a, b, c]
        then Right (a, b, c)
        else Left NonFiniteResult

{- | Real and dual quaternion pair, commonly representing a rigid rotation/translation.
Construct normalized rigid transforms with rigidTransform.
-}
data DualQuaternion a = DualQuaternion (Quaternion a) (Quaternion a) deriving (DualQuaternion a -> DualQuaternion a -> Bool
(DualQuaternion a -> DualQuaternion a -> Bool)
-> (DualQuaternion a -> DualQuaternion a -> Bool)
-> Eq (DualQuaternion a)
forall a. Eq a => DualQuaternion a -> DualQuaternion a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => DualQuaternion a -> DualQuaternion a -> Bool
== :: DualQuaternion a -> DualQuaternion a -> Bool
$c/= :: forall a. Eq a => DualQuaternion a -> DualQuaternion a -> Bool
/= :: DualQuaternion a -> DualQuaternion a -> Bool
Eq, Int -> DualQuaternion a -> ShowS
[DualQuaternion a] -> ShowS
DualQuaternion a -> String
(Int -> DualQuaternion a -> ShowS)
-> (DualQuaternion a -> String)
-> ([DualQuaternion a] -> ShowS)
-> Show (DualQuaternion a)
forall a. Show a => Int -> DualQuaternion a -> ShowS
forall a. Show a => [DualQuaternion a] -> ShowS
forall a. Show a => DualQuaternion a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> DualQuaternion a -> ShowS
showsPrec :: Int -> DualQuaternion a -> ShowS
$cshow :: forall a. Show a => DualQuaternion a -> String
show :: DualQuaternion a -> String
$cshowList :: forall a. Show a => [DualQuaternion a] -> ShowS
showList :: [DualQuaternion a] -> ShowS
Show)

{- | Multiply real/dual quaternion pairs as (a,b)*(c,d)=(a*c,a*d+b*c). Return the product
using host quaternion arithmetic; no validation or normalization is added.
-}
dualQuaternionMultiply :: (RealFloat a) => DualQuaternion a -> DualQuaternion a -> DualQuaternion a
dualQuaternionMultiply :: forall a.
RealFloat a =>
DualQuaternion a -> DualQuaternion a -> DualQuaternion a
dualQuaternionMultiply (DualQuaternion Quaternion a
a Quaternion a
b) (DualQuaternion Quaternion a
c Quaternion a
d) = Quaternion a -> Quaternion a -> DualQuaternion a
forall a. Quaternion a -> Quaternion a -> DualQuaternion a
DualQuaternion (Quaternion a
a Quaternion a -> Quaternion a -> Quaternion a
forall a. Num a => a -> a -> a
* Quaternion a
c) (Quaternion a
a Quaternion a -> Quaternion a -> Quaternion a
forall a. Num a => a -> a -> a
* Quaternion a
d Quaternion a -> Quaternion a -> Quaternion a
forall a. Num a => a -> a -> a
+ Quaternion a
b Quaternion a -> Quaternion a -> Quaternion a
forall a. Num a => a -> a -> a
* Quaternion a
c)

{- | Invert the real quaternion a, returning (inverse(a), -inverse(a)*b*inverse(a)). Invalid
real part propagates quaternionInverse failure.
-}
dualQuaternionInverse :: (RealFloat a) => DualQuaternion a -> Either MathError (DualQuaternion a)
dualQuaternionInverse :: forall a.
RealFloat a =>
DualQuaternion a -> Either MathError (DualQuaternion a)
dualQuaternionInverse (DualQuaternion Quaternion a
a Quaternion a
b) = do ai <- Quaternion a -> Either MathError (Quaternion a)
forall a.
RealFloat a =>
Quaternion a -> Either MathError (Quaternion a)
quaternionInverse Quaternion a
a; Right (DualQuaternion ai (negate (ai * b * ai)))

-- | Construct a rigid transform from a nonzero rotation quaternion and translation.
rigidTransform :: (RealFloat a) => Quaternion a -> (a, a, a) -> Either MathError (DualQuaternion a)
rigidTransform :: forall a.
RealFloat a =>
Quaternion a -> (a, a, a) -> Either MathError (DualQuaternion a)
rigidTransform Quaternion a
q (a
x, a
y, a
z) = do
    _ <- Quaternion a -> Either MathError (Quaternion a)
forall a.
RealFloat a =>
Quaternion a -> Either MathError (Quaternion a)
quaternionInverse Quaternion a
q
    let Quaternion a b c d = q
        scale = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs [a
a, a
b, a
c, a
d])
        qs = Quaternion a -> a -> Quaternion a
forall a. Fractional a => Quaternion a -> a -> Quaternion a
qDivide Quaternion a
q a
scale
        r = Quaternion a -> a -> Quaternion a
forall a. Fractional a => Quaternion a -> a -> Quaternion a
qDivide Quaternion a
qs (Quaternion a -> a
forall a. RealFloat a => Quaternion a -> a
quaternionNorm Quaternion a
qs)
        dual = a -> a -> a -> a -> Quaternion a
forall a. a -> a -> a -> a -> Quaternion a
Quaternion a
0 (a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2) (a
y a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2) (a
z a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2) Quaternion a -> Quaternion a -> Quaternion a
forall a. Num a => a -> a -> a
* Quaternion a
r
    if qFinite r && qFinite dual then Right (DualQuaternion r dual) else Left NonFiniteResult

{- | Given dual-quaternion transform then point tuple, rotate by its real part and add
translation from twice the dual part times the real inverse. Invalid inverse or
non-finite result fails.
-}
transformPoint :: (RealFloat a) => DualQuaternion a -> (a, a, a) -> Either MathError (a, a, a)
transformPoint :: forall a.
RealFloat a =>
DualQuaternion a -> (a, a, a) -> Either MathError (a, a, a)
transformPoint (DualQuaternion Quaternion a
r Quaternion a
d) (a, a, a)
v = do
    (a, b, c) <- Quaternion a -> (a, a, a) -> Either MathError (a, a, a)
forall a.
RealFloat a =>
Quaternion a -> (a, a, a) -> Either MathError (a, a, a)
rotateVector Quaternion a
r (a, a, a)
v
    ri <- quaternionInverse r
    let Quaternion _ x y z = qScale 2 (d * ri)
    if all (\a
u -> Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
u Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
u)) [a + x, b + y, c + z]
        then Right (a + x, b + y, c + z)
        else Left NonFiniteResult

{- | Sum a finite list using Kahan compensation; empty input yields zero. Arithmetic follows
Num and compensation does not promise exact rounding or reject non-finite values.
-}
kahanSum :: (Num a) => [a] -> a
kahanSum :: forall a. Num a => [a] -> a
kahanSum [a]
xs = (a, a) -> a
forall a b. (a, b) -> a
fst (((a, a) -> a -> (a, a)) -> (a, a) -> [a] -> (a, a)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
L.foldl' (a, a) -> a -> (a, a)
forall {b}. Num b => (b, b) -> b -> (b, b)
step (a
0, a
0) [a]
xs)
  where
    step :: (b, b) -> b -> (b, b)
step (b
s, b
c) b
x = let y :: b
y = b
x b -> b -> b
forall a. Num a => a -> a -> a
- b
c; t :: b
t = b
s b -> b -> b
forall a. Num a => a -> a -> a
+ b
y in (b
t, (b
t b -> b -> b
forall a. Num a => a -> a -> a
- b
s) b -> b -> b
forall a. Num a => a -> a -> a
- b
y)

{- | Sum a finite list with Neumaier compensation, correcting for differing term magnitudes.
Empty input gives zero; arithmetic follows the chosen ordered Num type.
-}
neumaierSum :: (Num a, Ord a) => [a] -> a
neumaierSum :: forall a. (Num a, Ord a) => [a] -> a
neumaierSum [a]
xs = let (a
s, a
c) = ((a, a) -> a -> (a, a)) -> (a, a) -> [a] -> (a, a)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
L.foldl' (a, a) -> a -> (a, a)
forall {b}. (Ord b, Num b) => (b, b) -> b -> (b, b)
step (a
0, a
0) [a]
xs in a
s a -> a -> a
forall a. Num a => a -> a -> a
+ a
c
  where
    step :: (b, b) -> b -> (b, b)
step (b
s, b
c) b
x =
        let t :: b
t = b
s b -> b -> b
forall a. Num a => a -> a -> a
+ b
x; correction :: b
correction = if b -> b
forall a. Num a => a -> a
abs b
s b -> b -> Bool
forall a. Ord a => a -> a -> Bool
>= b -> b
forall a. Num a => a -> a
abs b
x then (b
s b -> b -> b
forall a. Num a => a -> a -> a
- b
t) b -> b -> b
forall a. Num a => a -> a -> a
+ b
x else (b
x b -> b -> b
forall a. Num a => a -> a -> a
- b
t) b -> b -> b
forall a. Num a => a -> a -> a
+ b
s
         in (b
t, b
c b -> b -> b
forall a. Num a => a -> a -> a
+ b
correction)

{- | Recursively split and sum a finite list to reduce sequential rounding accumulation.
Empty input gives zero. Arithmetic remains that of the chosen Num type.
-}
pairwiseSum :: (Num a) => [a] -> a
pairwiseSum :: forall a. Num a => [a] -> a
pairwiseSum [] = a
0
pairwiseSum [a
x] = a
x
pairwiseSum [a]
xs = let ([a]
a, [a]
b) = Int -> [a] -> ([a], [a])
forall a. Int -> [a] -> ([a], [a])
splitAt ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2) [a]
xs in [a] -> a
forall a. Num a => [a] -> a
pairwiseSum [a]
a a -> a -> a
forall a. Num a => a -> a -> a
+ [a] -> a
forall a. Num a => [a] -> a
pairwiseSum [a]
b

{- | Empty input is -infinity; +infinity dominates finite/-infinite inputs.
NaN is rejected, and an all -infinity input is -infinity.
-}
logsumexp :: (RealFloat a) => [a] -> Either MathError a
logsumexp :: forall a. RealFloat a => [a] -> Either MathError a
logsumexp [] = a -> Either MathError a
forall a b. b -> Either a b
Right (-a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0)
logsumexp [a]
xs
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any a -> Bool
forall a. RealFloat a => a -> Bool
isNaN [a]
xs = MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"logsumexp NaN input")
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0) [a]
xs = a -> Either MathError a
forall a b. b -> Either a b
Right (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0)
    | a
m a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== -a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0 = a -> Either MathError a
forall a b. b -> Either a b
Right a
m
    | Bool
otherwise = a -> Either MathError a
forall a b. b -> Either a b
Right (a
m a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
log ([a] -> a
forall a. (Num a, Ord a) => [a] -> a
neumaierSum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\a
x -> a -> a
forall a. Floating a => a -> a
exp (a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
m)) [a]
xs)))
  where
    m :: a
m = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [a]
xs

{- | Conservative sequential sum bound gamma_(n-1) * sum abs x. Returns infinity
when n*u >= 1 or finite inputs overflow the bound. Does not claim bounds in
the presence of underflow; an extra (n-1)*minimumSubnormal term covers it.
-}
sumErrorBound :: (RealFloat a) => [a] -> Either MathError a
sumErrorBound :: forall a. RealFloat a => [a] -> Either MathError a
sumErrorBound [a]
xs
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\a
x -> a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x) [a]
xs = MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"non-finite summation input")
    | [a] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
xs = a -> Either MathError a
forall a b. b -> Either a b
Right a
0
    | Bool
otherwise =
        let lim :: NumericLimits a
lim = a -> NumericLimits a
forall a. RealFloat a => a -> NumericLimits a
numericLimits ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
xs)
            n :: Integer
n = Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
            u :: Rational
u = a -> Rational
forall a. Real a => a -> Rational
toRational (NumericLimits a -> a
forall a. NumericLimits a -> a
machineEpsilon NumericLimits a
lim) Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
2
            nu :: Rational
nu = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
n Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
u
            absoluteSum :: Rational
absoluteSum = [Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> Rational) -> [a] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Rational
forall a. Num a => a -> a
abs (Rational -> Rational) -> (a -> Rational) -> a -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> Rational
forall a. Real a => a -> Rational
toRational) [a]
xs)
            exactBound :: Rational
exactBound = Rational
nu Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ (Rational
1 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- Rational
nu) Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
absoluteSum Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
n Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* a -> Rational
forall a. Real a => a -> Rational
toRational (NumericLimits a -> a
forall a. NumericLimits a -> a
minimumSubnormal NumericLimits a
lim)
            rounded :: a
rounded = Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
exactBound
         in a -> Either MathError a
forall a b. b -> Either a b
Right
                ( if Rational
nu Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
1 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
rounded
                    then a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0
                    else if a -> Rational
forall a. Real a => a -> Rational
toRational a
rounded Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
exactBound then a -> a -> a
forall a. RealFloat a => a -> a -> a
nextAfter a
rounded (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0) else a
rounded
                )