1-----------------------------------------------------------------------------
2-- |
3-- Module    : Data.SBV.Core.Floating
4-- Copyright : (c) Levent Erkok
5-- License   : BSD3
6-- Maintainer: erkokl@gmail.com
7-- Stability : experimental
8--
9-- Implementation of floating-point operations mapping to SMT-Lib2 floats
10-----------------------------------------------------------------------------
11
12{-# LANGUAGE DefaultSignatures    #-}
13{-# LANGUAGE FlexibleContexts     #-}
14{-# LANGUAGE FlexibleInstances    #-}
15{-# LANGUAGE Rank2Types           #-}
16{-# LANGUAGE ScopedTypeVariables  #-}
17{-# LANGUAGE TypeApplications     #-}
18{-# LANGUAGE TypeFamilies         #-}
19{-# LANGUAGE TypeOperators        #-}
20{-# LANGUAGE UndecidableInstances #-}
21
22{-# OPTIONS_GHC -Wall -Werror -fno-warn-orphans #-}
23
24module Data.SBV.Core.Floating (
25         IEEEFloating(..), IEEEFloatConvertible(..)
26       , sFloatAsSWord32, sDoubleAsSWord64, sFloatingPointAsSWord
27       , sWord32AsSFloat, sWord64AsSDouble, sWordAsSFloatingPoint
28       , blastSFloat, blastSDouble,  blastSFloatingPoint
29       , sFloatAsComparableSWord32,  sDoubleAsComparableSWord64,  sFloatingPointAsComparableSWord
30       , sComparableSWord32AsSFloat, sComparableSWord64AsSDouble, sComparableSWordAsSFloatingPoint
31       ,
32       ) where
33
34import Data.Bits (testBit)
35import Data.Int  (Int8,  Int16,  Int32,  Int64)
36import Data.Word (Word8, Word16, Word32, Word64)
37
38import Data.Proxy
39
40import Data.SBV.Core.AlgReals (isExactRational)
41import Data.SBV.Core.Sized
42import Data.SBV.Core.SizedFloats
43
44import Data.SBV.Core.Data
45import Data.SBV.Core.Kind
46import Data.SBV.Core.Model
47import Data.SBV.Core.Symbolic (addSValOptGoal)
48
49import Data.SBV.Utils.Numeric
50
51import Data.Ratio
52
53import GHC.TypeLits
54
55import LibBF
56
57-- For doctest use only
58--
59-- $setup
60-- >>> :set -XTypeApplications
61-- >>> :set -XRankNTypes
62-- >>> :set -XScopedTypeVariables
63-- >>> import Data.SBV.Provers.Prover (prove)
64
65-- | A class of floating-point (IEEE754) operations, some of
66-- which behave differently based on rounding modes. Note that unless
67-- the rounding mode is concretely RoundNearestTiesToEven, we will
68-- not concretely evaluate these, but rather pass down to the SMT solver.
69class (SymVal a, RealFloat a) => IEEEFloating a where
70  -- | Compute the floating point absolute value.
71  fpAbs             ::                  SBV a -> SBV a
72
73  -- | Compute the unary negation. Note that @0 - x@ is not equivalent to @-x@ for floating-point, since @-0@ and @0@ are different.
74  fpNeg             ::                  SBV a -> SBV a
75
76  -- | Add two floating point values, using the given rounding mode
77  fpAdd             :: SRoundingMode -> SBV a -> SBV a -> SBV a
78
79  -- | Subtract two floating point values, using the given rounding mode
80  fpSub             :: SRoundingMode -> SBV a -> SBV a -> SBV a
81
82  -- | Multiply two floating point values, using the given rounding mode
83  fpMul             :: SRoundingMode -> SBV a -> SBV a -> SBV a
84
85  -- | Divide two floating point values, using the given rounding mode
86  fpDiv             :: SRoundingMode -> SBV a -> SBV a -> SBV a
87
88  -- | Fused-multiply-add three floating point values, using the given rounding mode. @fpFMA x y z = x*y+z@ but with only
89  -- one rounding done for the whole operation; not two. Note that we will never concretely evaluate this function since
90  -- Haskell lacks an FMA implementation.
91  fpFMA             :: SRoundingMode -> SBV a -> SBV a -> SBV a -> SBV a
92
93  -- | Compute the square-root of a float, using the given rounding mode
94  fpSqrt            :: SRoundingMode -> SBV a -> SBV a
95
96  -- | Compute the remainder: @x - y * n@, where @n@ is the truncated integer nearest to x/y. The rounding mode
97  -- is implicitly assumed to be @RoundNearestTiesToEven@.
98  fpRem             ::                  SBV a -> SBV a -> SBV a
99
100  -- | Round to the nearest integral value, using the given rounding mode.
101  fpRoundToIntegral :: SRoundingMode -> SBV a -> SBV a
102
103  -- | Compute the minimum of two floats, respects @infinity@ and @NaN@ values
104  fpMin             ::                  SBV a -> SBV a -> SBV a
105
106  -- | Compute the maximum of two floats, respects @infinity@ and @NaN@ values
107  fpMax             ::                  SBV a -> SBV a -> SBV a
108
109  -- | Are the two given floats exactly the same. That is, @NaN@ will compare equal to itself, @+0@ will /not/ compare
110  -- equal to @-0@ etc. This is the object level equality, as opposed to the semantic equality. (For the latter, just use '.=='.)
111  fpIsEqualObject   ::                  SBV a -> SBV a -> SBool
112
113  -- | Is the floating-point number a normal value. (i.e., not denormalized.)
114  fpIsNormal :: SBV a -> SBool
115
116  -- | Is the floating-point number a subnormal value. (Also known as denormal.)
117  fpIsSubnormal :: SBV a -> SBool
118
119  -- | Is the floating-point number 0? (Note that both +0 and -0 will satisfy this predicate.)
120  fpIsZero :: SBV a -> SBool
121
122  -- | Is the floating-point number infinity? (Note that both +oo and -oo will satisfy this predicate.)
123  fpIsInfinite :: SBV a -> SBool
124
125  -- | Is the floating-point number a NaN value?
126  fpIsNaN ::  SBV a -> SBool
127
128  -- | Is the floating-point number negative? Note that -0 satisfies this predicate but +0 does not.
129  fpIsNegative :: SBV a -> SBool
130
131  -- | Is the floating-point number positive? Note that +0 satisfies this predicate but -0 does not.
132  fpIsPositive :: SBV a -> SBool
133
134  -- | Is the floating point number -0?
135  fpIsNegativeZero :: SBV a -> SBool
136
137  -- | Is the floating point number +0?
138  fpIsPositiveZero :: SBV a -> SBool
139
140  -- | Is the floating-point number a regular floating point, i.e., not NaN, nor +oo, nor -oo. Normals or denormals are allowed.
141  fpIsPoint :: SBV a -> SBool
142
143  -- Default definitions. Minimal complete definition: None! All should be taken care by defaults
144  -- Note that we never evaluate FMA concretely, as there's no fma operator in Haskell
145  fpAbs              = lift1  FP_Abs             (Just abs)                Nothing
146  fpNeg              = lift1  FP_Neg             (Just negate)             Nothing
147  fpAdd              = lift2  FP_Add             (Just (+))                . Just
148  fpSub              = lift2  FP_Sub             (Just (-))                . Just
149  fpMul              = lift2  FP_Mul             (Just (*))                . Just
150  fpDiv              = lift2  FP_Div             (Just (/))                . Just
151  fpFMA              = lift3  FP_FMA             Nothing                   . Just
152  fpSqrt             = lift1  FP_Sqrt            (Just sqrt)               . Just
153  fpRem              = lift2  FP_Rem             (Just fpRemH)             Nothing
154  fpRoundToIntegral  = lift1  FP_RoundToIntegral (Just fpRoundToIntegralH) . Just
155  fpMin              = liftMM FP_Min             (Just fpMinH)             Nothing
156  fpMax              = liftMM FP_Max             (Just fpMaxH)             Nothing
157  fpIsEqualObject    = lift2B FP_ObjEqual        (Just fpIsEqualObjectH)   Nothing
158  fpIsNormal         = lift1B FP_IsNormal        fpIsNormalizedH
159  fpIsSubnormal      = lift1B FP_IsSubnormal     isDenormalized
160  fpIsZero           = lift1B FP_IsZero          (== 0)
161  fpIsInfinite       = lift1B FP_IsInfinite      isInfinite
162  fpIsNaN            = lift1B FP_IsNaN           isNaN
163  fpIsNegative       = lift1B FP_IsNegative      (\x -> x < 0 ||       isNegativeZero x)
164  fpIsPositive       = lift1B FP_IsPositive      (\x -> x >= 0 && not (isNegativeZero x))
165  fpIsNegativeZero x = fpIsZero x .&& fpIsNegative x
166  fpIsPositiveZero x = fpIsZero x .&& fpIsPositive x
167  fpIsPoint        x = sNot (fpIsNaN x .|| fpIsInfinite x)
168
169-- | SFloat instance
170instance IEEEFloating Float
171
172-- | SDouble instance
173instance IEEEFloating Double
174
175-- | Capture convertability from/to FloatingPoint representations.
176--
177-- Conversions to float: 'toSFloat' and 'toSDouble' simply return the
178-- nearest representable float from the given type based on the rounding
179-- mode provided. Similarly, 'toSFloatingPoint' converts to a generalized
180-- floating-point number with specified exponent and significand bith widths.
181--
182-- Conversions from float: 'fromSFloat', 'fromSDouble', 'fromSFloatingPoint' functions do
183-- the reverse conversion. However some care is needed when given values
184-- that are not representable in the integral target domain. For instance,
185-- converting an 'SFloat' to an 'SInt8' is problematic. The rules are as follows:
186--
187-- If the input value is a finite point and when rounded in the given rounding mode to an
188-- integral value lies within the target bounds, then that result is returned.
189-- (This is the regular interpretation of rounding in IEEE754.)
190--
191-- Otherwise (i.e., if the integral value in the float or double domain) doesn't
192-- fit into the target type, then the result is unspecified. Note that if the input
193-- is @+oo@, @-oo@, or @NaN@, then the result is unspecified.
194--
195-- Due to the unspecified nature of conversions, SBV will never constant fold
196-- conversions from floats to integral values. That is, you will always get a
197-- symbolic value as output. (Conversions from floats to other floats will be
198-- constant folded. Conversions from integral values to floats will also be
199-- constant folded.)
200--
201-- Note that unspecified really means unspecified: In particular, SBV makes
202-- no guarantees about matching the behavior between what you might get in
203-- Haskell, via SMT-Lib, or the C-translation. If the input value is out-of-bounds
204-- as defined above, or is @NaN@ or @oo@ or @-oo@, then all bets are off. In particular
205-- C and SMTLib are decidedly undefine this case, though that doesn't mean they do the
206-- same thing! Same goes for Haskell, which seems to convert via Int64, but we do
207-- not model that behavior in SBV as it doesn't seem to be intentional nor well documented.
208--
209-- You can check for @NaN@, @oo@ and @-oo@, using the predicates 'fpIsNaN', 'fpIsInfinite',
210-- and 'fpIsPositive', 'fpIsNegative' predicates, respectively; and do the proper conversion
211-- based on your needs. (0 is a good choice, as are min/max bounds of the target type.)
212--
213-- Currently, SBV provides no predicates to check if a value would lie within range for a
214-- particular conversion task, as this depends on the rounding mode and the types involved
215-- and can be rather tricky to determine. (See <http://github.com/LeventErkok/sbv/issues/456>
216-- for a discussion of the issues involved.) In a future release, we hope to be able to
217-- provide underflow and overflow predicates for these conversions as well.
218class SymVal a => IEEEFloatConvertible a where
219  -- | Convert from an IEEE74 single precision float.
220  fromSFloat :: SRoundingMode -> SFloat -> SBV a
221  fromSFloat = genericFromFloat
222
223  -- | Convert to an IEEE-754 Single-precision float.
224  --
225  -- >>> :{
226  -- roundTrip :: forall a. (Eq a, IEEEFloatConvertible a) => SRoundingMode -> SBV a -> SBool
227  -- roundTrip m x = fromSFloat m (toSFloat m x) .== x
228  -- :}
229  --
230  -- >>> prove $ roundTrip @Int8
231  -- Q.E.D.
232  -- >>> prove $ roundTrip @Word8
233  -- Q.E.D.
234  -- >>> prove $ roundTrip @Int16
235  -- Q.E.D.
236  -- >>> prove $ roundTrip @Word16
237  -- Q.E.D.
238  -- >>> prove $ roundTrip @Int32
239  -- Falsifiable. Counter-example:
240  --   s0 = RoundTowardNegative :: RoundingMode
241  --   s1 =           -44297675 :: Int32
242  --
243  -- Note how we get a failure on `Int32`. The counter-example value is not representable exactly as a single precision float:
244  --
245  -- >>> toRational (-44297675 :: Float)
246  -- (-44297676) % 1
247  --
248  -- Note how the numerator is different, it is off by 1. This is hardly surprising, since floats become sparser as
249  -- the magnitude increases to be able to cover all the integer values representable.
250  toSFloat :: SRoundingMode -> SBV a -> SFloat
251
252  -- default definition if we have an integral like
253  default toSFloat :: Integral a => SRoundingMode -> SBV a -> SFloat
254  toSFloat = genericToFloat (onlyWhenRNE (Just . fromRational . fromIntegral))
255
256  -- | Convert from an IEEE74 double precision float.
257  fromSDouble :: SRoundingMode -> SDouble -> SBV a
258  fromSDouble = genericFromFloat
259
260  -- | Convert to an IEEE-754 Double-precision float.
261  --
262  -- >>> :{
263  -- roundTrip :: forall a. (Eq a, IEEEFloatConvertible a) => SRoundingMode -> SBV a -> SBool
264  -- roundTrip m x = fromSDouble m (toSDouble m x) .== x
265  -- :}
266  --
267  -- >>> prove $ roundTrip @Int8
268  -- Q.E.D.
269  -- >>> prove $ roundTrip @Word8
270  -- Q.E.D.
271  -- >>> prove $ roundTrip @Int16
272  -- Q.E.D.
273  -- >>> prove $ roundTrip @Word16
274  -- Q.E.D.
275  -- >>> prove $ roundTrip @Int32
276  -- Q.E.D.
277  -- >>> prove $ roundTrip @Word32
278  -- Q.E.D.
279  -- >>> prove $ roundTrip @Int64
280  -- Falsifiable. Counter-example:
281  --   s0 = RoundNearestTiesToAway :: RoundingMode
282  --   s1 =   -3871901667041025751 :: Int64
283  --
284  -- Just like in the `SFloat` case, once we reach 64-bits, we no longer can exactly represent the
285  -- integer value for all possible values:
286  --
287  -- >>>  toRational (-3871901667041025751 :: Double)
288  -- (-3871901667041025536) % 1
289  --
290  -- In this case the numerator is off by 15.
291  toSDouble :: SRoundingMode -> SBV a -> SDouble
292
293  -- default definition if we have an integral like
294  default toSDouble :: Integral a => SRoundingMode -> SBV a -> SDouble
295  toSDouble = genericToFloat (onlyWhenRNE (Just . fromRational . fromIntegral))
296
297  -- | Convert from an arbitrary floating point.
298  fromSFloatingPoint :: ValidFloat eb sb => SRoundingMode -> SFloatingPoint eb sb -> SBV a
299  fromSFloatingPoint = genericFromFloat
300
301  -- | Convert to an arbitrary floating point.
302  toSFloatingPoint :: ValidFloat eb sb => SRoundingMode -> SBV a -> SFloatingPoint eb sb
303
304  -- -- default definition if we have an integral like
305  default toSFloatingPoint :: (Integral a, ValidFloat eb sb) => SRoundingMode -> SBV a -> SFloatingPoint eb sb
306  toSFloatingPoint = genericToFloat (const (Just . fromRational . fromIntegral))
307
308-- Run the function if the conversion is in RNE. Otherwise return Nothing.
309onlyWhenRNE :: (a -> Maybe b) -> RoundingMode -> a -> Maybe b
310onlyWhenRNE f RoundNearestTiesToEven v = f v
311onlyWhenRNE _ _                      _ = Nothing
312
313-- | A generic from-float converter. Note that this function does no constant folding since
314-- it's behavior is undefined when the input float is out-of-bounds or not a point.
315genericFromFloat :: forall a r. (IEEEFloating a, IEEEFloatConvertible r)
316                 => SRoundingMode            -- Rounding mode
317                 -> SBV a                    -- Input float/double
318                 -> SBV r
319genericFromFloat rm f = SBV (SVal kTo (Right (cache r)))
320  where kFrom = kindOf f
321        kTo   = kindOf (Proxy @r)
322        r st  = do msv <- sbvToSV st rm
323                   xsv <- sbvToSV st f
324                   newExpr st kTo (SBVApp (IEEEFP (FP_Cast kFrom kTo msv)) [xsv])
325
326-- | A generic to-float converter, which will constant-fold as necessary, but only in the sRNE mode for regular floats.
327genericToFloat :: forall a r. (IEEEFloatConvertible a, IEEEFloating r)
328               => (RoundingMode -> a -> Maybe r)     -- How to convert concretely, if possible
329               -> SRoundingMode                      -- Rounding mode
330               -> SBV a                              -- Input convertible
331               -> SBV r
332genericToFloat converter rm i
333  | Just w <- unliteral i, Just crm <- unliteral rm, Just result <- converter crm w
334  = literal result
335  | True
336  = SBV (SVal kTo (Right (cache r)))
337  where kFrom = kindOf i
338        kTo   = kindOf (Proxy @r)
339        r st  = do msv <- sbvToSV st rm
340                   xsv <- sbvToSV st i
341                   newExpr st kTo (SBVApp (IEEEFP (FP_Cast kFrom kTo msv)) [xsv])
342
343instance IEEEFloatConvertible Int8
344instance IEEEFloatConvertible Int16
345instance IEEEFloatConvertible Int32
346instance IEEEFloatConvertible Int64
347instance IEEEFloatConvertible Word8
348instance IEEEFloatConvertible Word16
349instance IEEEFloatConvertible Word32
350instance IEEEFloatConvertible Word64
351instance IEEEFloatConvertible Integer
352
353-- For float and double, skip the conversion if the same and do the constant folding, unlike all others.
354instance IEEEFloatConvertible Float where
355  toSFloat  _ f = f
356  toSDouble     = genericToFloat (onlyWhenRNE (Just . fp2fp))
357
358  toSFloatingPoint rm f = toSFloatingPoint rm $ toSDouble rm f
359
360  fromSFloat  _  f = f
361  fromSDouble rm f
362    | Just RoundNearestTiesToEven <- unliteral rm
363    , Just fv                     <- unliteral f
364    = literal (fp2fp fv)
365    | True
366    = genericFromFloat rm f
367
368instance IEEEFloatConvertible Double where
369  toSFloat      = genericToFloat (onlyWhenRNE (Just . fp2fp))
370  toSDouble _ d = d
371
372  toSFloatingPoint rm sd
373    | Just d <- unliteral sd, Just brm <- rmToRM rm
374    = literal $ FloatingPoint $ FP ei si $ fst (bfRoundFloat (mkBFOpts ei si brm) (bfFromDouble d))
375    | True
376    = res
377    where (k, ei, si) = case kindOf res of
378                         kr@(KFP eb sb) -> (kr, eb, sb)
379                         kr             -> error $ "Unexpected kind in toSFloatingPoint: " ++ show (kr, rm, sd)
380          res = SBV $ SVal k $ Right $ cache r
381          r st = do msv <- sbvToSV st rm
382                    xsv <- sbvToSV st sd
383                    newExpr st k (SBVApp (IEEEFP (FP_Cast KDouble k msv)) [xsv])
384
385  fromSDouble _  d = d
386  fromSFloat  rm d
387    | Just RoundNearestTiesToEven <- unliteral rm
388    , Just dv                     <- unliteral d
389    = literal (fp2fp dv)
390    | True
391    = genericFromFloat rm d
392
393convertWhenExactRational :: Fractional a => AlgReal -> Maybe a
394convertWhenExactRational r
395  | isExactRational r = Just (fromRational (toRational r))
396  | True              = Nothing
397
398-- For AlgReal; be careful to only process exact rationals concretely
399instance IEEEFloatConvertible AlgReal where
400  toSFloat         = genericToFloat (onlyWhenRNE convertWhenExactRational)
401  toSDouble        = genericToFloat (onlyWhenRNE convertWhenExactRational)
402  toSFloatingPoint = genericToFloat (const       convertWhenExactRational)
403
404-- Arbitrary floats can handle all rounding modes in concrete mode
405instance ValidFloat eb sb => IEEEFloatConvertible (FloatingPoint eb sb) where
406  toSFloat rm i
407    | Just (FloatingPoint (FP _ _ v)) <- unliteral i, Just brm <- rmToRM rm
408    = literal $ fp2fp $ fst (bfToDouble brm (fst (bfRoundFloat (mkBFOpts ei si brm) v)))
409    | True
410    = genericToFloat (\_ _ -> Nothing) rm i
411    where ei = intOfProxy (Proxy @eb)
412          si = intOfProxy (Proxy @sb)
413
414  fromSFloat rm i
415    | Just f <- unliteral i, Just brm <- rmToRM rm
416    = literal $ FloatingPoint $ FP ei si $ fst (bfRoundFloat (mkBFOpts ei si brm) (bfFromDouble (fp2fp f :: Double)))
417    | True
418    = genericFromFloat rm i
419    where ei = intOfProxy (Proxy @eb)
420          si = intOfProxy (Proxy @sb)
421
422  toSDouble rm i
423    | Just (FloatingPoint (FP _ _ v)) <- unliteral i, Just brm <- rmToRM rm
424    = literal $ fst (bfToDouble brm (fst (bfRoundFloat (mkBFOpts ei si brm) v)))
425    | True
426    = genericToFloat (\_ _ -> Nothing) rm i
427    where ei = intOfProxy (Proxy @eb)
428          si = intOfProxy (Proxy @sb)
429
430  fromSDouble rm i
431    | Just f <- unliteral i, Just brm <- rmToRM rm
432    = literal $ FloatingPoint $ FP ei si $ fst (bfRoundFloat (mkBFOpts ei si brm) (bfFromDouble f))
433    | True
434    = genericFromFloat rm i
435    where ei = intOfProxy (Proxy @eb)
436          si = intOfProxy (Proxy @sb)
437
438  toSFloatingPoint rm i
439    | Just (FloatingPoint (FP _ _ v)) <- unliteral i, Just brm <- rmToRM rm
440    = literal $ FloatingPoint $ FP ei si $ fst (bfRoundFloat (mkBFOpts ei si brm) v)
441    | True
442    = genericToFloat (\_ _ -> Nothing) rm i
443    where ei = intOfProxy (Proxy @eb)
444          si = intOfProxy (Proxy @sb)
445
446  -- From and To are the same when the source is an arbitrary float!
447  fromSFloatingPoint = toSFloatingPoint
448
449-- | Concretely evaluate one arg function, if rounding mode is RoundNearestTiesToEven and we have enough concrete data
450concEval1 :: SymVal a => Maybe (a -> a) -> Maybe SRoundingMode -> SBV a -> Maybe (SBV a)
451concEval1 mbOp mbRm a = do op <- mbOp
452                           v  <- unliteral a
453                           case unliteral =<< mbRm of
454                                   Nothing                     -> (Just . literal) (op v)
455                                   Just RoundNearestTiesToEven -> (Just . literal) (op v)
456                                   _                           -> Nothing
457
458-- | Concretely evaluate two arg function, if rounding mode is RoundNearestTiesToEven and we have enough concrete data
459concEval2 :: SymVal a => Maybe (a -> a -> a) -> Maybe SRoundingMode -> SBV a -> SBV a -> Maybe (SBV a)
460concEval2 mbOp mbRm a b = do op <- mbOp
461                             v1 <- unliteral a
462                             v2 <- unliteral b
463                             case unliteral =<< mbRm of
464                                     Nothing                     -> (Just . literal) (v1 `op` v2)
465                                     Just RoundNearestTiesToEven -> (Just . literal) (v1 `op` v2)
466                                     _                           -> Nothing
467
468-- | Concretely evaluate a bool producing two arg function, if rounding mode is RoundNearestTiesToEven and we have enough concrete data
469concEval2B :: SymVal a => Maybe (a -> a -> Bool) -> Maybe SRoundingMode -> SBV a -> SBV a -> Maybe SBool
470concEval2B mbOp mbRm a b = do op <- mbOp
471                              v1 <- unliteral a
472                              v2 <- unliteral b
473                              case unliteral =<< mbRm of
474                                      Nothing                     -> (Just . literal) (v1 `op` v2)
475                                      Just RoundNearestTiesToEven -> (Just . literal) (v1 `op` v2)
476                                      _                           -> Nothing
477
478-- | Concretely evaluate two arg function, if rounding mode is RoundNearestTiesToEven and we have enough concrete data
479concEval3 :: SymVal a => Maybe (a -> a -> a -> a) -> Maybe SRoundingMode -> SBV a -> SBV a -> SBV a -> Maybe (SBV a)
480concEval3 mbOp mbRm a b c = do op <- mbOp
481                               v1 <- unliteral a
482                               v2 <- unliteral b
483                               v3 <- unliteral c
484                               case unliteral =<< mbRm of
485                                       Nothing                     -> (Just . literal) (op v1 v2 v3)
486                                       Just RoundNearestTiesToEven -> (Just . literal) (op v1 v2 v3)
487                                       _                           -> Nothing
488
489-- | Add the converted rounding mode if given as an argument
490addRM :: State -> Maybe SRoundingMode -> [SV] -> IO [SV]
491addRM _  Nothing   as = return as
492addRM st (Just rm) as = do svm <- sbvToSV st rm
493                           return (svm : as)
494
495-- | Lift a 1 arg FP-op
496lift1 :: SymVal a => FPOp -> Maybe (a -> a) -> Maybe SRoundingMode -> SBV a -> SBV a
497lift1 w mbOp mbRm a
498  | Just cv <- concEval1 mbOp mbRm a
499  = cv
500  | True
501  = SBV $ SVal k $ Right $ cache r
502  where k    = kindOf a
503        r st = do sva  <- sbvToSV st a
504                  args <- addRM st mbRm [sva]
505                  newExpr st k (SBVApp (IEEEFP w) args)
506
507-- | Lift an FP predicate
508lift1B :: SymVal a => FPOp -> (a -> Bool) -> SBV a -> SBool
509lift1B w f a
510   | Just v <- unliteral a = literal $ f v
511   | True                  = SBV $ SVal KBool $ Right $ cache r
512   where r st = do sva <- sbvToSV st a
513                   newExpr st KBool (SBVApp (IEEEFP w) [sva])
514
515
516-- | Lift a 2 arg FP-op
517lift2 :: SymVal a => FPOp -> Maybe (a -> a -> a) -> Maybe SRoundingMode -> SBV a -> SBV a -> SBV a
518lift2 w mbOp mbRm a b
519  | Just cv <- concEval2 mbOp mbRm a b
520  = cv
521  | True
522  = SBV $ SVal k $ Right $ cache r
523  where k    = kindOf a
524        r st = do sva  <- sbvToSV st a
525                  svb  <- sbvToSV st b
526                  args <- addRM st mbRm [sva, svb]
527                  newExpr st k (SBVApp (IEEEFP w) args)
528
529-- | Lift min/max: Note that we protect against constant folding if args are alternating sign 0's, since
530-- SMTLib is deliberately nondeterministic in this case
531liftMM :: (SymVal a, RealFloat a) => FPOp -> Maybe (a -> a -> a) -> Maybe SRoundingMode -> SBV a -> SBV a -> SBV a
532liftMM w mbOp mbRm a b
533  | Just v1 <- unliteral a
534  , Just v2 <- unliteral b
535  , not ((isN0 v1 && isP0 v2) || (isP0 v1 && isN0 v2))          -- If not +0/-0 or -0/+0
536  , Just cv <- concEval2 mbOp mbRm a b
537  = cv
538  | True
539  = SBV $ SVal k $ Right $ cache r
540  where isN0   = isNegativeZero
541        isP0 x = x == 0 && not (isN0 x)
542        k    = kindOf a
543        r st = do sva  <- sbvToSV st a
544                  svb  <- sbvToSV st b
545                  args <- addRM st mbRm [sva, svb]
546                  newExpr st k (SBVApp (IEEEFP w) args)
547
548-- | Lift a 2 arg FP-op, producing bool
549lift2B :: SymVal a => FPOp -> Maybe (a -> a -> Bool) -> Maybe SRoundingMode -> SBV a -> SBV a -> SBool
550lift2B w mbOp mbRm a b
551  | Just cv <- concEval2B mbOp mbRm a b
552  = cv
553  | True
554  = SBV $ SVal KBool $ Right $ cache r
555  where r st = do sva  <- sbvToSV st a
556                  svb  <- sbvToSV st b
557                  args <- addRM st mbRm [sva, svb]
558                  newExpr st KBool (SBVApp (IEEEFP w) args)
559
560-- | Lift a 3 arg FP-op
561lift3 :: SymVal a => FPOp -> Maybe (a -> a -> a -> a) -> Maybe SRoundingMode -> SBV a -> SBV a -> SBV a -> SBV a
562lift3 w mbOp mbRm a b c
563  | Just cv <- concEval3 mbOp mbRm a b c
564  = cv
565  | True
566  = SBV $ SVal k $ Right $ cache r
567  where k    = kindOf a
568        r st = do sva  <- sbvToSV st a
569                  svb  <- sbvToSV st b
570                  svc  <- sbvToSV st c
571                  args <- addRM st mbRm [sva, svb, svc]
572                  newExpr st k (SBVApp (IEEEFP w) args)
573
574-- | Convert an 'SFloat' to an 'SWord32', preserving the bit-correspondence. Note that since the
575-- representation for @NaN@s are not unique, this function will return a symbolic value when given a
576-- concrete @NaN@.
577--
578-- Implementation note: Since there's no corresponding function in SMTLib for conversion to
579-- bit-representation due to partiality, we use a translation trick by allocating a new word variable,
580-- converting it to float, and requiring it to be equivalent to the input. In code-generation mode, we simply map
581-- it to a simple conversion.
582sFloatAsSWord32 :: SFloat -> SWord32
583sFloatAsSWord32 fVal
584  | Just f <- unliteral fVal, not (isNaN f)
585  = literal (floatToWord f)
586  | True
587  = SBV (SVal w32 (Right (cache y)))
588  where w32  = KBounded False 32
589        y st = do cg <- isCodeGenMode st
590                  if cg
591                     then do f <- sbvToSV st fVal
592                             newExpr st w32 (SBVApp (IEEEFP (FP_Reinterpret KFloat w32)) [f])
593                     else do n   <- internalVariable st w32
594                             ysw <- newExpr st KFloat (SBVApp (IEEEFP (FP_Reinterpret w32 KFloat)) [n])
595                             internalConstraint st False [] $ unSBV $ fVal `fpIsEqualObject` SBV (SVal KFloat (Right (cache (\_ -> return ysw))))
596                             return n
597
598-- | Convert an 'SDouble' to an 'SWord64', preserving the bit-correspondence. Note that since the
599-- representation for @NaN@s are not unique, this function will return a symbolic value when given a
600-- concrete @NaN@.
601--
602-- See the implementation note for 'sFloatAsSWord32', as it applies here as well.
603sDoubleAsSWord64 :: SDouble -> SWord64
604sDoubleAsSWord64 fVal
605  | Just f <- unliteral fVal, not (isNaN f)
606  = literal (doubleToWord f)
607  | True
608  = SBV (SVal w64 (Right (cache y)))
609  where w64  = KBounded False 64
610        y st = do cg <- isCodeGenMode st
611                  if cg
612                     then do f <- sbvToSV st fVal
613                             newExpr st w64 (SBVApp (IEEEFP (FP_Reinterpret KDouble w64)) [f])
614                     else do n   <- internalVariable st w64
615                             ysw <- newExpr st KDouble (SBVApp (IEEEFP (FP_Reinterpret w64 KDouble)) [n])
616                             internalConstraint st False [] $ unSBV $ fVal `fpIsEqualObject` SBV (SVal KDouble (Right (cache (\_ -> return ysw))))
617                             return n
618
619-- | Extract the sign\/exponent\/mantissa of a single-precision float. The output will have
620-- 8 bits in the second argument for exponent, and 23 in the third for the mantissa.
621blastSFloat :: SFloat -> (SBool, [SBool], [SBool])
622blastSFloat = extract . sFloatAsSWord32
623 where extract x = (sTestBit x 31, sExtractBits x [30, 29 .. 23], sExtractBits x [22, 21 .. 0])
624
625-- | Extract the sign\/exponent\/mantissa of a single-precision float. The output will have
626-- 11 bits in the second argument for exponent, and 52 in the third for the mantissa.
627blastSDouble :: SDouble -> (SBool, [SBool], [SBool])
628blastSDouble = extract . sDoubleAsSWord64
629 where extract x = (sTestBit x 63, sExtractBits x [62, 61 .. 52], sExtractBits x [51, 50 .. 0])
630
631-- | Extract the sign\/exponent\/mantissa of an arbitrary precision float. The output will have
632-- @eb@ bits in the second argument for exponent, and @sb-1@ bits in the third for mantissa.
633blastSFloatingPoint :: forall eb sb. (ValidFloat eb sb, KnownNat (eb + sb), BVIsNonZero (eb + sb))
634                    => SFloatingPoint eb sb -> (SBool, [SBool], [SBool])
635blastSFloatingPoint = extract . sFloatingPointAsSWord
636  where ei = intOfProxy (Proxy @eb)
637        si = intOfProxy (Proxy @sb)
638        extract x = (sTestBit x (ei + si - 1), sExtractBits x [ei + si - 2, ei + si - 3 .. si - 1], sExtractBits x [si - 2, si - 3 .. 0])
639
640-- | Reinterpret the bits in a 32-bit word as a single-precision floating point number
641sWord32AsSFloat :: SWord32 -> SFloat
642sWord32AsSFloat fVal
643  | Just f <- unliteral fVal = literal $ wordToFloat f
644  | True                     = SBV (SVal KFloat (Right (cache y)))
645  where y st = do xsv <- sbvToSV st fVal
646                  newExpr st KFloat (SBVApp (IEEEFP (FP_Reinterpret (kindOf fVal) KFloat)) [xsv])
647
648-- | Reinterpret the bits in a 32-bit word as a single-precision floating point number
649sWord64AsSDouble :: SWord64 -> SDouble
650sWord64AsSDouble dVal
651  | Just d <- unliteral dVal = literal $ wordToDouble d
652  | True                     = SBV (SVal KDouble (Right (cache y)))
653  where y st = do xsv <- sbvToSV st dVal
654                  newExpr st KDouble (SBVApp (IEEEFP (FP_Reinterpret (kindOf dVal) KDouble)) [xsv])
655
656-- | Convert a float to a comparable 'SWord32'. The trick is to ignore the
657-- sign of -0, and if it's a negative value flip all the bits, and otherwise
658-- only flip the sign bit. This is known as the lexicographic ordering on floats
659-- and it works as long as you do not have a @NaN@.
660sFloatAsComparableSWord32 :: SFloat -> SWord32
661sFloatAsComparableSWord32 f = ite (fpIsNegativeZero f) (sFloatAsComparableSWord32 0) (fromBitsBE $ sNot sb : ite sb (map sNot rest) rest)
662  where (sb : rest) = blastBE $ sFloatAsSWord32 f
663
664-- | Inverse transformation to 'sFloatAsComparableSWord32'. Note that this isn't a perfect inverse, since @-0@ maps to @0@ and back to @0@.
665-- Otherwise, it's faithful:
666--
667-- >>> prove  $ \x -> let f = sComparableSWord32AsSFloat x in fpIsNaN f .|| fpIsNegativeZero f .|| sFloatAsComparableSWord32 f .== x
668-- Q.E.D.
669-- >>> prove $ \x -> fpIsNegativeZero x .|| sComparableSWord32AsSFloat (sFloatAsComparableSWord32 x) `fpIsEqualObject` x
670-- Q.E.D.
671sComparableSWord32AsSFloat :: SWord32 -> SFloat
672sComparableSWord32AsSFloat w = sWord32AsSFloat $ ite sb (fromBitsBE $ sFalse : rest) (fromBitsBE $ map sNot allBits)
673  where allBits@(sb : rest) = blastBE w
674
675-- | Convert a double to a comparable 'SWord64'. The trick is to ignore the
676-- sign of -0, and if it's a negative value flip all the bits, and otherwise
677-- only flip the sign bit. This is known as the lexicographic ordering on doubles
678-- and it works as long as you do not have a @NaN@.
679sDoubleAsComparableSWord64 :: SDouble -> SWord64
680sDoubleAsComparableSWord64 d = ite (fpIsNegativeZero d) (sDoubleAsComparableSWord64 0) (fromBitsBE $ sNot sb : ite sb (map sNot rest) rest)
681  where (sb : rest) = blastBE $ sDoubleAsSWord64 d
682
683-- | Inverse transformation to 'sDoubleAsComparableSWord64'. Note that this isn't a perfect inverse, since @-0@ maps to @0@ and back to @0@.
684-- Otherwise, it's faithful:
685--
686-- >>> prove  $ \x -> let d = sComparableSWord64AsSDouble x in fpIsNaN d .|| fpIsNegativeZero d .|| sDoubleAsComparableSWord64 d .== x
687-- Q.E.D.
688-- >>> prove $ \x -> fpIsNegativeZero x .|| sComparableSWord64AsSDouble (sDoubleAsComparableSWord64 x) `fpIsEqualObject` x
689-- Q.E.D.
690sComparableSWord64AsSDouble :: SWord64 -> SDouble
691sComparableSWord64AsSDouble w = sWord64AsSDouble $ ite sb (fromBitsBE $ sFalse : rest) (fromBitsBE $ map sNot allBits)
692  where allBits@(sb : rest) = blastBE w
693
694-- | 'Float' instance for 'Metric' goes through the lexicographic ordering on 'Word32'.
695-- It implicitly makes sure that the value is not @NaN@.
696instance Metric Float where
697
698   type MetricSpace Float = Word32
699   toMetricSpace          = sFloatAsComparableSWord32
700   fromMetricSpace        = sComparableSWord32AsSFloat
701
702   msMinimize nm o = do constrain $ sNot $ fpIsNaN o
703                        addSValOptGoal $ unSBV `fmap` Minimize nm (toMetricSpace o)
704
705   msMaximize nm o = do constrain $ sNot $ fpIsNaN o
706                        addSValOptGoal $ unSBV `fmap` Maximize nm (toMetricSpace o)
707
708-- | 'Double' instance for 'Metric' goes through the lexicographic ordering on 'Word64'.
709-- It implicitly makes sure that the value is not @NaN@.
710instance Metric Double where
711
712   type MetricSpace Double = Word64
713   toMetricSpace           = sDoubleAsComparableSWord64
714   fromMetricSpace         = sComparableSWord64AsSDouble
715
716   msMinimize nm o = do constrain $ sNot $ fpIsNaN o
717                        addSValOptGoal $ unSBV `fmap` Minimize nm (toMetricSpace o)
718
719   msMaximize nm o = do constrain $ sNot $ fpIsNaN o
720                        addSValOptGoal $ unSBV `fmap` Maximize nm (toMetricSpace o)
721
722-- | Real instance for FloatingPoint. NB. The methods haven't been subjected to much testing, so beware of any floating-point snafus here.
723instance ValidFloat eb sb => Real (FloatingPoint eb sb) where
724  toRational (FloatingPoint (FP _ _ r)) = case bfToRep r of
725                                            BFNaN     -> toRational (0/0 :: Double)
726                                            BFRep s n -> case n of
727                                                           Zero    -> 0 % 1
728                                                           Inf     -> (if s == Neg then -1 else 1) % 0
729                                                           Num x y -> -- The value here is x * 2^y
730                                                                      let v :: Integer
731                                                                          v   = 2 ^ abs (fromIntegral y :: Integer)
732                                                                          sgn = if s == Neg then ((-1) *) else id
733                                                                      in if y > 0
734                                                                            then sgn $ x * v % 1
735                                                                            else sgn $ x % v
736
737-- | RealFrac instance for FloatingPoint. NB. The methods haven't been subjected to much testing, so beware of any floating-point snafus here.
738instance ValidFloat eb sb => RealFrac (FloatingPoint eb sb) where
739  properFraction (FloatingPoint f) = (a, FloatingPoint b)
740     where (a, b) = properFraction f
741
742-- | RealFloat instance for FloatingPoint. NB. The methods haven't been subjected to much testing, so beware of any floating-point snafus here.
743instance ValidFloat eb sb => RealFloat (FloatingPoint eb sb) where
744  floatRadix     (FloatingPoint f) = floatRadix     f
745  floatDigits    (FloatingPoint f) = floatDigits    f
746  floatRange     (FloatingPoint f) = floatRange     f
747  isNaN          (FloatingPoint f) = isNaN          f
748  isInfinite     (FloatingPoint f) = isInfinite     f
749  isDenormalized (FloatingPoint f) = isDenormalized f
750  isNegativeZero (FloatingPoint f) = isNegativeZero f
751  isIEEE         (FloatingPoint f) = isIEEE         f
752  decodeFloat    (FloatingPoint f) = decodeFloat    f
753
754  encodeFloat m n = res
755     where res = FloatingPoint $ fpEncodeFloat ei si m n
756           ei = intOfProxy (Proxy @eb)
757           si = intOfProxy (Proxy @sb)
758
759-- | Convert a float to the word containing the corresponding bit pattern
760sFloatingPointAsSWord :: forall eb sb. (ValidFloat eb sb, KnownNat (eb + sb), BVIsNonZero (eb + sb)) => SFloatingPoint eb sb -> SWord (eb + sb)
761sFloatingPointAsSWord fVal
762  | Just f@(FloatingPoint (FP eb sb v)) <- unliteral fVal, not (isNaN f)
763  = fromIntegral $ bfToBits (mkBFOpts eb sb NearEven) v
764  | True
765  = SBV (SVal kTo (Right (cache y)))
766  where ieb   = intOfProxy (Proxy @eb)
767        isb   = intOfProxy (Proxy @sb)
768        kFrom = KFP ieb isb
769        kTo   = KBounded False (ieb + isb)
770        y st = do cg <- isCodeGenMode st
771                  if cg
772                     then do f <- sbvToSV st fVal
773                             newExpr st kTo (SBVApp (IEEEFP (FP_Reinterpret kFrom kTo)) [f])
774                     else do n   <- internalVariable st kTo
775                             ysw <- newExpr st kFrom (SBVApp (IEEEFP (FP_Reinterpret kTo kFrom)) [n])
776                             internalConstraint st False [] $ unSBV $ fVal `fpIsEqualObject` SBV (SVal kFrom (Right (cache (\_ -> return ysw))))
777                             return n
778
779-- | Convert a float to the correct size word, that can be used in lexicographic ordering. Used in optimization.
780sFloatingPointAsComparableSWord :: forall eb sb. (ValidFloat eb sb, KnownNat (eb + sb), BVIsNonZero (eb + sb)) => SFloatingPoint eb sb -> SWord (eb + sb)
781sFloatingPointAsComparableSWord f = ite (fpIsNegativeZero f) posZero (fromBitsBE $ sNot sb : ite sb (map sNot rest) rest)
782  where posZero     = sFloatingPointAsComparableSWord (0 :: SFloatingPoint eb sb)
783        (sb : rest) = blastBE (sFloatingPointAsSWord f :: SWord (eb + sb))
784
785-- | Inverse transformation to 'sFloatingPointAsComparableSWord'. Note that this isn't a perfect inverse, since @-0@ maps to @0@ and back to @0@.
786-- Otherwise, it's faithful:
787--
788-- >>> prove  $ \x -> let d = sComparableSWordAsSFloatingPoint x in fpIsNaN d .|| fpIsNegativeZero d .|| sFloatingPointAsComparableSWord (d :: SFPHalf) .== x
789-- Q.E.D.
790-- >>> prove $ \x -> fpIsNegativeZero x .|| sComparableSWordAsSFloatingPoint (sFloatingPointAsComparableSWord x) `fpIsEqualObject` (x :: SFPHalf)
791-- Q.E.D.
792sComparableSWordAsSFloatingPoint :: forall eb sb. (KnownNat (eb + sb), BVIsNonZero (eb + sb), ValidFloat eb sb) => SWord (eb + sb) -> SFloatingPoint eb sb
793sComparableSWordAsSFloatingPoint w = sWordAsSFloatingPoint $ ite signBit (fromBitsBE $ sFalse : rest) (fromBitsBE $ map sNot allBits)
794  where allBits@(signBit : rest) = blastBE w
795
796-- | Convert a word to an arbitrary float, by reinterpreting the bits of the word as the corresponding bits of the float.
797sWordAsSFloatingPoint :: forall eb sb. (KnownNat (eb + sb), BVIsNonZero (eb + sb), ValidFloat eb sb) => SWord (eb + sb) -> SFloatingPoint eb sb
798sWordAsSFloatingPoint sw
799   | Just (f :: WordN (eb + sb)) <- unliteral sw
800   = let ext i = f `testBit` i
801         exts  = map ext
802         (s, ebits, sigbits) = (ext (ei + si - 1), exts [ei + si - 2, ei + si - 3 .. si - 1], exts [si - 2, si - 3 .. 0])
803
804         cvt :: [Bool] -> Integer
805         cvt = foldr (\b sofar -> 2 * sofar + if b then 1 else 0) 0 . reverse
806
807         eIntV = cvt ebits
808         sIntV = cvt sigbits
809         fp    = fpFromRawRep s (eIntV, ei) (sIntV, si)
810     in literal $ FloatingPoint fp
811   | True
812   = SBV (SVal kTo (Right (cache y)))
813   where ei   = intOfProxy (Proxy @eb)
814         si   = intOfProxy (Proxy @sb)
815         kTo  = KFP ei si
816         y st = do xsv <- sbvToSV st sw
817                   newExpr st kTo (SBVApp (IEEEFP (FP_Reinterpret (kindOf sw) kTo)) [xsv])
818
819instance (BVIsNonZero (eb + sb), KnownNat (eb + sb), ValidFloat eb sb) => Metric (FloatingPoint eb sb) where
820
821   type MetricSpace (FloatingPoint eb sb) = WordN (eb + sb)
822   toMetricSpace                          = sFloatingPointAsComparableSWord
823   fromMetricSpace                        = sComparableSWordAsSFloatingPoint
824
825   msMinimize nm o = do constrain $ sNot $ fpIsNaN o
826                        addSValOptGoal $ unSBV `fmap` Minimize nm (toMetricSpace o)
827
828   msMaximize nm o = do constrain $ sNot $ fpIsNaN o
829                        addSValOptGoal $ unSBV `fmap` Maximize nm (toMetricSpace o)
830
831-- Map SBV's rounding modes to LibBF's
832rmToRM :: SRoundingMode -> Maybe RoundMode
833rmToRM srm = cvt <$> unliteral srm
834  where cvt RoundNearestTiesToEven = NearEven
835        cvt RoundNearestTiesToAway = NearAway
836        cvt RoundTowardPositive    = ToPosInf
837        cvt RoundTowardNegative    = ToNegInf
838        cvt RoundTowardZero        = ToZero
839
840-- | Lift a 1 arg Big-float op
841lift1FP :: forall eb sb. ValidFloat eb sb =>
842           (BFOpts -> BigFloat -> (BigFloat, Status))
843        -> (Maybe SRoundingMode -> SFloatingPoint eb sb -> SFloatingPoint eb sb)
844        -> SRoundingMode
845        -> SFloatingPoint eb sb
846        -> SFloatingPoint eb sb
847lift1FP bfOp mkDef rm a
848  | Just (FloatingPoint (FP _ _ v)) <- unliteral a
849  , Just brm <- rmToRM rm
850  = literal $ FloatingPoint (FP ei si (fst (bfOp (mkBFOpts ei si brm) v)))
851  | True
852  = mkDef (Just rm) a
853  where ei = intOfProxy (Proxy @eb)
854        si = intOfProxy (Proxy @sb)
855
856-- | Lift a 2 arg Big-float op
857lift2FP :: forall eb sb. ValidFloat eb sb =>
858           (BFOpts -> BigFloat -> BigFloat -> (BigFloat, Status))
859        -> (Maybe SRoundingMode -> SFloatingPoint eb sb -> SFloatingPoint eb sb -> SFloatingPoint eb sb)
860        -> SRoundingMode
861        -> SFloatingPoint eb sb
862        -> SFloatingPoint eb sb
863        -> SFloatingPoint eb sb
864lift2FP bfOp mkDef rm a b
865  | Just (FloatingPoint (FP _ _ v1)) <- unliteral a
866  , Just (FloatingPoint (FP _ _ v2)) <- unliteral b
867  , Just brm <- rmToRM rm
868  = literal $ FloatingPoint (FP ei si (fst (bfOp (mkBFOpts ei si brm) v1 v2)))
869  | True
870  = mkDef (Just rm) a b
871  where ei = intOfProxy (Proxy @eb)
872        si = intOfProxy (Proxy @sb)
873
874-- | Lift a 3 arg Big-float op
875lift3FP :: forall eb sb. ValidFloat eb sb =>
876           (BFOpts -> BigFloat -> BigFloat -> BigFloat -> (BigFloat, Status))
877        -> (Maybe SRoundingMode -> SFloatingPoint eb sb -> SFloatingPoint eb sb -> SFloatingPoint eb sb -> SFloatingPoint eb sb)
878        -> SRoundingMode
879        -> SFloatingPoint eb sb
880        -> SFloatingPoint eb sb
881        -> SFloatingPoint eb sb
882        -> SFloatingPoint eb sb
883lift3FP bfOp mkDef rm a b c
884  | Just (FloatingPoint (FP _ _ v1)) <- unliteral a
885  , Just (FloatingPoint (FP _ _ v2)) <- unliteral b
886  , Just (FloatingPoint (FP _ _ v3)) <- unliteral c
887  , Just brm <- rmToRM rm
888  = literal $ FloatingPoint (FP ei si (fst (bfOp (mkBFOpts ei si brm) v1 v2 v3)))
889  | True
890  = mkDef (Just rm) a b c
891  where ei = intOfProxy (Proxy @eb)
892        si = intOfProxy (Proxy @sb)
893
894-- Sized-floats have a special instance, since it can handle arbitrary rounding modes when it matters.
895instance ValidFloat eb sb => IEEEFloating (FloatingPoint eb sb) where
896  fpAdd  = lift2FP bfAdd      (lift2 FP_Add  (Just (+)))
897  fpSub  = lift2FP bfSub      (lift2 FP_Sub  (Just (-)))
898  fpMul  = lift2FP bfMul      (lift2 FP_Mul  (Just (*)))
899  fpDiv  = lift2FP bfDiv      (lift2 FP_Div  (Just (/)))
900  fpFMA  = lift3FP bfFMA      (lift3 FP_FMA  Nothing)
901  fpSqrt = lift1FP bfSqrt     (lift1 FP_Sqrt (Just sqrt))
902
903  fpRoundToIntegral rm a
904    | Just (FloatingPoint (FP ei si v)) <- unliteral a
905    , Just brm <- rmToRM rm
906    = literal $ FloatingPoint (FP ei si (fst (bfRoundInt brm v)))
907    | True
908    = lift1 FP_RoundToIntegral (Just fpRoundToIntegralH) (Just rm) a
909
910  -- All other operations are agnostic to the rounding mode, hence the defaults are sufficient:
911  --
912  --       fpAbs            :: SBV a -> SBV a
913  --       fpNeg            :: SBV a -> SBV a
914  --       fpRem            :: SBV a -> SBV a -> SBV a
915  --       fpMin            :: SBV a -> SBV a -> SBV a
916  --       fpMax            :: SBV a -> SBV a -> SBV a
917  --       fpIsEqualObject  :: SBV a -> SBV a -> SBool
918  --       fpIsNormal       :: SBV a -> SBool
919  --       fpIsSubnormal    :: SBV a -> SBool
920  --       fpIsZero         :: SBV a -> SBool
921  --       fpIsInfinite     :: SBV a -> SBool
922  --       fpIsNaN          :: SBV a -> SBool
923  --       fpIsNegative     :: SBV a -> SBool
924  --       fpIsPositive     :: SBV a -> SBool
925  --       fpIsNegativeZero :: SBV a -> SBool
926  --       fpIsPositiveZero :: SBV a -> SBool
927  --       fpIsPoint        :: SBV a -> SBool
928
929{-# ANN module ("HLint: ignore Reduce duplication" :: String) #-}
930