{-# 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), 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), Double -> Double -> (Double, Double) #-}