what4-1.8: src/What4/SFloat.hs

{-# Language DataKinds #-}
{-# Language FlexibleContexts #-}
{-# Language GADTs #-}
{-# Language RankNTypes #-}
{-# Language TypeApplications #-}
{-# Language TypeOperators #-}

-- | Working with floats of dynamic sizes.
module What4.SFloat
  ( -- * Interface
    SFloat(..)
  , fpReprOf
  , fpSize
  , fpRepr
  , fpAsLit
  , fpIte

    -- * Constants
  , fpFresh
  , fpNaN
  , fpPosZero
  , fpNegZero
  , fpPosInf
  , fpNegInf
  , fpFromLit
  , fpFromRationalLit

    -- * Interchange formats
  , fpFromBinary
  , fpToBinary

    -- * Relations
  , SFloatRel
  , fpEq
  , fpNe
  , fpEqIEEE
  , fpNeIEEE
  , fpLeIEEE
  , fpLtIEEE
  , fpGeIEEE
  , fpGtIEEE
  , fpUnordered

    -- * Arithmetic
  , SFloatBinArith
  , fpNeg
  , fpAbs
  , fpSqrt
  , fpAdd
  , fpSub
  , fpMul
  , fpDiv
  , fpRem
  , fpMin
  , fpMax
  , fpFMA

    -- * Conversions
  , fpRound
  , fpToReal
  , fpFromReal
  , fpFromRational
  , fpToRational
  , fpFromInteger
  , fpCast
  , fpFromBV
  , fpFromSBV
  , fpToBV
  , fpToSBV

    -- * Queries
  , fpIsInf
  , fpIsNaN
  , fpIsZero
  , fpIsNeg
  , fpIsPos
  , fpIsSubnorm
  , fpIsNorm

  -- * Exceptions
  , UnsupportedFloat(..)
  , FPTypeError(..)
  ) where

import Control.Exception
import LibBF (BigFloat)
import Numeric.Natural

import Data.Parameterized.Some
import Data.Parameterized.NatRepr

import What4.BaseTypes
import What4.Panic(panic)
import What4.SWord
import What4.Interface

-- | Symbolic floating point numbers.
data SFloat sym where
  SFloat :: IsExpr (SymExpr sym) => SymFloat sym fpp -> SFloat sym



--------------------------------------------------------------------------------

-- | This exception is thrown if the operations try to create a
-- floating point value we do not support
data UnsupportedFloat =
  UnsupportedFloat { fpWho :: String, exponentBits, precisionBits :: Integer }
  deriving Show


-- | Throw 'UnsupportedFloat' exception
unsupported ::
  String  {- ^ Label -} ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  IO a
unsupported l e p =
  throwIO UnsupportedFloat { fpWho         = l
                           , exponentBits  = e
                           , precisionBits = p }

instance Exception UnsupportedFloat

-- | This exceptoin is throws if the types don't match.
data FPTypeError =
  FPTypeError { fpExpected :: Some BaseTypeRepr
              , fpActual   :: Some BaseTypeRepr
              }
    deriving Show

instance Exception FPTypeError

fpTypeMismatch :: BaseTypeRepr t1 -> BaseTypeRepr t2 -> IO a
fpTypeMismatch expect actual =
  throwIO FPTypeError { fpExpected = Some expect
                      , fpActual   = Some actual
                      }
fpTypeError :: FloatPrecisionRepr t1 -> FloatPrecisionRepr t2 -> IO a
fpTypeError t1 t2 =
  fpTypeMismatch (BaseFloatRepr t1) (BaseFloatRepr t2)


--------------------------------------------------------------------------------
-- | Construct the 'FloatPrecisionRepr' with the given parameters.
fpRepr ::
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  Maybe (Some FloatPrecisionRepr)
fpRepr iE iP =
  do Some e    <- someNat iE
     LeqProof  <- testLeq (knownNat @2) e
     Some p    <- someNat iP
     LeqProof  <- testLeq (knownNat @2) p
     pure (Some (FloatingPointPrecisionRepr e p))

fpReprOf ::
  IsExpr (SymExpr sym) => sym -> SymFloat sym fpp -> FloatPrecisionRepr fpp
fpReprOf _ e =
  case exprType e of
    BaseFloatRepr r -> r

fpSize :: SFloat sym -> (Integer,Integer)
fpSize (SFloat f) =
  case exprType f of
    BaseFloatRepr (FloatingPointPrecisionRepr e p) -> (intValue e, intValue p)

-- | See 'asFloat'.
fpAsLit :: SFloat sym -> Maybe BigFloat
fpAsLit (SFloat f) = asFloat f

--------------------------------------------------------------------------------
-- Constants

-- | A fresh variable of the given type (see 'freshConstant').
fpFresh ::
  IsSymExprBuilder sym =>
  sym ->
  Integer ->
  Integer ->
  IO (SFloat sym)
fpFresh sym e p
  | Just (Some fpp) <- fpRepr e p =
    SFloat <$> freshConstant sym emptySymbol (BaseFloatRepr fpp)
  | otherwise = unsupported "fpFresh" e p

-- | Not a number (see 'floatNaN').
fpNaN ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  IO (SFloat sym)
fpNaN sym e p
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatNaN sym fpp
  | otherwise = unsupported "fpNaN" e p

-- | Positive zero (see 'floatPZero').
fpPosZero ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  IO (SFloat sym)
fpPosZero sym e p
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatPZero sym fpp
  | otherwise = unsupported "fpPosZero" e p

-- | Negative zero (see 'floatNZero').
fpNegZero ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  IO (SFloat sym)
fpNegZero sym e p
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatNZero sym fpp
  | otherwise = unsupported "fpNegZero" e p

-- | Positive infinity (see 'floatPInf').
fpPosInf ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  IO (SFloat sym)
fpPosInf sym e p
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatPInf sym fpp
  | otherwise = unsupported "fpPosInf" e p

-- | Negative infinity (see 'floatNInf').
fpNegInf ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  IO (SFloat sym)
fpNegInf sym e p
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatNInf sym fpp
  | otherwise = unsupported "fpNegInf" e p


-- | A floating point number corresponding to the given BigFloat (see
-- 'floatLit').
fpFromLit ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  BigFloat ->
  IO (SFloat sym)
fpFromLit sym e p f
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatLit sym fpp f
  | otherwise = unsupported "fpFromLit" e p

-- | A floating point number corresponding to the given rational (see
-- 'floatLitRational').
fpFromRationalLit ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  Rational ->
  IO (SFloat sym)
fpFromRationalLit sym e p r
  | Just (Some fpp) <- fpRepr e p = SFloat <$> floatLitRational sym fpp r
  | otherwise = unsupported "fpFromRationalLit" e p


-- | Make a floating point number with the given bit representation (see
-- 'floatFromBinary').
fpFromBinary ::
  IsExprBuilder sym =>
  sym ->
  Integer {- ^ Exponent width -} ->
  Integer {- ^ Precision width -} ->
  SWord sym ->
  IO (SFloat sym)
fpFromBinary sym e p swe
  | DBV sw <- swe
  , Just (Some fpp) <- fpRepr e p
  , FloatingPointPrecisionRepr ew pw <- fpp
  , let expectW = addNat ew pw
  , actual@(BaseBVRepr actualW)  <- exprType sw =
    case testEquality expectW actualW of
      Just Refl -> SFloat <$> floatFromBinary sym fpp sw
      Nothing -- we want to report type correct type errors! :-)
        | Just LeqProof <- testLeq (knownNat @1) expectW ->
                fpTypeMismatch (BaseBVRepr expectW) actual
        | otherwise -> panic "fpFromBits" [ "1 >= 2" ]
  | otherwise = unsupported "fpFromBits" e p

-- | See 'floatToBinary'.
fpToBinary :: IsExprBuilder sym => sym -> SFloat sym -> IO (SWord sym)
fpToBinary sym (SFloat f)
  | FloatingPointPrecisionRepr e p <- fpReprOf sym f
  , Just LeqProof <- testLeq (knownNat @1) (addNat e p)
    = DBV <$> floatToBinary sym f
  | otherwise = panic "fpToBinary" [ "we messed up the types" ]


--------------------------------------------------------------------------------
-- Arithmetic

-- | See 'floatNeg'.
fpNeg :: IsExprBuilder sym => sym -> SFloat sym -> IO (SFloat sym)
fpNeg sym (SFloat fl) = SFloat <$> floatNeg sym fl

-- | See 'floatAbs'.
fpAbs :: IsExprBuilder sym => sym -> SFloat sym -> IO (SFloat sym)
fpAbs sym (SFloat fl) = SFloat <$> floatAbs sym fl

-- | See 'floatSqrt'.
fpSqrt :: IsExprBuilder sym => sym -> RoundingMode -> SFloat sym -> IO (SFloat sym)
fpSqrt sym r (SFloat fl) = SFloat <$> floatSqrt sym r fl

fpBinArith ::
  IsExprBuilder sym =>
  (forall t.
      sym ->
      RoundingMode ->
      SymFloat sym t ->
      SymFloat sym t ->
      IO (SymFloat sym t)
  ) ->
  sym -> RoundingMode -> SFloat sym -> SFloat sym -> IO (SFloat sym)
fpBinArith fun sym r (SFloat x) (SFloat y) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
  in
  case testEquality t1 t2 of
    Just Refl -> SFloat <$> fun sym r x y
    _         -> fpTypeError t1 t2

type SFloatBinArith sym =
  sym -> RoundingMode -> SFloat sym -> SFloat sym -> IO (SFloat sym)

-- | See 'floatAdd'.
fpAdd :: IsExprBuilder sym => SFloatBinArith sym
fpAdd = fpBinArith floatAdd

-- | See 'floatSub'.
fpSub :: IsExprBuilder sym => SFloatBinArith sym
fpSub = fpBinArith floatSub

-- | See 'floatMul'.
fpMul :: IsExprBuilder sym => SFloatBinArith sym
fpMul = fpBinArith floatMul

-- | See 'floatDiv'.
fpDiv :: IsExprBuilder sym => SFloatBinArith sym
fpDiv = fpBinArith floatDiv

-- | See 'floatRem'.
fpRem :: IsExprBuilder sym => sym -> SFloat sym -> SFloat sym -> IO (SFloat sym)
fpRem sym (SFloat x) (SFloat y) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
  in
  case testEquality t1 t2 of
    Just Refl -> SFloat <$> floatRem sym x y
    _         -> fpTypeError t1 t2

-- | See 'floatMin'.
fpMin :: IsExprBuilder sym => sym -> SFloat sym -> SFloat sym -> IO (SFloat sym)
fpMin sym (SFloat x) (SFloat y) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
  in
  case testEquality t1 t2 of
    Just Refl -> SFloat <$> floatMin sym x y
    _         -> fpTypeError t1 t2

-- | See 'floatMax'.
fpMax :: IsExprBuilder sym => sym -> SFloat sym -> SFloat sym -> IO (SFloat sym)
fpMax sym (SFloat x) (SFloat y) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
  in
  case testEquality t1 t2 of
    Just Refl -> SFloat <$> floatMax sym x y
    _         -> fpTypeError t1 t2

-- | See 'floatRMA'.
fpFMA :: IsExprBuilder sym =>
  sym -> RoundingMode -> SFloat sym -> SFloat sym -> SFloat sym -> IO (SFloat sym)
fpFMA sym r (SFloat x) (SFloat y) (SFloat z) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
      t3 = sym `fpReprOf` z
   in
   case (testEquality t1 t2, testEquality t2 t3) of
     (Just Refl, Just Refl) -> SFloat <$> floatFMA sym r x y z
     (Nothing, _) -> fpTypeError t1 t2
     (_, Nothing) -> fpTypeError t2 t3

-- | See 'floatIte'.
fpIte :: IsExprBuilder sym =>
  sym -> Pred sym -> SFloat sym -> SFloat sym -> IO (SFloat sym)
fpIte sym p (SFloat x) (SFloat y) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
  in
  case testEquality t1 t2 of
    Just Refl -> SFloat <$> floatIte sym p x y
    _         -> fpTypeError t1 t2

--------------------------------------------------------------------------------

fpRel ::
  IsExprBuilder sym =>
  (forall t.
    sym ->
    SymFloat sym t ->
    SymFloat sym t ->
    IO (Pred sym)
  ) ->
  sym -> SFloat sym -> SFloat sym -> IO (Pred sym)
fpRel fun sym (SFloat x) (SFloat y) =
  let t1 = sym `fpReprOf` x
      t2 = sym `fpReprOf` y
  in
  case testEquality t1 t2 of
    Just Refl -> fun sym x y
    _         -> fpTypeError t1 t2




type SFloatRel sym =
  sym -> SFloat sym -> SFloat sym -> IO (Pred sym)

-- | See 'floatEq'.
fpEq :: IsExprBuilder sym => SFloatRel sym
fpEq = fpRel floatEq

-- | See 'floatNe'.
fpNe :: IsExprBuilder sym => SFloatRel sym
fpNe = fpRel floatNe

-- | See 'floatFpEq'.
fpEqIEEE :: IsExprBuilder sym => SFloatRel sym
fpEqIEEE = fpRel floatFpEq

-- | See 'floatFpApart'.
fpNeIEEE :: IsExprBuilder sym => SFloatRel sym
fpNeIEEE = fpRel floatFpApart

-- | See 'floatLe'.
fpLeIEEE :: IsExprBuilder sym => SFloatRel sym
fpLeIEEE = fpRel floatLe

-- | See 'floatLt'.
fpLtIEEE :: IsExprBuilder sym => SFloatRel sym
fpLtIEEE = fpRel floatLt

-- | See 'floatGe'.
fpGeIEEE :: IsExprBuilder sym => SFloatRel sym
fpGeIEEE = fpRel floatGe

-- | See 'floatGt'.
fpGtIEEE :: IsExprBuilder sym => SFloatRel sym
fpGtIEEE = fpRel floatGt

-- | See 'floatFpUnordered'.
fpUnordered :: IsExprBuilder sym => SFloatRel sym
fpUnordered = fpRel floatFpUnordered

--------------------------------------------------------------------------------
-- | See 'floatRound'.
fpRound ::
  IsExprBuilder sym => sym -> RoundingMode -> SFloat sym -> IO (SFloat sym)
fpRound sym r (SFloat x) = SFloat <$> floatRound sym r x

-- | See 'floatToReal'. This is undefined on "special" values (NaN,infinity)
fpToReal :: IsExprBuilder sym => sym -> SFloat sym -> IO (SymReal sym)
fpToReal sym (SFloat x) = floatToReal sym x

-- | See 'realToFloat'.
fpFromReal ::
  IsExprBuilder sym =>
  sym -> Integer -> Integer -> RoundingMode -> SymReal sym -> IO (SFloat sym)
fpFromReal sym e p r x
  | Just (Some repr) <- fpRepr e p = SFloat <$> realToFloat sym repr r x
  | otherwise = unsupported "fpFromReal" e p


fpFromInteger ::
  IsExprBuilder sym =>
  sym -> Integer -> Integer -> RoundingMode -> SymInteger sym -> IO (SFloat sym)
fpFromInteger sym e p r x = fpFromReal sym e p r =<< integerToReal sym x


fpFromRational ::
  IsExprBuilder sym =>
  sym -> Integer -> Integer -> RoundingMode ->
  SymInteger sym -> SymInteger sym -> IO (SFloat sym)
fpFromRational sym e p r x y =
  do num <- integerToReal sym x
     den <- integerToReal sym y
     res <- realDiv sym num den
     fpFromReal sym e p r res

{- | Returns a predicate and two integers, @x@ and @y@.
If the the predicate holds, then @x / y@ is a rational representing
the floating point number. Assumes the FP number is not one of the
special ones that has no real representation. -}
fpToRational ::
  IsSymExprBuilder sym =>
  sym ->
  SFloat sym ->
  IO (Pred sym, SymInteger sym, SymInteger sym)
fpToRational sym fp =
  do r    <- fpToReal sym fp
     x    <- freshConstant sym emptySymbol BaseIntegerRepr
     y    <- freshConstant sym emptySymbol BaseIntegerRepr
     one  <- intLit sym 1
     -- Avoid rationals with zero denominators, which are invalid.
     yPos <- intLt sym one y
     num  <- integerToReal sym x
     den  <- integerToReal sym y
     res  <- realDiv sym num den
     same <- realEq sym r res
     rel  <- andPred sym yPos same
     pure (rel, x, y)

-- | Change the precision of a floating point number (see 'floatCast').
fpCast ::
  IsExprBuilder sym =>
  sym -> Integer -> Integer -> RoundingMode -> SFloat sym -> IO (SFloat sym)
fpCast sym e p r (SFloat x)
  | Just (Some repr) <- fpRepr e p = SFloat <$> floatCast sym repr r x
  | otherwise = unsupported "fpFromReal" e p

-- | Convert a unsigned bitvector to a floating point number (see 'bvToFloat').
fpFromBV ::
  IsExprBuilder sym =>
  sym -> Integer -> Integer -> RoundingMode -> SWord sym -> IO (SFloat sym)
fpFromBV sym e p r swe
  | DBV sw <- swe
  , Just (Some fpp) <- fpRepr e p =
    SFloat <$> bvToFloat sym fpp r sw
  | otherwise = unsupported "fpFromBV" e p

-- | Convert a signed bitvector to a floating point number (see 'sbvToFloat').
fpFromSBV ::
  IsExprBuilder sym =>
  sym -> Integer -> Integer -> RoundingMode -> SWord sym -> IO (SFloat sym)
fpFromSBV sym e p r swe
  | DBV sw <- swe
  , Just (Some fpp) <- fpRepr e p =
    SFloat <$> sbvToFloat sym fpp r sw
  | otherwise = unsupported "fpFromSBV" e p

-- | Convert a floating point number to a unsigned bitvector (see 'floatToBV').
-- Precondition: the supplied 'Natural' (the bit width of the returned
-- bitvector) is non-zero.
fpToBV ::
  IsExprBuilder sym =>
  sym -> Natural -> RoundingMode -> SFloat sym -> IO (SWord sym)
fpToBV sym n r (SFloat x) =
  case mkNatRepr n of
    Some nr ->
      case isPosNat nr of
        Nothing -> panic "fpToBV" ["bit width must be non-zero"]
        Just LeqProof -> DBV <$> floatToBV sym nr r x

-- | Convert a floating point number to a signed bitvector (see 'floatToSBV').
-- Precondition: the supplied 'Natural' (the bit width of the returned
-- bitvector) is non-zero.
fpToSBV ::
  IsExprBuilder sym =>
  sym -> Natural -> RoundingMode -> SFloat sym -> IO (SWord sym)
fpToSBV sym n r (SFloat x) =
  case mkNatRepr n of
    Some nr ->
      case isPosNat nr of
        Nothing -> panic "fpToSBV" ["bit width must be non-zero"]
        Just LeqProof -> DBV <$> floatToSBV sym nr r x

--------------------------------------------------------------------------------
-- | See 'floatIsInf'.
fpIsInf :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsInf sym (SFloat x) = floatIsInf sym x

-- | See 'floatIsNaN'.
fpIsNaN :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsNaN sym (SFloat x) = floatIsNaN sym x

-- | See 'floatIsZero'.
fpIsZero :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsZero sym (SFloat x) = floatIsZero sym x

-- | See 'floatIsNeg'.
fpIsNeg :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsNeg sym (SFloat x) = floatIsNeg sym x

-- | See 'floatIsPos'.
fpIsPos :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsPos sym (SFloat x) = floatIsPos sym x

-- | See 'floatIsSubnorm'.
fpIsSubnorm :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsSubnorm sym (SFloat x) = floatIsSubnorm sym x

-- | See 'floatIsNorm'.
fpIsNorm :: IsExprBuilder sym => sym -> SFloat sym -> IO (Pred sym)
fpIsNorm sym (SFloat x) = floatIsNorm sym x