diff --git a/Fp/Basic.lean b/Fp/Basic.lean index 4e567a3..73b5ba9 100644 --- a/Fp/Basic.lean +++ b/Fp/Basic.lean @@ -644,6 +644,8 @@ def normalize (uf : UnpackedFloat e s) : UnpackedFloat e s := sig := uf.sig <<< uf.sig.clz } + + @[bv_normalize] def toEUnpackedFloat (uf : UnpackedFloat e s) : EUnpackedFloat e s := EUnpackedFloat.mk .Number uf @@ -793,6 +795,8 @@ def minSubnormE5M2 := PackedFloat.ofBits 5 2 0b00000001#8 /-- info: { state := num, num := + 0x000000004#34 } -/ #guard_msgs in #eval minNormE5M2.toEFixed + + /- - (@bollu's thought): We may like to have `FixedPoint.toRat : FixedPoint → ℚ`, which interprets the FP as a rational. diff --git a/Fp/Division.lean b/Fp/Division.lean index 29b23b4..3dced23 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -268,7 +268,7 @@ def div_on_packedFloat (a b : PackedFloat e s) (mode : RoundingMode) : PackedFlo -- Do division, collapse remainder to a single sticky bit let quot_with_sticky := (dividend / divisor) ++ BitVec.ofBool ((dividend % divisor) ≠ 0) -- Calculate shifts - let divResult := fixedWidthDivideAtPrecision sig_a sig_b (prec := 2 * (s + 1)) (outw := 3 * (s + 1)) + -- let divResult := fixedWidthDivideAtPrecision sig_a sig_b (prec := 2 * (s + 1)) (outw := 3 * (s + 1)) -- let quot_with_sticky := divResult.quotWithSticky let expNumerator := if a.ex > 0 then a.ex - 1 else 0 let expDenominator := if b.ex > 0 then b.ex - 1 else 0 @@ -281,7 +281,7 @@ def div_on_packedFloat (a b : PackedFloat e s) (mode : RoundingMode) : PackedFlo num := { sign -- numerator is larger than denominator, so multiply by the correct amount. - val := quot_with_sticky.setWidth _ <<< (expNumerator - expDenominator) + val := quot_with_sticky.setWidth (2 ^ e + div_len + 1) <<< (expNumerator - expDenominator) hExOffset := by rewrite [Nat.add_lt_add_iff_right] apply Nat.lt_add_left @@ -290,6 +290,7 @@ def div_on_packedFloat (a b : PackedFloat e s) (mode : RoundingMode) : PackedFlo } round _ _ mode quot_lshift else + debug_assert! false let quot_rshift : EFixedPoint (2^e+div_len+1) (2^e+unit_pos+1) := { state := .Number num := { @@ -321,3 +322,135 @@ def div (a b : PackedFloat e s) (m : RoundingMode) : PackedFloat e s := { PackedFloat.getZero _ _ with sign := a.sign ^^ b.sign } else div_on_packedFloat a b m + +theorem div_one_is_id (a : PackedFloat 5 2) + : (div a oneE5M2 .RTZ).equal_denotation a := by + bv_decide + +theorem div_self_is_one (a : PackedFloat 5 2) + (h : ¬a.isNaN ∧ ¬a.isInfinite ∧ ¬a.isZero) + : (div a a .RTZ) = oneE5M2 := by + bv_decide + +/-- x ≥ y ↔ y ≤ x-/ +@[bv_normalize] +def BitVec.sge (x y : BitVec w) : Bool := y.sle x + +@[bv_normalize] +def div_on_unpackedFloat (a b : UnpackedFloat (exponentWidth e s) (s + 1)) (mode : RoundingMode) : PackedFloat e s := + -- a.toRat = (-1)^a.sign * a.sig.toNat * 2 ^ (a.ex.toInt) * 2 ^(-s) + -- a.toRat / b.toRat = (-1)^(a.sign⊕b.sign) * (a.sig.toNat /b.sig.toNat) * (2 ^ (a.ex.toInt) / 2 ^ (b.ex.toInt) *(2 ^(-s) / 2^(-s)) + -- = (-1)^(a.sign⊕b.sign) * (a.sig.toNat /b.sig.toNat) * 2 ^(a.ex.toInt - b.ex.toInt) + -- = (-1)^(a.sign⊕b.sign) * (a.sig.toNat /b.sig.toNat) * 2 ^(a.ex.toInt - b.ex.toInt) + let sign := a.sign ^^ b.sign + let sig_a := a.sig + let sig_b := b.sig + let div_len := 3*(s+1) -- (s + 1) because we add the bit. + let unit_pos := 2*(s+1) + let dividend := (sig_a.setWidth div_len <<< unit_pos) + let divisor := sig_b.setWidth div_len + -- Do division, collapse remainder to a single sticky bit + let quot_with_sticky := (dividend / divisor) ++ BitVec.ofBool ((dividend % divisor) ≠ 0) + -- Calculate shifts + -- let divResult := fixedWidthDivideAtPrecision sig_a sig_b (prec := 2 * (s + 1)) (outw := 3 * (s + 1)) + -- let quot_with_sticky := divResult.quotWithSticky + let expNumerator := a.ex -- if a.ex > 0 then a.ex - 1 else 0 + let expDenominator := b.ex -- if b.ex > 0 then b.ex - 1 else 0 + -- Shift and round + -- | TODO: For the rounding, we still expand out into fixed point. + -- We should instead use the cleverer rounder. + if expNumerator.sge expDenominator then + -- +1 for the sticky bit. + let quot_lshift : EFixedPoint (2^e+div_len+1) (unit_pos+1) := { + state := .Number + num := { + sign + -- numerator is larger than denominator, so multiply by the correct amount. + val := quot_with_sticky.setWidth (2 ^ e + div_len + 1) <<< (expNumerator - expDenominator) + hExOffset := by + rewrite [Nat.add_lt_add_iff_right] + apply Nat.lt_add_left + omega + } + } + round _ _ mode quot_lshift + -- PackedFloat.mk true (BitVec.ofNat _ 0xcafebabe) (BitVec.ofNat _ 0xcafebabe) + else + let quot_rshift : EFixedPoint (2^e+div_len+1) (2^e+unit_pos+1) := { + state := .Number + num := { + sign + -- denominator is larger than numerator, so divide by the correct amount, but first increase the + -- exponent to (2^e), so that the round function rounds correctly. + -- TODO: understand the bounds on our rounding. + -- we shift by '2^e' so we scale up by the same factor as we do with (2^e + unit_pos + 1), + -- such that we represent the same number. + -- 2^(a.ex.toInt - b.ex.toInt) = 2^(-E + E + a.ex.toIn - b.ex.toInt) + -- = 2^-E * 2^(E - b.ex.toInt + a.ex.toInt) [which is ≥ 1 ] + val := ((quot_with_sticky.setWidth _ <<< 2^e) >>> (expDenominator - expNumerator)) + hExOffset := by + rewrite [Nat.add_lt_add_iff_right, Nat.add_lt_add_iff_left] + omega + } + } + round _ _ mode quot_rshift + +/-- +Division of two floating-point numbers, rounded to a floating point number +using the provided rounding mode. +-/ +@[bv_normalize] +def div' (a b : PackedFloat e s) (m : RoundingMode) : PackedFloat e s := + let a := a.unpack + let b := b.unpack + -- a = 0, b = 0 + if a.isNaN ∨ b.isNaN ∨ (a.isInfinite ∧ b.isInfinite) ∨ (a.isZero ∧ b.isZero) then + PackedFloat.getNaN _ _ + -- a = ∞, b = 0 + else if a.isInfinite ∨ b.isZero then + PackedFloat.getInfinity _ _ (a.sign ^^ b.sign) + -- a = *, b = ∞ + else if b.isInfinite then + { PackedFloat.getZero _ _ with sign := a.sign ^^ b.sign } + else if a.isZero then -- note that b must be finite here. + { PackedFloat.getZero _ _ with sign := a.sign ^^ b.sign } + else + div_on_unpackedFloat a.num b.num m + + +/-- Dividing zero by a finite nonzero number gives a zero with the right sign. -/ +theorem zero_div_is_zero (a b : PackedFloat 5 2) (ha : a.isZero) (hb : b.isNormOrSubnorm) (hb : ¬b.isZero) + : (div a b .RTZ).isZero ∧ (div a b .RTZ).sign = a.sign ^^ b.sign := by + bv_decide + +def cex' : PackedFloat 5 2 where + ex := 0#5 + sig := 2#2 + sign := false + +#check EUnpackedFloat.toRat? +#eval cex'.toRat? -- 2^-15 +#eval cex'.unpack.toRat? -- 2^-15 +-- 0xfff2#16 = * 2 +/-- info: -14 -/ +#guard_msgs in #eval cex'.unpack.num.ex.toInt +/-- info: 0 -/ +#guard_msgs in #eval oneE5M2.unpack.num.ex.toInt + +/-- info: { sign := +, ex := 0x00#5, sig := 0x2#2 } -/ +#guard_msgs in #eval div_on_unpackedFloat cex'.unpack.num oneE5M2.unpack.num .RTZ + +set_option debugAssertions true in +theorem div_one_is_id_on_numerator_ex_larger (a : PackedFloat 5 2) (h : a.unpack.num.ex.sge oneE5M2.unpack.num.ex) + : (div' a oneE5M2 .RTZ).equal_denotation a := by + bv_decide + +set_option debugAssertions true in +theorem div_one_is_id_on_numerator_ex_smaller (a : PackedFloat 5 2) (h : ¬ a.unpack.num.ex.sge oneE5M2.unpack.num.ex) + : (div' a oneE5M2 .RTZ).equal_denotation a := by + bv_decide + +theorem div_self_is_one' (a : PackedFloat 5 2) + (h : ¬a.isNaN ∧ ¬a.isInfinite ∧ ¬a.isZero) + : (div' a a .RTZ) = oneE5M2 := by + bv_decide diff --git a/Fp/Packing.lean b/Fp/Packing.lean index 36e597c..31be504 100644 --- a/Fp/Packing.lean +++ b/Fp/Packing.lean @@ -33,14 +33,14 @@ def EUnpackedFloat.pack (uf : EUnpackedFloat (exponentWidth e s) (s + 1)) ex := bif uf.isNaN || uf.isInfinite then BitVec.allOnes e else if uf.isZero || !inNormalRange then - BitVec.zero _ + (0#_) else -- bif uf.isNorm then -- Truncate msbs used to normalize subnormals (uf.exp + BitVec.ofNat _ (bias e)).truncate _ sig := bif uf.isNaN then BitVec.ofNat _ (2 ^ (s - 1)) else bif uf.isInfinite || uf.isZero then - BitVec.zero _ + (0#_) else bif inNormalRange then uf.sig.truncate _ else -- bif uf.isSubnorm then @@ -48,6 +48,18 @@ def EUnpackedFloat.pack (uf : EUnpackedFloat (exponentWidth e s) (s + 1)) (uf.sig >>> shift).truncate _ } +theorem PackedFloat.isInfinite_pack_unpack (pf : PackedFloat 5 10) (hpf : pf.unpack.isInfinite) : + pf.unpack.pack.isInfinite = pf.isInfinite ∧ pf.unpack.pack.sign = pf.sign := by + bv_decide + +theorem PackedFloat.isNaN_pack_unpack (pf : PackedFloat 5 10) : + pf.unpack.pack.isNaN = pf.isNaN := by + bv_decide + +theorem PackedFloat.pack_unpack (pf : PackedFloat 5 10) (hpf : pf.isNormOrSubnorm) : + pf.unpack.pack = pf := by + bv_decide + /-- info: { sign := +, ex := 0x00#5, sig := 0x001#10 } -/ #guard_msgs in #eval { sign := false, ex := 0x00#5, sig := 0x001#10 : PackedFloat _ _ }.unpack.pack /-- info: { sign := +, ex := 0x0f#5, sig := 0x3ff#10 } -/