{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE NoImplicitPrelude #-}
#if defined(HAS_FMA_PRIM)
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE UnboxedTuples #-}
#endif
module Numeric.Floating.IEEE.Internal.FMA
  ( isMantissaEven
  , twoSum
  , addToOdd
  , split
  , twoProductFloat_viaDouble
  , twoProduct
  , twoProduct_nonscaling
  , twoProductFloat
  , twoProductDouble
  , fusedMultiplyAddFloat_viaDouble
  , fusedMultiplyAdd
  , fusedMultiplyAddFloat
  , fusedMultiplyAddDouble
  ) where
import           Control.Exception (assert)
import           Data.Bits
import           GHC.Float.Compat (castDoubleToWord64, castFloatToWord32,
                                   double2Float, float2Double)
import           MyPrelude
import           Numeric.Floating.IEEE.Internal.Base (isDoubleBinary64,
                                                      isFloatBinary32, (^!))
import           Numeric.Floating.IEEE.Internal.Classify (isFinite)
import           Numeric.Floating.IEEE.Internal.NextFloat (nextDown, nextUp)
#if defined(HAS_FMA_PRIM)
import           GHC.Exts
#endif

default ()

-- $setup
-- >>> :set -XScopedTypeVariables
-- >>> import Numeric.Floating.IEEE.Internal.FMA

-- Assumption: input is finite
isMantissaEven :: RealFloat a => a -> Bool
isMantissaEven :: forall a. RealFloat a => a -> Bool
isMantissaEven a
0 = Bool
True
isMantissaEven a
x = let !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a -> Bool
forall a. RealFloat a => a -> Bool
isFinite a
x) ()
                       (Integer
m,Int
n) = a -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat a
x
                       d :: Int
d = a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
x
                       !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a -> Integer
forall a. RealFloat a => a -> Integer
floatRadix a
x Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Bool -> Bool -> Bool
&& Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< a -> Integer
forall a. RealFloat a => a -> Integer
floatRadix a
x Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
d) ()
                       (Int
expMin, Int
_expMax) = a -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange a
x
                       s :: Int
s = Int
expMin Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d)
                       !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a -> Bool
forall a. RealFloat a => a -> Bool
isDenormalized a
x Bool -> Bool -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0)) ()
                   in if Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 then
                        Integer -> Bool
