FloatingPoint in Cryptol

FloatingPoint.cry

module FloatingPoint where

import Float

type FloatingPoint e m =
  { sign : Bit
  , exp : [e]
  , sfd : [m]
  }

type Single = FloatingPoint  8 23
type Double = FloatingPoint 11 52

// posZero, negZero

zeroFP : {e, m} Bit -> FloatingPoint e m
zeroFP sign =
  { sign = sign
  , exp = zero
  , sfd = zero
  }

posZero : {e, m} FloatingPoint e m
posZero = zeroFP False

posZeroFromFloat : Bit
property posZeroFromFloat = fromFloat fpPosZero == (posZero : Double)

posZeroToFloat : Bit
property posZeroToFloat = toFloat (posZero : Double) =.= fpPosZero

negZero : {e, m} FloatingPoint e m
negZero = zeroFP True

negZeroFromFloat : Bit
property negZeroFromFloat = fromFloat fpNegZero == (negZero : Double)

negZeroToFloat : Bit
property negZeroToFloat = toFloat (negZero : Double) =.= fpNegZero

// one

one : {e, m} (fin e) => Bit -> FloatingPoint e m
one sign =
  { sign = sign
  , exp = fromInteger (bias `e)
  , sfd = zero
  }

bias : Integer -> Integer
bias e = (2 ^^ (e - 1)) - 1

posOne : {e, m} (fin e) => FloatingPoint e m
posOne = one False

posOneFromFloat : Bit
property posOneFromFloat = fromFloat 1 == (posOne : Double)

posOneToFloat : Bit
property posOneToFloat = toFloat (posOne : Double) =.= 1

negOne : {e, m} (fin e) => FloatingPoint e m
negOne = one True

negOneFromFloat : Bit
property negOneFromFloat = fromFloat (-1) == (negOne : Double)

negOneToFloat : Bit
property negOneToFloat = toFloat (negOne : Double) =.= (-1)

// infinity

infinity : {e, m} Bit -> FloatingPoint e m
infinity sign =
  { sign = sign
  , exp = complement zero
  , sfd = zero
  }

posInf : {e, m} FloatingPoint e m
posInf = infinity False

posInfFromFloat : Bit
property posInfFromFloat = fromFloat fpPosInf == (posInf : Double)

posInfToFloat : Bit
property posInfToFloat = toFloat (posInf : Double) =.= fpPosInf

negInf : {e, m} FloatingPoint e m
negInf = infinity True

negInfFromFloat : Bit
property negInfFromFloat = fromFloat fpNegInf == (negInf : Double)

negInfToFloat : Bit
property negInfToFloat = toFloat (negInf : Double) =.= fpNegInf

// qnan

// LibBF encodes all NaN values as quiet NaNs.
qnan : {e, m} (m >= 1) => FloatingPoint e m
qnan =
  { sign = False
  , exp = complement zero
  , sfd = 0b1 # zero
  }

qnanFromFloat : Bit
qnanFromFloat = fromFloat fpNaN == (qnan : Double)

qnanToFloat : Bit
qnanToFloat = toFloat (qnan : Double) =.= fpNaN

// isInf

isNaNOrInf : {e, m} (fin e) => FloatingPoint e m -> Bit
isNaNOrInf fp = and fp.exp

isInf : {e, m} (fin e, fin m) => FloatingPoint e m -> Bit
isInf fp = and fp.exp /\ fp.sfd == 0

isInfFromFloat : Float64 -> Bit
isInfFromFloat f = isInf (fromFloat f) == fpIsInf f

isInfToFloat : Double -> Bit
isInfToFloat fp = fpIsInf (toFloat fp) == isInf fp

// isNaN

isNaN : {e, m} (fin e, fin m) => FloatingPoint e m -> Bit
isNaN fp = isNaNOrInf fp /\ ~(isInf fp)

isNaNFromFloat : Float64 -> Bit
isNaNFromFloat f = isNaN (fromFloat f) == fpIsNaN f

isNaNToFloat : Double -> Bit
isNaNToFloat fp = fpIsNaN (toFloat fp) == isNaN fp

// isSubNormal

// NB: This definition treats zero as a subnormal value, but many other FP
// implementations (e.g., LibBF) do not. For this reason, we also offer
// `isSubNormalStrict`, which excludes zero.
isSubNormal : {e, m} (fin e) => FloatingPoint e m -> Bit
isSubNormal fp = fp.exp == 0

// Like `isSubNormal`, but return `False` for zero values.
isSubNormalStrict : {e, m} (fin e, fin m) => FloatingPoint e m -> Bit
isSubNormalStrict fp = fp.exp == 0 /\ fp.sfd != 0

isSubNormalFromFloat : Float64 -> Bit
isSubNormalFromFloat f = isSubNormalStrict (fromFloat f) == fpIsSubnormal f

isSubNormalToFloat : Double -> Bit
isSubNormalToFloat fp = fpIsSubnormal (toFloat fp) == isSubNormalStrict fp

// isZero

isZero : {e, m} (fin e, fin m) => FloatingPoint e m -> Bit
isZero fp = isSubNormal fp /\ fp.sfd == 0

isZeroFromFloat : Float64 -> Bit
isZeroFromFloat f = isZero (fromFloat f) == fpIsZero f

isZeroToFloat : Double -> Bit
isZeroToFloat fp = fpIsZero (toFloat fp) == isZero fp

// compare

enum Disorder
  = LT
  | EQ
  | GT
  | UO

compare : {e, m} (fin e, fin m) => FloatingPoint e m -> FloatingPoint e m -> Disorder
compare x y =
  if isNaN x \/ isNaN y            then UO
   | isZero x /\ isZero y          then EQ
   | x.sign /\ ~y.sign             then LT
   | ~x.sign /\ y.sign             then GT
   | x.sign then
     (if expLT \/ (expEQ /\ sfdLT) then LT
       | expGT \/ (expEQ /\ sfdGT) then GT
      else                              EQ)
  else
     (if expGT \/ (expEQ /\ sfdGT) then LT
       | expLT \/ (expEQ /\ sfdLT) then GT
      else                              EQ)
  where
    expLT  = x.exp < y.exp
    expEQ  = x.exp == y.exp
    expGT  = x.exp > y.exp
    sfdLT  = x.sfd < y.sfd
    sfdGT  = x.sfd > y.sfd
    sfdEQ  = x.sfd == y.sfd

// eq

eq : {e, m} (fin e, fin m) => FloatingPoint e m -> FloatingPoint e m -> Bit
eq x y =
  case compare x y of
    EQ -> True
    _  -> False

eqFromFloat : Float64 -> Float64 -> Bit
property eqFromFloat x y = eq (fromFloat x) (fromFloat y) == (x == y)

eqToFloat : Double -> Double -> Bit
property eqToFloat x y = (toFloat x == toFloat y) == eq x y

// eqBits

// A more descriptive name for how (==) works on FloatingPoint values.
eqBits : {e, m} (fin e, fin m) => FloatingPoint e m -> FloatingPoint e m -> Bit
eqBits x y = x == y

eqBitsFromFloat : Float64 -> Float64 -> Bit
property eqBitsFromFloat x y = eqBits (fromFloat x) (fromFloat y) == (x =.= y)

// eqBitsModuloNaNs

// LibBF encodes all NaN values as quiet NaNs, but this Cryptol encoding can
// express many more NaN values. This notion of equality is like `eq`, but where
// all NaN values are considered equal.
eqBitsModuloNaNs : {e, m} (fin e, fin m) => FloatingPoint e m -> FloatingPoint e m -> Bit
eqBitsModuloNaNs x y =
  if isNaN x
  then isNaN y
  else eqBits x y

eqBitsModuloNaNsFromFloat : Float64 -> Float64 -> Bit
property eqBitsModuloNaNsFromFloat x y = eqBitsModuloNaNs (fromFloat x) (fromFloat y) == (x =.= y)

eqBitsModuloNaNsToFloat : Double -> Double -> Bit
property eqBitsModuloNaNsToFloat x y = (toFloat x =.= toFloat y) == eqBitsModuloNaNs x y

// fromBits, toBits

fromBits : {e, m} (fin e) => [1 + e + m] -> FloatingPoint e m
fromBits bits =
  { sign = sign
  , exp = exp
  , sfd = sfd
  }
  where
    ([sign] # exp # sfd) = bits

toBits : {e, m} (fin e) => FloatingPoint e m -> [1 + e + m]
toBits fp = [fp.sign] # fp.exp # fp.sfd

fromToBits : Double -> Bit
fromToBits fp = fromBits (toBits fp) == fp

toFromBits : [64] -> Bit
toFromBits bits = toBits (fromBits bits : Double) == bits

// fromFloat, toFloat

fromFloat : {e, m} (ValidFloat e (1 + m)) => Float e (1 + m) -> FloatingPoint e m
fromFloat f = fromBits (fpToBits f)

toFloat : {e, m} (ValidFloat e (1 + m)) => FloatingPoint e m -> Float e (1 + m)
toFloat fp = fpFromBits (toBits fp)

fromToFloat : Double -> Bit
property fromToFloat fp = eqBitsModuloNaNs fp (fromFloat (toFloat fp))

toFromFloat : Float64 -> Bit
property toFromFloat f = toFloat (fromFloat f) =.= f

// absFP

absFP : {e, m} FloatingPoint e m -> FloatingPoint e m
absFP fp = { fp | sign = False }

absFromFloat : Float64 -> Bit
property absFromFloat f = absFP (fromFloat f) == fromFloat (fpAbs f)

absToFloat : Double -> Bit
property absToFloat fp = fpAbs (toFloat fp) =.= toFloat (absFP fp)

// RoundMode

enum RoundMode
  = Rnd_Nearest_Even
  | Rnd_Nearest_Away_Zero
  | Rnd_Plus_Inf
  | Rnd_Minus_Inf
  | Rnd_Zero

eqRoundMode : RoundMode -> RoundMode -> Bit
eqRoundMode rm1 rm2 =
  case rm1 of
    Rnd_Nearest_Even ->
      case rm2 of
        Rnd_Nearest_Even -> True
        _ -> False
    Rnd_Nearest_Away_Zero ->
      case rm2 of
        Rnd_Nearest_Away_Zero -> True
        _ -> False
    Rnd_Plus_Inf ->
      case rm2 of
        Rnd_Plus_Inf -> True
        _ -> False
    Rnd_Minus_Inf ->
      case rm2 of
        Rnd_Minus_Inf -> True
        _ -> False
    Rnd_Zero ->
      case rm2 of
        Rnd_Zero -> True
        _ -> False

eqRoundModeReflexive : RoundMode -> Bit
property eqRoundModeReflexive x = eqRoundMode x x

eqRoundModeSymmetric : RoundMode -> RoundMode -> Bit
property eqRoundModeSymmetric x y =
  eqRoundMode x y ==> eqRoundMode y x

eqRoundModeTransitive : RoundMode -> RoundMode -> RoundMode -> Bit
property eqRoundModeTransitive x y z =
  (eqRoundMode x y /\ eqRoundMode y z) ==> eqRoundMode x z

// fromRoundingMode, toRoundingMode

fromRoundingMode : RoundingMode -> RoundMode
fromRoundingMode rm =
  if rm == rne then Rnd_Nearest_Even
   | rm == rna then Rnd_Nearest_Away_Zero
   | rm == rtp then Rnd_Plus_Inf
   | rm == rtn then Rnd_Minus_Inf
   | rm == rtz then Rnd_Zero
  else error "Unknown rounding mode"

toRoundingMode : RoundMode -> RoundingMode
toRoundingMode rm =
  case rm of
    Rnd_Nearest_Even      -> rne
    Rnd_Nearest_Away_Zero -> rna
    Rnd_Plus_Inf          -> rtp
    Rnd_Minus_Inf         -> rtn
    Rnd_Zero              -> rtz

fromToRoundingMode : RoundMode -> Bit
property fromToRoundingMode rm = eqRoundMode (fromRoundingMode (toRoundingMode rm)) rm

toFromRoundingMode : RoundingMode -> Bit
property toFromRoundingMode rm =
  validRoundingMode rm ==> toRoundingMode (fromRoundingMode rm) == rm

validRoundingMode : RoundingMode -> Bit
validRoundingMode rm = elem rm [rne, rna, rtp, rtn, rtz]

// Exception

type Exception =
  { invalid_op : Bit
  , divide_0   : Bit
  , overflow   : Bit
  , underflow  : Bit
  , inexact    : Bit
  }

// round

type LimbT  = [64]
type SLimbT = [64]

LIMB_LOG2_BITS = 6
LIMB_BITS = 1 << LIMB_LOG2_BITS

limbMask : [32] -> [32] -> LimbT
limbMask start last_ =
  if n == LIMB_BITS
    then -1
    else ((1 << n) - 1) << start
  where
    n = last_ - start + 1

// Return True if one bit between 0 and bit_pos inclusive is not zero.
scanBitNz :
  {m} (fin m, m >= 1) =>
  [m] -> [m] -> Bit
scanBitNz r bit_pos =
  if bit_pos <$ 0
    then False
    else any (\(i : [m]) -> i <= bit_pos /\ r!i != 0) [0 .. m]

/*
bfGetRndAdd :
  {e, m} (fin m, 64 >= width m) =>
  FloatingPoint e m -> RoundMode -> ([32], [32])
bfGetRndAdd r rnd_mode =
  where
    bit0 = scanBitBz r.sfd
*/

FloatingPoint.icry

:set prover=w4-z3
:load FloatingPoint.cry

// posZero
:prove posZeroFromFloat
:prove posZeroToFloat

// negZero
:prove negZeroFromFloat
:prove negZeroToFloat

// posOne
:prove posOneFromFloat
:prove posOneToFloat

// negOne
:prove negOneFromFloat
:prove negOneToFloat

// posInf
:prove posInfFromFloat
:prove posInfToFloat

// negInf
:prove negInfFromFloat
:prove negInfToFloat

// qnan
:prove qnanFromFloat
:prove qnanToFloat

// isInf
:prove isInfFromFloat
:prove isInfToFloat

// isNaN
:prove isNaNFromFloat
:prove isNaNToFloat

// isSubNormal
:prove isSubNormalFromFloat
:prove isSubNormalToFloat

// isZero
:prove isZeroFromFloat
:prove isZeroToFloat

// eq
:prove eqFromFloat
:prove eqToFloat

// eqBits
:prove eqBitsFromFloat

// eqBitsModuloNaNs
:prove eqBitsModuloNaNsFromFloat
:prove eqBitsModuloNaNsToFloat

// fromBits, toBits
:prove fromToBits
:prove toFromBits

// fromFloat, toFloat
:prove fromToFloat
:prove toFromFloat

// absFP
:prove absFromFloat
:prove absToFloat

// eqRoundMode
:prove eqRoundModeReflexive
:prove eqRoundModeSymmetric
:prove eqRoundModeTransitive

// fromRoundingMode, toRoundingMode
:prove fromToRoundingMode
:prove toFromRoundingMode