    {-# LANGUAGE BangPatterns #-}
    {-# LANGUAGE NoImplicitPrelude #-}
    module Numeric.Floating.IEEE.Internal.Augmented where
    import           Control.Exception (assert)
    import           MyPrelude
    import           Numeric.Floating.IEEE.Internal.FMA (isMantissaEven,
                                                         twoProduct_nonscaling,
                                                         twoSum)
    import           Numeric.Floating.IEEE.Internal.NextFloat (nextDown,
                                                               nextTowardZero,
                                                               nextUp)

    default ()

    -- $setup
    -- >>> :set -XDeriveFunctor
    -- >>> import Numeric.Floating.IEEE.Internal.Augmented
    -- >>> import Numeric.Floating.IEEE.Internal.Rounding.Common
    -- >>> import Numeric.Floating.IEEE.Internal.Rounding.Rational
    -- >>> :{
    -- newtype RoundTiesTowardZero a = RoundTiesTowardZero { roundTiesTowardZero :: a }
    --   deriving (Functor)
    -- instance RoundingStrategy RoundTiesTowardZero where
    --   exact = RoundTiesTowardZero
    --   inexact o _neg _parity zero away = RoundTiesTowardZero $ case o of
    --                                                              LT -> zero
    --                                                              EQ -> zero
    --                                                              GT -> away
    --   doRound _exact o _neg _parity zero away = RoundTiesTowardZero $ case o of
    --     LT -> zero
    --     EQ -> zero
    --     GT -> away
    -- :}

    -- |
    -- IEEE 754 @augmentedAddition@ operation.
    --
    -- The first return value is the approximation of the sum, and the second return value is the error.
    --
    -- prop> fst (augmentedAddition x y) == roundTiesTowardZero (fromRationalR (toRational x + toRational y)) `const` (x :: Double)
    -- prop> let (u, v) = augmentedAddition x y in toRational u + toRational v == toRational x + toRational y `const` (x :: Double)
    augmentedAddition :: RealFloat a => a -> a -> (a, a)
    augmentedAddition :: forall a. RealFloat a => a -> a -> (a, a)
    augmentedAddition !a
    x !a
    y
      | 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 Bool -> Bool -> Bool
    || a -> Bool
    forall a. RealFloat a => a -> Bool
    isNaN a
    y Bool -> Bool -> Bool
    || a -> Bool
    forall a. RealFloat a => a -> Bool
    isInfinite a
    y = let !result :: a
    result = a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    y in (a
    result, a
    result)
      | Bool
    otherwise = let (a
    u1, a
    u2) = a -> a -> (a, a)
    forall a. RealFloat a => a -> a -> (a, a)
    twoSum a
    x a
    y
                        ulpTowardZero :: a
    ulpTowardZero = a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a -> a
    forall a. RealFloat a => a -> a
    nextTowardZero a
    u1
                    in if a -> Bool
    forall a. RealFloat a => a -> Bool
    isNaN a
    u2 then
                         -- Handle undue overflow: e.g. 0x1.ffff_ffff_ffff_f8p1023
                         (a, a)
    handleUndueOverflow
                       else
                         if a
    u2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    0 then
                           (a
    u1, a
    0 a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    u1) -- signed zero
                         else
                           if (-a
    2) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    u2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    ulpTowardZero then
                             (a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    ulpTowardZero, a
    ulpTowardZero a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    u2)
                           else
                             (a
    u1, a
    u2)
      where
        handleUndueOverflow :: (a, a)
    handleUndueOverflow =
          -- The exponents of inputs should be close enough so that neither x' nor y' underflow.
          let e :: Int
    e = Int -> Int -> Int
    forall a. Ord a => a -> a -> a
    max (a -> Int
    forall a. RealFloat a => a -> Int
    exponent a
    x) (a -> Int
    forall a. RealFloat a => a -> Int
    exponent a
    y)
              x' :: a
    x' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat (- Int
    e) a
    x
              y' :: a
    y' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat (- Int
    e) a
    y
              (a
    u1, a
    u2) = a -> a -> (a, a)
    forall a. RealFloat a => a -> a -> (a, a)
    twoSum a
    x' a
    y'
              ulpTowardZero :: a
    ulpTowardZero = a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a -> a
    forall a. RealFloat a => a -> a
    nextTowardZero a
    u1
              (a
    v1, a
    v2) | (-a
    2) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    u2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    ulpTowardZero = (a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    ulpTowardZero, a
    ulpTowardZero a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    u2)
                       | Bool
    otherwise = (a
    u1, a
    u2)
              r1 :: a
    r1 = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    e a
    v1
              r2 :: a
    r2 = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    e a
    v2
          in if a -> Bool
    forall a. RealFloat a => a -> Bool
    isInfinite a
    r1 then
               (a
    r1, a
    r1) -- unavoidable overflow
             else
               Bool -> (a, a) -> (a, a)
    forall a. (?callStack::CallStack) => Bool -> a -> a
    assert (a
    r2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    /= a
    0) (a
    r1, a
    r2)
    {-# SPECIALIZE augmentedAddition :: Float -> Float -> (Float, Float) #-}
    {-# SPECIALIZE augmentedAddition :: Double -> Double -> (Double, Double) #-}

    -- |
    -- IEEE 754 @augmentedSubtraction@ operation.
    --
    -- The first return value is the approximation of the difference, and the second return value is the error.
    --
    -- prop> fst (augmentedSubtraction x y) == roundTiesTowardZero (fromRationalR (toRational x - toRational y)) `const` (x :: Double)
    -- prop> let (u, v) = augmentedSubtraction x y in toRational u + toRational v == toRational x - toRational y `const` (x :: Double)
    augmentedSubtraction :: RealFloat a => a -> a -> (a, a)
    augmentedSubtraction :: forall a. RealFloat a => a -> a -> (a, a)
    augmentedSubtraction a
    x a
    y = a -> a -> (a, a)
    forall a. RealFloat a => a -> a -> (a, a)
    augmentedAddition a
    x (a -> a
    forall a. Num a => a -> a
    negate a
    y)

    -- |
    -- IEEE 754 @augmentedMultiplication@ operation.
    --
    -- The first return value is the approximation of the product, and the second return value is the error.
    --
    -- prop> fst (augmentedMultiplication x y) == roundTiesTowardZero (fromRationalR (toRational x * toRational y)) `const` (x :: Double)
    -- prop> let (u, v) = augmentedMultiplication x y in toRational u + toRational v == toRational x * toRational y `const` (x :: Double)
    augmentedMultiplication :: RealFloat a => a -> a -> (a, a)
    augmentedMultiplication :: forall a. RealFloat a => a -> a -> (a, a)
    augmentedMultiplication !a
    x !a
    y
      | 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 Bool -> Bool -> Bool
    || a -> Bool
    forall a. RealFloat a => a -> Bool
    isNaN a
    y Bool -> Bool -> Bool
    || a -> Bool
    forall a. RealFloat a => a -> Bool
    isInfinite a
    y Bool -> Bool -> Bool
    || a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    y a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    0 = let !result :: a
    result = a
    x a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    y in (a
    result, a
    result)
      | Bool
    otherwise = let exy :: Int
    exy = a -> Int
    forall a. RealFloat a => a -> Int
    exponent a
    x Int -> Int -> Int
    forall a. Num a => a -> a -> a
    + a -> Int
    forall a. RealFloat a => a -> Int
    exponent a
    y
                        x' :: a
    x' = a -> a
    forall a. RealFloat a => a -> a
    significand a
    x
                        y' :: a
    y' = a -> a
    forall a. RealFloat a => a -> a
    significand a
    y
                        (a
    u1, a
    u2) = a -> a -> (a, a)
    forall a. RealFloat a => a -> a -> (a, a)
    twoProduct_nonscaling a
    x' a
    y'
                        !()
    _ = Bool -> () -> ()
    forall a. (?callStack::CallStack) => Bool -> a -> a
    assert (a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    x' Rational -> Rational -> Rational
    forall a. Num a => a -> a -> a
    * a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    y' Rational -> Rational -> Bool
    forall a. Eq a => a -> a -> Bool
    == a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    u1 Rational -> Rational -> Rational
    forall a. Num a => a -> a -> a
    + a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    u2) ()
                        -- The product is subnormal <=> exy + exponent u1 < expMin
                        -- The product is inexact => exy + exponent u1 < expMin + d
                    in if Int
    exy Int -> Int -> Int
    forall a. Num a => a -> a -> a
    + a -> Int
    forall a. RealFloat a => a -> Int
    exponent a
    u1 Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    >= Int
    expMin then
                         -- The result is exact
                         let ulpTowardZero :: a
    ulpTowardZero = a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a -> a
    forall a. RealFloat a => a -> a
    nextTowardZero a
    u1
                             !()
    _ = Bool -> () -> ()
    forall a. (?callStack::CallStack) => Bool -> a -> a
    assert (case a
    u1 of
                                            a
    0.5  -> - a
    ulpTowardZero a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <= a
    2 a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    u2 Bool -> Bool -> Bool
    && a
    u2 a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <= a
    ulpTowardZero
                                            -0.5 -> a
    ulpTowardZero a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <= a
    u2 Bool -> Bool -> Bool
    && a
    2 a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    u2 a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <= - a
    ulpTowardZero
                                            a
    _    -> a
    2 a -> a -> a
    forall a. Num a => a -> a -> a
    * a -> a
    forall a. Num a => a -> a
    abs a
    u2 a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    <= a -> a
    forall a. Num a => a -> a
    abs a
    ulpTowardZero) ()
                             (a
    v1, a
    v2) = if (-a
    2) a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    u2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    ulpTowardZero then
                                          (a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    ulpTowardZero, a
    ulpTowardZero a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    u2)
                                        else
                                          (a
    u1, a
    u2)
                             !()
    _ = Bool -> () -> ()
    forall a. (?callStack::CallStack) => Bool -> a -> a
    assert (a
    v1 a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    v2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    u2) ()
                             r1 :: a
    r1 = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    exy a
    v1
                             -- !_ = assert (r1 == roundTiesTowardZero (fromRationalR (toRational x * toRational y))) ()
                         in if a -> Bool
    forall a. RealFloat a => a -> Bool
    isInfinite a
    r1 then
                              (a
    r1, a
    r1)
                            else
                              if a
    v2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    0 then
                                (a
    r1, a
    0 a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    r1) -- signed zero
                              else
                                if Int
    exy Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    >= Int
    expMin Int -> Int -> Int
    forall a. Num a => a -> a -> a
    + Int
    d then
                                  -- The result is exact
                                  let r2 :: a
    r2 = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    exy a
    v2
                                  in (a
    r1, a
    r2)
                                else
                                  -- The upper part is normal, the lower is subnormal (and inexact)
                                  -- Compute 'scaleFloat exy v2' with roundTiesTowardZero
                                  let !r2 :: a
    r2 = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloatIntoSubnormalTiesTowardZero Int
    exy a
    v2
                                      -- !_ = assert (r2 == roundTiesTowardZero (fromRationalR (toRational x * toRational y - toRational r1))) ()
                                  in (a
    r1, a
    r2)
                       else
                         -- The upper part is subnormal (possibly inexact), and the lower is signed zero (possibly inexact)
                         if a
    u2 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    0 then
                           -- u1 is exact
                           let !()
    _ = Bool -> () -> ()
    forall a. (?callStack::CallStack) => Bool -> a -> a
    assert (a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    x' Rational -> Rational -> Rational
    forall a. Num a => a -> a -> a
    * a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    y' Rational -> Rational -> Bool
    forall a. Eq a => a -> a -> Bool
    == a -> Rational
    forall a. Real a => a -> Rational
    toRational a
    u1) ()
                               r1 :: a
    r1 = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloatIntoSubnormalTiesTowardZero Int
    exy a
    u1
                               r1' :: a
    r1' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat (-Int
    exy) a
    r1
                           in if a
    u1 a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    r1' then
                                (a
    r1, a
    0 a -> a -> a
    forall a. Num a => a -> a -> a
    * a
    r1)
                              else
                                (a
    r1, a
    0 a -> a -> a
    forall a. Num a => a -> a -> a
    * (a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    r1'))
                         else
                           let u1' :: a
    u1' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    exy a
    u1
                               v1' :: a
    v1' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    exy (if a
    u2 a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    > a
    0 then a -> a
    forall a. RealFloat a => a -> a
    nextUp a
    u1 else a -> a
    forall a. RealFloat a => a -> a
    nextDown a
    u1)
                               r1 :: a
    r1 = if a
    u1' a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    v1' Bool -> Bool -> Bool
    || Bool -> Bool
    not (a -> Bool
    forall a. RealFloat a => a -> Bool
    isMantissaEven a
    u1') then
                                      a
    u1'
                                    else
                                      a
    v1'
                               r1' :: a
    r1' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat (-Int
    exy) a
    r1
                           in (a
    r1, a
    0 a -> a -> a
    forall a. Num a => a -> a -> a
    * (a
    u1 a -> a -> a
    forall a. Num a => a -> a -> a
    - a
    r1' a -> a -> a
    forall a. Num a => a -> a -> a
    + a
    u2))
      where
        d :: Int
    d = a -> Int
    forall a. RealFloat a => a -> Int
    floatDigits a
    x
        (Int
    expMin,Int
    _expMax) = a -> (Int, Int)
    forall a. RealFloat a => a -> (Int, Int)
    floatRange a
    x

        -- Compute 'scaleFloat e z' with roundTiesTowardZero
        scaleFloatIntoSubnormalTiesTowardZero :: Int -> a -> a
    scaleFloatIntoSubnormalTiesTowardZero Int
    e a
    z =
          let z' :: a
    z' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    e a
    z
              w' :: a
    w' = Int -> a -> a
    forall a. RealFloat a => Int -> a -> a
    scaleFloat Int
    e (a -> a
    forall a. RealFloat a => a -> a
    nextTowardZero a
    z)
          in if a
    z' a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    w' Bool -> Bool -> Bool
    || Bool -> Bool
    not (a -> Bool
    forall a. RealFloat a => a -> Bool
    isMantissaEven a
    z') then
               a
    z'
             else
               a
    w'
    {-# SPECIALIZE augmentedMultiplication :: Float -> Float -> (Float, Float) #-}
    {-# SPECIALIZE augmentedMultiplication :: Double -> Double -> (Double, Double) #-}