forall a. Integral a => a -> Bool
even (Integer
m Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shiftR` Int
s)
                      else
                        Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
m
{-# NOINLINE [1] isMantissaEven #-}
{-# RULES
"isMantissaEven/Double"
  isMantissaEven = \x -> even (castDoubleToWord64 x)
"isMantissaEven/Float"
  isMantissaEven = \x -> even (castFloatToWord32 x)
  #-}

-- |
-- Returns @x := a + b@ and @x - \<the exact value of (a + b)\>@.
--
-- This function does not avoid undue overflow;
-- For example, the second component of
-- @twoSum (0x1.017bd555b0b1fp1022) (-0x1.fffffffffffffp1023)@
-- is a NaN.
--
-- prop> \(a :: Double) (b :: Double) -> let (_,expMax) = floatRange a in max (exponent a) (exponent b) < expMax ==> let (x, y) = twoSum a b in a + b == x && toRational a + toRational b == toRational x + toRational y
twoSum :: RealFloat a => a -> a -> (a, a)
twoSum :: forall a. RealFloat a => a -> a -> (a, a)
twoSum a
a a
b =
  let x :: a
x = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
b
      t :: a
t = a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
a
      y :: a
y = (a
a a -> a -> a
forall a. Num a => a -> a -> a
- (a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
t)) a -> a -> a
forall a. Num a => a -> a -> a
+ (a
b a -> a -> a
forall a. Num a => a -> a -> a
- a
t)
      {-
        Alternative:
         y = if abs b <= abs a then
               b - (x - a)
             else
               a - (x - b)
      -}
  in (a
x, a
y)
{-# SPECIALIZE twoSum :: Float -> Float -> (Float, Float), Double -> Double -> (Double, Double) #-}

-- |
-- Addition, with round to nearest odd floating-point number.
-- Like 'twoSum', this function does not handle undue overflow.
addToOdd :: RealFloat a => a -> a -> a
addToOdd :: forall a. RealFloat a => a -> a -> a
addToOdd a
x a
y = let (a
u, a
v) = a -> a -> (a, a)
forall a. RealFloat a => a -> a -> (a, a)
twoSum a
x a
y
                   result :: a
result | a -> Bool
forall a. RealFloat a => a -> Bool
isMantissaEven a
u Bool -> Bool -> Bool
&& a
v a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = a -> a
forall a. RealFloat a => a -> a
nextDown a
u
                          | a -> Bool
forall a. RealFloat a => a -> Bool
isMantissaEven a
u Bool -> Bool -> Bool
&& a
v a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = a -> a
forall a. RealFloat a => a -> a
nextUp a
u
                          | a -> Bool
forall a. RealFloat a => a -> Bool
isMantissaEven a
u Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
v Bool -> Bool -> Bool
&& Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
u) =
                              let v' :: a
v' = if a -> a
forall a. Num a => a -> a
abs a
y a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a -> a
forall a. Num a => a -> a
abs a
x then
                                         a
y a -> a -> a
forall a. Num a => a -> a -> a
- (a
u a -> a -> a
forall a. Num a => a -> a -> a
- a
x)
                                       else
                                         a
x a -> a -> a
forall a. Num a => a -> a -> a
- (a
u a -> a -> a
forall a. Num a => a -> a -> a
- a
y)
                              in if a
v' a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 then
                                   a -> a
forall a. RealFloat a => a -> a
nextDown a
u
                                 else if a
v' a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 then
                                        a -> a
forall a. RealFloat a => a -> a
nextUp a
u
                                      else
                                        a
u
                          | Bool
otherwise = a
u
                   !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
u Bool -> Bool -> Bool
|| a -> Rational
forall a. Real a => a -> Rational
toRational a
u Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== 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 Bool -> Bool -> Bool
|| Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isMantissaEven a
result)) ()
               in a
result
{-# SPECIALIZE addToOdd :: Float -> Float -> Float, Double -> Double -> Double #-}

-- This function doesn't handle overflow or underflow
split :: RealFloat a => a -> (a, a)
split :: forall a. RealFloat a => a -> (a, a)
split a
a =
  let c :: a
c = a
factor a -> a -> a
forall a. Num a => a -> a -> a
* a
a
      x :: a
x = a
c a -> a -> a
forall a. Num a => a -> a -> a
- (a
c a -> a -> a
forall a. Num a => a -> a -> a
- a
a)
      y :: a
y = a
a a -> a -> a
forall a. Num a => a -> a -> a
- a
x
  in (a
x, a
y)
  where factor :: a
factor = Integer -> a
forall a. Num a => Integer -> a
fromInteger (Integer -> a) -> Integer -> a
forall a b. (a -> b) -> a -> b
$ Integer
1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ a -> Integer
forall a. RealFloat a => a -> Integer
floatRadix a
a Integer -> Int -> Integer
^! ((a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
a 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)
  -- factor == 134217729 for Double, 4097 for Float
{-# SPECIALIZE split :: Float -> (Float, Float), Double -> (Double, Double) #-}

-- This function will be rewritten into fastTwoProduct{Float,Double} if fast FMA is available; the rewriting may change behavior regarding overflow.
-- TODO: subnormal behavior?
-- |
-- prop> \(a :: Double) (b :: Double) -> let (x, y) = twoProduct a b in a * b == x && fromRational (toRational a * toRational b - toRational x) == y
twoProduct :: RealFloat a => a -> a -> (a, a)
twoProduct :: forall a. RealFloat a => a -> a -> (a, a)
twoProduct a
a a
b =
  let eab :: Int
eab = a -> Int
forall a. RealFloat a => a -> Int
exponent a
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ a -> Int
forall a. RealFloat a => a -> Int
exponent a
b
      a' :: a
a' = a -> a
forall a. RealFloat a => a -> a
significand a
a
      b' :: a
b' = a -> a
forall a. RealFloat a => a -> a
significand a
b
      (a
ah, a
al) = a -> (a, a)
forall a. RealFloat a => a -> (a, a)
split a
a'
      (a
bh, a
bl) = a -> (a, a)
forall a. RealFloat a => a -> (a, a)
split a
b'
      x :: a
x = a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b -- Since 'significand' doesn't honor the sign of zero, we can't use @a' * b'@
      y' :: a
y' = a
al a -> a -> a
forall a. Num a => a -> a -> a
* a
bl a -> a -> a
forall a. Num a => a -> a -> a
- (Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (-Int
eab) a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
ah a -> a -> a
forall a. Num a => a -> a -> a
* a
bh a -> a -> a
forall a. Num a => a -> a -> a
- a
al a -> a -> a
forall a. Num a => a -> a -> a
* a
bh a -> a -> a
forall a. Num a => a -> a -> a
- a
ah a -> a -> a
forall a. Num a => a -> a -> a
* a
bl)
  in (a
x, Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
eab a
y')
{-# INLINABLE [1] twoProduct #-}

twoProductFloat_viaDouble :: Float -> Float -> (Float, Float)
twoProductFloat_viaDouble :: Float -> Float -> (Float, Float)
twoProductFloat_viaDouble Float
a Float
b =
  let x, y :: Float
      a', b', x' :: Double
      a' :: Double
a' = Float -> Double
float2Double Float
a
      b' :: Double
b' = Float -> Double
float2Double Float
b
      x' :: Double
x' = Double
a' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b'
      x :: Float
x = Double -> Float
double2Float Double
x'
      y :: Float
y = Double -> Float
double2Float (Double
x' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Float -> Double
float2Double Float
x)
  in (Float
x, Float
y)

-- This function will be rewritten into fastTwoProduct{Float,Double} if fast FMA is available; the rewriting may change behavior regarding overflow.
twoProduct_nonscaling :: RealFloat a => a -> a -> (a, a)
twoProduct_nonscaling :: forall a. RealFloat a => a -> a -> (a, a)
twoProduct_nonscaling a
a a
b =
  let (a
ah, a
al) = a -> (a, a)
forall a. RealFloat a => a -> (a, a)
split a
a
      (a
bh, a
bl) = a -> (a, a)
forall a. RealFloat a => a -> (a, a)
split a
b
      x :: a
x = a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b
      y :: a
y = a
al a -> a -> a
forall a. Num a => a -> a -> a
* a
bl a -> a -> a
forall a. Num a => a -> a -> a
- (a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
ah a -> a -> a
forall a. Num a => a -> a -> a
* a
bh a -> a -> a
forall a. Num a => a -> a -> a
- a
al a -> a -> a
forall a. Num a => a -> a -> a
* a
bh a -> a -> a
forall a. Num a => a -> a -> a
- a
ah a -> a -> a
forall a. Num a => a -> a -> a
* a
bl)
  in (a
x, a
y)
{-# NOINLINE [1] twoProduct_nonscaling #-}

twoProductFloat :: Float -> Float -> (Float, Float)
twoProductDouble :: Double -> Double -> (Double, Double)

#if defined(HAS_FMA_PRIM) && 0
-- Disabled for now: https://gitlab.haskell.org/ghc/ghc/-/issues/24160

twoProductFloat# :: Float# -> Float# -> (# Float#, Float# #)
twoProductFloat# x y = let !r = x `timesFloat#` y
                           !s = fmsubFloat# x y r
                       in (# r, s #)

twoProductDouble# :: Double# -> Double# -> (# Double#, Double# #)
twoProductDouble# x y = let !r = x *## y
                            !s = fmsubDouble# x y r
                        in (# r, s #)

#if defined(DONT_INLINE_FMA_PRIM)
{-# NOINLINE twoProductFloat# #-}
{-# NOINLINE twoProductDouble# #-}
#else
{-# INLINE twoProductFloat# #-}
{-# INLINE twoProductDouble# #-}
#endif

twoProductFloat (F# x) (F# y) = case twoProductFloat# x y of
                                  (# r, s #) -> (F# r, F# s)

twoProductDouble (D# x) (D# y) = case twoProductDouble# x y of
                                   (# r, s #) -> (D# r, D# s)

{-# INLINE twoProductFloat #-}
{-# INLINE twoProductDouble #-}

{-# RULES
"twoProduct/Float" twoProduct = twoProductFloat
"twoProduct/Double" twoProduct = twoProductDouble
"twoProduct_nonscaling/Float" twoProduct_nonscaling = twoProductFloat
"twoProduct_nonscaling/Double" twoProduct_nonscaling = twoProductDouble
  #-}

#elif defined(HAS_FAST_FMA) || defined(HAS_FMA_PRIM)

twoProductFloat x y = let !r = x * y
                          !s = fusedMultiplyAddFloat x y (-r)
                      in (r, s)

twoProductDouble x y = let !r = x * y
                           !s = fusedMultiplyAddDouble x y (-r)
                       in (r, s)

{-# RULES
"twoProduct/Float" twoProduct = twoProductFloat
"twoProduct/Double" twoProduct = twoProductDouble
"twoProduct_nonscaling/Float" twoProduct_nonscaling = twoProductFloat
"twoProduct_nonscaling/Double" twoProduct_nonscaling = twoProductDouble
  #-}

#else

twoProductFloat :: Float -> Float -> (Float, Float)
twoProductFloat = Float -> Float -> (Float, Float)
twoProductFloat_viaDouble
{-# INLINE twoProductFloat #-}

twoProductDouble :: Double -> Double -> (Double, Double)
twoProductDouble = Double -> Double -> (Double, Double)
forall a. RealFloat a => a -> a -> (a, a)
twoProduct
{-# INLINE twoProductDouble #-}

{-# RULES
"twoProduct/Float" twoProduct = twoProductFloat_viaDouble
"twoProduct_nonscaling/Float" twoProduct_nonscaling = twoProductFloat_viaDouble
  #-}
{-# SPECIALIZE twoProduct :: Double -> Double -> (Double, Double) #-}
{-# SPECIALIZE twoProduct_nonscaling :: Double -> Double -> (Double, Double) #-}

#endif

-- |
-- @'fusedMultiplyAdd' a b c@ computes @a * b + c@ as a single, ternary operation.
-- Rounding is done only once.
--
-- May make use of hardware FMA instructions if the target architecture has it; set @fma3@ package flag on x86 systems.
--
-- IEEE 754 @fusedMultiplyAdd@ operation.
--
-- prop> \(a :: Double) (b :: Double) (c :: Double) -> fusedMultiplyAdd a b c == fromRational (toRational a * toRational b + toRational c)
fusedMultiplyAdd :: RealFloat a => a -> a -> a -> a
fusedMultiplyAdd :: forall a. RealFloat a => a -> a -> a -> a
fusedMultiplyAdd a
a a
b a
c
  | a -> Bool
forall a. RealFloat a => a -> Bool
isFinite a
a Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
isFinite a
b Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
isFinite a
c =
    let eab :: Int
eab | a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
|| a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (a -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange a
a) Int -> Int -> Int
forall a. Num a => a -> a -> a
- a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
a -- reasonably small
            | Bool
otherwise = a -> Int
forall a. RealFloat a => a -> Int
exponent a
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ a -> Int
forall a. RealFloat a => a -> Int
exponent a
b
        ec :: Int
ec | a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (a -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange a
c) Int -> Int -> Int
forall a. Num a => a -> a -> a
- a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
c
           | Bool
otherwise = a -> Int
forall a. RealFloat a => a -> Int
exponent a
c

        -- Avoid overflow in twoProduct
        a' :: a
a' = a -> a
forall a. RealFloat a => a -> a
significand a
a
        b' :: a
b' = a -> a
forall a. RealFloat a => a -> a
significand a
b
        (a
x', a
y') = a -> a -> (a, a)
forall a. RealFloat a => a -> a -> (a, a)
twoProduct_nonscaling a
a' a
b'
        !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a -> Rational
forall a. Real a => a -> Rational
toRational a
a' Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* a -> Rational
forall a. Real a => a -> Rational
toRational a
b' Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== 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') ()

        -- Avoid overflow in twoSum
        e :: Int
e = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
eab Int
ec
        x :: a
x = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (Int
eab Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
e) a
x'
        y :: a
y = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (Int
eab Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
e) a
y'
        c'' :: a
c'' = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max ((Int, Int) -> Int
forall a b. (a, b) -> a
fst (a -> (Int, Int)
forall a. RealFloat a => a -> (Int, Int)
floatRange a
c) Int -> Int -> Int
forall a. Num a => a -> a -> a
- a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Int
ec Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
e) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
ec) a
c -- may be inexact

        (a
u1,a
u2) = a -> a -> (a, a)
forall a. RealFloat a => a -> a -> (a, a)
twoSum a
y a
c''
        (a
v1,a
v2) = a -> a -> (a, a)
forall a. RealFloat a => a -> a -> (a, a)
twoSum a
u1 a
x
        w :: a
w = a -> a -> a
forall a. RealFloat a => a -> a -> a
addToOdd a
u2 a
v2
        result0 :: a
result0 = a
v1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
w
        !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a
result0 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> a
forall a. Fractional a => Rational -> a
fromRational (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 -> Rational
forall a. Num a => a -> a -> a
+ a -> Rational
forall a. Real a => a -> Rational
toRational a
c'')) ()
        result :: a
result = Int -> a -> a
forall a. RealFloat a => Int -> a -> a
scaleFloat Int
e a
result0
        !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (a
result a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> a
forall a. Fractional a => Rational -> a
fromRational (a -> Rational
forall a. Real a => a -> Rational
toRational a
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* a -> Rational
forall a. Real a => a -> Rational
toRational a
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ a -> Rational
forall a. Real a => a -> Rational
toRational a
c) Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isDenormalized a
result) ()
    in if a
result0 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then
         -- We need to handle the sign of zero
         if a
c a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
&& a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0 Bool -> Bool -> Bool
&& a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0 then
           a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b -- let a * b underflow
         else
           a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
c -- -0 if both a * b and c are -0
       else
         if a -> Bool
forall a. RealFloat a => a -> Bool
isDenormalized a
result then
           -- The rounding in 'scaleFloat e result0' may yield an incorrect result.
           -- Take the slow path.
           case a -> Rational
forall a. Real a => a -> Rational
toRational a
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* a -> Rational
forall a. Real a => a -> Rational
toRational a
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ a -> Rational
forall a. Real a => a -> Rational
toRational a
c of
             Rational
0 -> a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
c -- This should be exact
             Rational
r -> Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
r
         else
           a
result
  | a -> Bool
forall a. RealFloat a => a -> Bool
isFinite a
a Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
isFinite a
b = a
c a -> a -> a
forall a. Num a => a -> a -> a
+ a
c -- c is +-Infinity or NaN
  | Bool
otherwise = a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
c -- Infinity or NaN
{-# INLINABLE [1] fusedMultiplyAdd #-} -- May be rewritten into a more efficient one

fusedMultiplyAddFloat_viaDouble :: Float -> Float -> Float -> Float
fusedMultiplyAddFloat_viaDouble :: Float -> Float -> Float -> Float
fusedMultiplyAddFloat_viaDouble Float
a Float
b Float
c
  | Float -> Bool
forall a. RealFloat a => a -> Bool
isFinite Float
a Bool -> Bool -> Bool
&& Float -> Bool
forall a. RealFloat a => a -> Bool
isFinite Float
b Bool -> Bool -> Bool
&& Float -> Bool
forall a. RealFloat a => a -> Bool
isFinite Float
c =
    let a', b', c' :: Double
        a' :: Double
a' = Float -> Double
float2Double Float
a
        b' :: Double
b' = Float -> Double
float2Double Float
b
        c' :: Double
c' = Float -> Double
float2Double Float
c
        ab :: Double
ab = Double
a' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b' -- exact
        !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Double -> Rational
forall a. Real a => a -> Rational
toRational Double
ab Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== 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
b') ()
        result :: Float
result = Double -> Float
double2Float (Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
addToOdd Double
ab Double
c')
        !()
_ = Bool -> () -> ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Float
result Float -> Float -> Bool
forall a. Eq a => a -> a -> Bool
== Rational -> Float
forall a. Fractional a => Rational -> a
fromRational (Float -> Rational
forall a. Real a => a -> Rational
toRational Float
a Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Float -> Rational
forall a. Real a => a -> Rational
toRational Float
b Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ Float -> Rational
forall a. Real a => a -> Rational
toRational Float
c)) ()
    in Float
result
  | Float -> Bool
forall a. RealFloat a => a -> Bool
isFinite Float
a Bool -> Bool -> Bool
&& Float -> Bool
forall a. RealFloat a => a -> Bool
isFinite Float
b = Float
c Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
c -- a * b is finite, but c is Infinity or NaN
  | Bool
otherwise = Float
a Float -> Float -> Float
forall a. Num a => a -> a -> a
* Float
b Float -> Float -> Float
forall a. Num a => a -> a -> a
+ Float
c
  where
    !() = if Bool
isFloatBinary32 then () else [Char] -> ()
forall a. (?callStack::CallStack) => [Char] -> a
error [Char]
"fusedMultiplyAdd/Float: Float must be IEEE binary32"
    !() = if Bool
isDoubleBinary64 then () else [Char] -> ()
forall a. (?callStack::CallStack) => [Char] -> a
error [Char]
"fusedMultiplyAdd/Float: Double must be IEEE binary64"

#if defined(HAS_FMA_PRIM)

#if defined(DONT_INLINE_FMA_PRIM)

fusedMultiplyAddFloat# :: Float# -> Float# -> Float# -> Float#
fusedMultiplyAddFloat# x y z = fmaddFloat# x y z
{-# NOINLINE fusedMultiplyAddFloat# #-}

fusedMultiplyAddDouble# :: Double# -> Double# -> Double# -> Double#
fusedMultiplyAddDouble# x y z = fmaddDouble# x y z
{-# NOINLINE fusedMultiplyAddDouble# #-}

fusedMultiplyAddFloat :: Float -> Float -> Float -> Float
fusedMultiplyAddFloat (F# x) (F# y) (F# z) = F# (fusedMultiplyAddFloat# x y z)

fusedMultiplyAddDouble :: Double -> Double -> Double -> Double
fusedMultiplyAddDouble (D# x) (D# y) (D# z) = D# (fusedMultiplyAddDouble# x y z)

#else

fusedMultiplyAddFloat :: Float -> Float -> Float -> Float
fusedMultiplyAddFloat (F# x) (F# y) (F# z) = F# (fmaddFloat# x y z)

fusedMultiplyAddDouble :: Double -> Double -> Double -> Double
fusedMultiplyAddDouble (D# x) (D# y) (D# z) = D# (fmaddDouble# x y z)

#endif

{-# INLINE fusedMultiplyAddFloat #-}
{-# INLINE fusedMultiplyAddDouble #-}

{-# RULES
"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat
"fusedMultiplyAdd/Double" fusedMultiplyAdd = fusedMultiplyAddDouble
  #-}

#elif defined(HAS_FAST_FMA)

foreign import ccall unsafe "hs_fusedMultiplyAddFloat"
  fusedMultiplyAddFloat :: Float -> Float -> Float -> Float
foreign import ccall unsafe "hs_fusedMultiplyAddDouble"
  fusedMultiplyAddDouble :: Double -> Double -> Double -> Double

{-# RULES
"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat
"fusedMultiplyAdd/Double" fusedMultiplyAdd = fusedMultiplyAddDouble
  #-}

#elif defined(USE_C99_FMA)

-- libm's fma might be implemented with hardware
foreign import ccall unsafe "fmaf"
  fusedMultiplyAddFloat :: Float -> Float -> Float -> Float
foreign import ccall unsafe "fma"
  fusedMultiplyAddDouble :: Double -> Double -> Double -> Double

{-# RULES
"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat
"fusedMultiplyAdd/Double" fusedMultiplyAdd = fusedMultiplyAddDouble
  #-}

#else

fusedMultiplyAddFloat :: Float -> Float -> Float -> Float
fusedMultiplyAddFloat = fusedMultiplyAddFloat_viaDouble
{-# INLINE fusedMultiplyAddFloat #-}

fusedMultiplyAddDouble :: Double -> Double -> Double -> Double
fusedMultiplyAddDouble = fusedMultiplyAdd -- generic implementation
{-# INLINE fusedMultiplyAddDouble #-}

{-# RULES
"fusedMultiplyAdd/Float" fusedMultiplyAdd = fusedMultiplyAddFloat_viaDouble
  #-}
{-# SPECIALIZE fusedMultiplyAdd :: Double -> Double -> Double -> Double #-}

#endif