From 41f2c4493d24198937f1aa3e6fe4238abb9c7dfa Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Fri, 19 Dec 2025 18:35:10 +0000 Subject: [PATCH 1/9] chore: new rounder --- Fp/Basic.lean | 10 ++++++++++ Fp/Division.lean | 36 ++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/Fp/Basic.lean b/Fp/Basic.lean index 3469fc0..0b92107 100644 --- a/Fp/Basic.lean +++ b/Fp/Basic.lean @@ -805,10 +805,20 @@ def Nat.ceilLog2 (n : Nat) : Nat := def bias (e : Nat) : Nat := 2 ^ (e - 1) - 1 +/-- The minimum value the exponent can take when unbiased. -/ @[bv_normalize] def minNormalExp (e : Nat) : Int := -(bias e - 1) +/-- The max value the exponent can take when unbiased. -/ +@[bv_normalize] +def maxNormalExp (e : Nat) : Int := (bias e) + +/-- The value the subnormal exponent can take. -/ +@[bv_normalize] +def subnormalExp (e : Nat) : Int := + minNormalExp e - 1 + -- This is a simpler (but less tight) bound than `exponentWidth`. -- It's logarithmically larger. @[bv_normalize] diff --git a/Fp/Division.lean b/Fp/Division.lean index ecdcf7d..63fd6e8 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -259,6 +259,42 @@ def BitVec.monus (a : BitVec w) (b : BitVec w) : BitVec w := @[bv_normalize] def BitVec.sge (x y : BitVec w) : Bool := y.sle x +#check round + +-- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) +-- choose e = 0, s = 1 then see that after normalization, we e. +-- If our exponent is +def roundFast (up : UnpackedFloat e s) : EUnpackedFloat e s := + let up := up.normalize + if up.sig = 0 + then EUnpackedFloat.mkNumber up + else + -- we are in the <1.-----> case. + -- let exOffset' := 2^(exWidth - 1) + sigWidth - 2 + -- trim bitvector + -- can absorb 2^s component into the significand, but any more + -- cannot be absorbed. + if up.ex > BitVec.ofInt e (maxNormalExp e) + BitVec.ofInt e s + then + EUnpackedFloat.mkInfinity up.sign + else if up.ex < BitVec.ofInt e (minNormalExp e) - BitVec.ofInt e s + then + -- | too small to be subnormal. + EUnpackedFloat.mkZero up.sign + else if up.ex < BitVec.ofInt e (minNormalExp e) + then + -- | subnormal case. + let shift := (BitVec.ofInt e (minNormalExp e) - up.ex).toNat + let sig_shifted := up.sig >>> shift + sorry + -- let round_result := roundSignificandAndStickyBit up.sign sig_shifted s RoundingMode.RTZ + -- EUnpackedFloat.mkSubnormal round_result + else + -- | normal case. + -- let round_result := roundSignificandAndStickyBit up.sign up.sig s RoundingMode.RTZ + -- EUnpackedFloat.mkNormal (up.ex) round_result + sorry + @[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) From 509d812d55a7a2031dc2affa2e377f99db05db23 Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Tue, 23 Dec 2025 12:14:58 +0100 Subject: [PATCH 2/9] chore: add new notes --- Fp/Division.lean | 85 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 72 insertions(+), 13 deletions(-) diff --git a/Fp/Division.lean b/Fp/Division.lean index 63fd6e8..d2a7b2d 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -261,39 +261,98 @@ def BitVec.sge (x y : BitVec w) : Bool := y.sle x #check round + +/- + +The correct way to think of an UnpackedFloat is as a fixed point number, +written in scientific notation. + +So, consider fixed point numbers with three digits, dot right after the first digit (as in normal floating point). +This gives us: + +→ 0.00 0.01 0.10 0.11 +→ 1.00 1.01 1.10 1.11 + +→ 0.00 = _ | 0.01 = 1.00 * 10^(-2) | 0.10 = 1.00 * 10^(-1) | 0.11 = 1.10 * 10^(-1) +→ 1.00 = 1.00 * 10^(0) | 1.01 = 1.01 * 10^(0) | 1.10 = 1.10 * 10^(0) | 1.11 = 1.11 * 10^(0) + +Next, to encode the `1.00`, we write this as: + +→ 0.00 = _ | 0.01 = (100 * 10^(-2)) * 10^(-2) | 0.10 = (100 * 10^(-2)) * 10^(-1) | 0.11 = (110 * 10^(-2)) * 10^(-1) +→ 1.00 = (100 * 10^(-2)) * 10^(0) | 1.01 = (101 * 10^(-2)) * 10^(0) | 1.10 = (110*10^(-2)) * 10^(0) | 1.11 = (111 * 10^(-2)) * 10^(0) + +In this way, we can see our number as a composition of scientific notation with fixed-width exponent. +In the UnpackedFloat, we only store the exponent, not the constant fixed-point shift (10^(-2)), +since this is common to all numbers, and also, morally, it is absorbed into the significand as a fixed-point. +-/ + -- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) -- choose e = 0, s = 1 then see that after normalization, we e. -- If our exponent is -def roundFast (up : UnpackedFloat e s) : EUnpackedFloat e s := +def roundFast (up : UnpackedFloat e s) : EUnpackedFloat e' s' := let up := up.normalize if up.sig = 0 - then EUnpackedFloat.mkNumber up + then EUnpackedFloat.mkNumber { + sign := up.sign + ex := 0#e' + sig := 0#s' + } else -- we are in the <1.-----> case. -- let exOffset' := 2^(exWidth - 1) + sigWidth - 2 -- trim bitvector -- can absorb 2^s component into the significand, but any more -- cannot be absorbed. - if up.ex > BitVec.ofInt e (maxNormalExp e) + BitVec.ofInt e s - then + if hExTooBig : up.ex > BitVec.ofInt e (maxNormalExp e') + BitVec.ofInt e s' then EUnpackedFloat.mkInfinity up.sign - else if up.ex < BitVec.ofInt e (minNormalExp e) - BitVec.ofInt e s - then + else if hExTooSmall : up.ex < BitVec.ofInt e (minNormalExp e') - BitVec.ofInt e s' then -- | too small to be subnormal. EUnpackedFloat.mkZero up.sign - else if up.ex < BitVec.ofInt e (minNormalExp e) - then - -- | subnormal case. - let shift := (BitVec.ofInt e (minNormalExp e) - up.ex).toNat - let sig_shifted := up.sig >>> shift - sorry + -- else if up.ex < BitVec.ofInt e (minNormalExp e) then + -- -- | subnormal case. + -- let shift := (BitVec.ofInt e (minNormalExp e) - up.ex).toNat + -- let sig_shifted := up.sig >>> shift + -- sorry -- let round_result := roundSignificandAndStickyBit up.sign sig_shifted s RoundingMode.RTZ -- EUnpackedFloat.mkSubnormal round_result else + let adjustedExp := + if up.ex > BitVec.ofInt e (maxNormalExp e') then + BitVec.ofInt e (maxNormalExp e) + else if up.ex < BitVec.ofInt e (minNormalExp e') then + BitVec.ofInt e (minNormalExp e') + else + up.ex + let truncatedSig : BitVec s' := + if adjustedExp = up.ex then + up.sig.truncate _ + else if adjustedExp < up.ex then + -- we decreased exponent, so increase significand + (up.sig <<< (up.ex - adjustedExp)).truncate _ + else + -- we increased exponent, so decrease significand + let shift := (adjustedExp - up.ex) + (up.sig >>> shift).truncate _ + let remainder : BitVec (s - s') := + if adjustedExp = up.ex then + up.sig.drop s' + else if adjustedExp < up.ex then + (up.sig <<< (up.ex - adjustedExp)).drop s' + else + let shift := (adjustedExp - up.ex) + (up.sig >>> shift).drop s' + sorry + + + -- I need to ensure that exponent is in + -- [minNormalExp-1, maxNormalExp] + -- normal plus subnormal case. + -- drop the number of bits that can fit into the exponent. + + -- let shouldRoundAway := shouldRoundAway .RNE -- | normal case. -- let round_result := roundSignificandAndStickyBit up.sign up.sig s RoundingMode.RTZ -- EUnpackedFloat.mkNormal (up.ex) round_result - sorry @[bv_normalize] def div_on_unpackedFloat (a b : UnpackedFloat (exponentWidth e s) (s + 1)) (mode : RoundingMode) : PackedFloat e s := From a669dbd905b8a2e5adb05707ac1c0615c6813c19 Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 08:35:23 +0100 Subject: [PATCH 3/9] chore: write about guard and sticky --- Fp/Division.lean | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/Fp/Division.lean b/Fp/Division.lean index d2a7b2d..d603670 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -262,6 +262,21 @@ def BitVec.sge (x y : BitVec w) : Bool := y.sle x #check round +def BitVec.dropLsbs (bv : BitVec w) (n : Nat) : BitVec (w - n) := + (bv >>> n).setWidth _ + +theorem getMsbD_dropLsbs {w n : Nat} (h : n < w) (bv : BitVec w) : + bv.getMsbD i = (bv.dropLsbs n).getMsbD i := by + simp [BitVec.dropLsbs, BitVec.getMsbD] + by_cases hi : i < w + · simp [hi] + sorry + · simp [hi] + sorry + +def BitVec.dropMsb (bv : BitVec w) (n : Nat) : BitVec (w - n) := + bv.setWidth _ + /- The correct way to think of an UnpackedFloat is as a fixed point number, @@ -286,10 +301,32 @@ In the UnpackedFloat, we only store the exponent, not the constant fixed-point s since this is common to all numbers, and also, morally, it is absorbed into the significand as a fixed-point. -/ +/- +#### Why guard and sticky bits: + +Suppose we want to round to 0 bits of precision, using round to nearest even. Then, if we have: + +- 10.1 -> 10 +- 10.5???, we need to know if the ??? is zero or not. + + 10.5000 -> 10 + * If it is zero, then we round to the nearest even, which is 10. + + 10.5001 -> 11 + * If it is non-zero, we *always* round up to 11, since it's strictly closer to 11 than 10. +- 10.9 -> 11 +- 11.5????? + + 11.5000 -> 12 + * If it is zero, then we round to the nearest even, which is 12. + + 11.5001 -> 12 + * If it is non-zero, we *always* round up to 12, since it's strictly closer to 12 than 11. + +Thus, see that just having one sticky bit is enough, since it tells us whether we are at the exact halfway point or not, and that's all we need to know. + +-/ + -- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) -- choose e = 0, s = 1 then see that after normalization, we e. -- If our exponent is -def roundFast (up : UnpackedFloat e s) : EUnpackedFloat e' s' := +def roundRNEFast (up : UnpackedFloat e s) : EUnpackedFloat e' s' := let up := up.normalize if up.sig = 0 then EUnpackedFloat.mkNumber { From 1f53fe54d0013c7577fd8fe89bb3b8950b20e84e Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 09:02:48 +0100 Subject: [PATCH 4/9] chore: add rounder --- Fp/Division.lean | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/Fp/Division.lean b/Fp/Division.lean index d603670..5ed194f 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -321,6 +321,11 @@ Suppose we want to round to 0 bits of precision, using round to nearest even. Th Thus, see that just having one sticky bit is enough, since it tells us whether we are at the exact halfway point or not, and that's all we need to know. + +We will return an EUnpackedFloat, +which can signal NaN/Infinity/Zero, in addition to normal numbers. + + -/ -- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) @@ -342,7 +347,7 @@ def roundRNEFast (up : UnpackedFloat e s) : EUnpackedFloat e' s' := -- cannot be absorbed. if hExTooBig : up.ex > BitVec.ofInt e (maxNormalExp e') + BitVec.ofInt e s' then EUnpackedFloat.mkInfinity up.sign - else if hExTooSmall : up.ex < BitVec.ofInt e (minNormalExp e') - BitVec.ofInt e s' then + else if hExTooSmall : up.ex < BitVec.ofInt e (minNormalExp e') - BitVec.ofInt e s' then -- | too small to be subnormal. EUnpackedFloat.mkZero up.sign -- else if up.ex < BitVec.ofInt e (minNormalExp e) then @@ -372,13 +377,20 @@ def roundRNEFast (up : UnpackedFloat e s) : EUnpackedFloat e' s' := (up.sig >>> shift).truncate _ let remainder : BitVec (s - s') := if adjustedExp = up.ex then - up.sig.drop s' + up.sig.dropMsb s' else if adjustedExp < up.ex then - (up.sig <<< (up.ex - adjustedExp)).drop s' + (up.sig <<< (up.ex - adjustedExp)).dropMsb s' else let shift := (adjustedExp - up.ex) - (up.sig >>> shift).drop s' - sorry + (up.sig >>> shift).dropMsb s' + let guardBit := remainder.getMsbD 0 + let stickyBit : Bool := (remainder.dropMsb 1) ≠ 0 + let isOdd := truncatedSig.getMsbD 0 + let shouldRoundAway : Bool := guardBit && (stickyBit || isOdd) + if shouldRoundAway then + sorry + else + sorry -- I need to ensure that exponent is in From d30cdc81a6400b7feffb907abc8253c87bddc365 Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 09:59:09 +0100 Subject: [PATCH 5/9] chore: implement rounder, but this will need refactoring / fixing We also develop a bunch of infra to allow us to compute on BVs without having to worry about extension, by writing custom functions which extend appropriately. Their `toInt` lemmas will not have modulos, which will do the right thing when performing arithemtic. --- Fp/Division.lean | 172 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 124 insertions(+), 48 deletions(-) diff --git a/Fp/Division.lean b/Fp/Division.lean index 5ed194f..fdf0ffb 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -255,6 +255,37 @@ def BitVec.monus (a : BitVec w) (b : BitVec w) : BitVec w := if a ≤ b then 0#w else a - b +def BitVec.addExtending (a : BitVec v) (b : BitVec w) : BitVec (max v w + 1) := + let a' := a.signExtend (max v w + 1) + let b' := b.signExtend (max v w + 1) + a' + b' + +def BitVec.subExtending (a : BitVec v) (b : BitVec w) : BitVec (max v w + 1) := + let a' := a.signExtend (max v w + 1) + let b' := b.signExtend (max v w + 1) + a' - b' + +def BitVec.eqExtending (a : BitVec v) (b : BitVec w) : Prop := + let a' := a.signExtend (max v w) + let b' := b.signExtend (max v w) + a' = b' + +def BitVec.sltExtending (a : BitVec v) (b : BitVec w) : Prop := + let a' := a.signExtend (max v w) + let b' := b.signExtend (max v w) + a'.slt b' + + +instance : Decidable (BitVec.eqExtending a b) := by + unfold BitVec.eqExtending + simp + infer_instance + +instance : Decidable (BitVec.sltExtending a b) := by + unfold BitVec.sltExtending + simp + infer_instance + /-- x ≥ y ↔ y ≤ x-/ @[bv_normalize] def BitVec.sge (x y : BitVec w) : Bool := y.sle x @@ -277,6 +308,33 @@ theorem getMsbD_dropLsbs {w n : Nat} (h : n < w) (bv : BitVec w) : def BitVec.dropMsb (bv : BitVec w) (n : Nat) : BitVec (w - n) := bv.setWidth _ +def BitVec.takeMsb (bv : BitVec w) (n : Nat) : BitVec n := + (bv >>> (w - n)).setWidth _ + +def BitVec.splitAtMsbs (bv : BitVec w) (n : Nat) : BitVec n × BitVec (w - n) := + (bv.takeMsb n, bv.dropMsb n) + + +def EUnpackedFloat.incrSignificand {e s : Nat} (eu : EUnpackedFloat e s) : EUnpackedFloat e s := + match eu.state with + | .NaN => eu + | .Infinity => eu + | .Number => + if eu.num.sig = BitVec.allOnes _ then + -- we overflowed the significand, so we need to adjust exponent. + let newExp := eu.num.ex + 1 + EUnpackedFloat.mkNumber { + sign := eu.num.sign + ex := newExp + sig := 1#s -- since we overflowed, the new significand is 1.000... (is that right?) + } + else + EUnpackedFloat.mkNumber { + sign := eu.num.sign + ex := eu.num.ex + sig := eu.num.sig + 1 + } + /- The correct way to think of an UnpackedFloat is as a fixed point number, @@ -326,71 +384,89 @@ We will return an EUnpackedFloat, which can signal NaN/Infinity/Zero, in addition to normal numbers. +Rounding normalized numbers +--------------------------- + +Suppose we have 3 digits of precision, and the exponent can go from -9 to 7. +Then, suppose we have an unnormalized number: + +##### Smallest Exponent +- 0.001 * 10^7 + This upon normalization becomes: 1.000 * 10^(7-3) = 1.000 * 10^4 + +- 0.001 * 10^(-9) + This upon normalization becomes: 1.000 * 10^(-9-3) = 1.000 * 10^(-12) + + +So, even though our scale can go from -9 to 7, after normalization, we can have exponents from -12 to 7. +If our exponent is smaller than -12, then we cannot fit it. + +##### Largest Exponent + +- 1.111 * 10^7 | This upon normalization stays the same: 1.111 * 10^7. -/ -- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) -- choose e = 0, s = 1 then see that after normalization, we e. -- If our exponent is -def roundRNEFast (up : UnpackedFloat e s) : EUnpackedFloat e' s' := - let up := up.normalize - if up.sig = 0 - then EUnpackedFloat.mkNumber { - sign := up.sign - ex := 0#e' - sig := 0#s' - } +def roundRNEFast (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := + let inUf := inUf.normalize + -- Great, we have a number of the form <1.xxxxx> * 2^(ex) or, 0 * 2^(ex). + if inUf.sig = 0 + -- we are in the <0> case. + then EUnpackedFloat.mkNumber { + sign := inUf.sign + ex := 0#e' + sig := 0#s' + } else -- we are in the <1.-----> case. - -- let exOffset' := 2^(exWidth - 1) + sigWidth - 2 - -- trim bitvector - -- can absorb 2^s component into the significand, but any more - -- cannot be absorbed. - if hExTooBig : up.ex > BitVec.ofInt e (maxNormalExp e') + BitVec.ofInt e s' then - EUnpackedFloat.mkInfinity up.sign - else if hExTooSmall : up.ex < BitVec.ofInt e (minNormalExp e') - BitVec.ofInt e s' then - -- | too small to be subnormal. - EUnpackedFloat.mkZero up.sign - -- else if up.ex < BitVec.ofInt e (minNormalExp e) then - -- -- | subnormal case. - -- let shift := (BitVec.ofInt e (minNormalExp e) - up.ex).toNat - -- let sig_shifted := up.sig >>> shift - -- sorry - -- let round_result := roundSignificandAndStickyBit up.sign sig_shifted s RoundingMode.RTZ - -- EUnpackedFloat.mkSubnormal round_result + if hExTooBig : inUf.ex > BitVec.ofInt e (maxNormalExp e') then -- + BitVec.ofInt e s' then + EUnpackedFloat.mkInfinity inUf.sign + -- | TODO: why does this computation fit? In particular, why does adding `s'` not overflow? + -- We may need an extension here. + else if hExTooSmall : (inUf.ex).slt (BitVec.ofInt e (minNormalExp e') + BitVec.ofInt e s') then + -- See example, where `0.001 * 10^-5 -> 1.00 * 10^-7`, which is the smallest exponent we can use. + EUnpackedFloat.mkZero inUf.sign else - let adjustedExp := - if up.ex > BitVec.ofInt e (maxNormalExp e') then - BitVec.ofInt e (maxNormalExp e) - else if up.ex < BitVec.ofInt e (minNormalExp e') then - BitVec.ofInt e (minNormalExp e') + -- Adjust the exponent to be in range. + let adjustedExp : BitVec e' := + if inUf.ex > BitVec.ofInt e (maxNormalExp e') then + BitVec.ofInt e' (maxNormalExp e') + else if inUf.ex < BitVec.ofInt e (minNormalExp e') then + BitVec.ofInt e' (minNormalExp e') else - up.ex - let truncatedSig : BitVec s' := - if adjustedExp = up.ex then - up.sig.truncate _ - else if adjustedExp < up.ex then + inUf.ex.signExtend e' + -- Adjust the significand. + let (truncatedSig, remainder) : BitVec s' × BitVec (s - s') := + if BitVec.eqExtending adjustedExp inUf.ex then + (inUf.sig.splitAtMsbs s') + else if BitVec.sltExtending adjustedExp inUf.ex then -- we decreased exponent, so increase significand - (up.sig <<< (up.ex - adjustedExp)).truncate _ + let sig' := (inUf.sig <<< (BitVec.subExtending inUf.ex adjustedExp)) + sig'.splitAtMsbs s' else -- we increased exponent, so decrease significand - let shift := (adjustedExp - up.ex) - (up.sig >>> shift).truncate _ - let remainder : BitVec (s - s') := - if adjustedExp = up.ex then - up.sig.dropMsb s' - else if adjustedExp < up.ex then - (up.sig <<< (up.ex - adjustedExp)).dropMsb s' - else - let shift := (adjustedExp - up.ex) - (up.sig >>> shift).dropMsb s' + let shift := BitVec.subExtending adjustedExp inUf.ex + let sig' := (inUf.sig >>> shift).truncate _ + sig'.splitAtMsbs s' let guardBit := remainder.getMsbD 0 let stickyBit : Bool := (remainder.dropMsb 1) ≠ 0 - let isOdd := truncatedSig.getMsbD 0 + let isOdd := truncatedSig.getLsbD 0 let shouldRoundAway : Bool := guardBit && (stickyBit || isOdd) + let ufOut : UnpackedFloat e' s' := + { + sign := inUf.sign + ex := adjustedExp + sig := truncatedSig + } + let eufOut : EUnpackedFloat e' s' := + EUnpackedFloat.mkNumber ufOut if shouldRoundAway then - sorry + -- we should round up. + eufOut.incrSignificand else - sorry + eufOut -- I need to ensure that exponent is in From 1fccc9cdb91f00757f3fa203f6fc30129c49d75c Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 10:32:05 +0100 Subject: [PATCH 6/9] chore: think about over/underflow in rounding --- Fp/Division.lean | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/Fp/Division.lean b/Fp/Division.lean index fdf0ffb..73d1267 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -296,6 +296,7 @@ def BitVec.sge (x y : BitVec w) : Bool := y.sle x def BitVec.dropLsbs (bv : BitVec w) (n : Nat) : BitVec (w - n) := (bv >>> n).setWidth _ +@[bv_normalize] theorem getMsbD_dropLsbs {w n : Nat} (h : n < w) (bv : BitVec w) : bv.getMsbD i = (bv.dropLsbs n).getMsbD i := by simp [BitVec.dropLsbs, BitVec.getMsbD] @@ -305,15 +306,27 @@ theorem getMsbD_dropLsbs {w n : Nat} (h : n < w) (bv : BitVec w) : · simp [hi] sorry +@[bv_normalize] def BitVec.dropMsb (bv : BitVec w) (n : Nat) : BitVec (w - n) := bv.setWidth _ +@[bv_normalize] def BitVec.takeMsb (bv : BitVec w) (n : Nat) : BitVec n := (bv >>> (w - n)).setWidth _ +@[bv_normalize] def BitVec.splitAtMsbs (bv : BitVec w) (n : Nat) : BitVec n × BitVec (w - n) := (bv.takeMsb n, bv.dropMsb n) +/-- Shift left 'bv' by 'shAmt', extending 'bv' to width 'v' first. -/ +@[bv_normalize] +def BitVec.shlExtending (bv : BitVec w) (shAmt : BitVec v) : BitVec v := + (bv.zeroExtend v) <<< shAmt + +@[bv_normalize] +def BitVec.shrExtending (bv : BitVec w) (shAmt : BitVec v) : BitVec v := + (bv.zeroExtend v) >>> shAmt + def EUnpackedFloat.incrSignificand {e s : Nat} (eu : EUnpackedFloat e s) : EUnpackedFloat e s := match eu.state with @@ -441,12 +454,16 @@ def roundRNEFast (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := let (truncatedSig, remainder) : BitVec s' × BitVec (s - s') := if BitVec.eqExtending adjustedExp inUf.ex then (inUf.sig.splitAtMsbs s') + -- e' < e else if BitVec.sltExtending adjustedExp inUf.ex then -- we decreased exponent, so increase significand + -- | TODO: can this overflow? let sig' := (inUf.sig <<< (BitVec.subExtending inUf.ex adjustedExp)) sig'.splitAtMsbs s' else + -- e' > e -- we increased exponent, so decrease significand + -- | TODO: can this overflow? let shift := BitVec.subExtending adjustedExp inUf.ex let sig' := (inUf.sig >>> shift).truncate _ sig'.splitAtMsbs s' From 24c6eabfb997b2bd320f66f451ae96e12b90a1d3 Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 10:49:14 +0100 Subject: [PATCH 7/9] chore: make bv_decide internal panic info: Fp/Division.lean:600:0: PANIC at List.get!Internal Init.GetElem:328:18: invalid index backtrace: 0 libleanshared.dylib 0x00000001149692fc _ZN4leanL15lean_panic_implEPKcmb + 268 1 libleanshared.dylib 0x00000001149697d0 lean_panic_fn + 40 2 libleanshared.dylib 0x0000000110276ed4 l___private_Lean_Elab_Tactic_BVDecide_Frontend_BVDecide_0__Lean_Elab_Tactic_BVDecide_Frontend_DiagnosisM_diagnose_transformEquation + 7064 3 libleanshared.dylib 0x000000011027bf34 l___private_Init_Data_Array_Basic_0__Array_forIn_x27Unsafe_loop___at___00Lean_Elab_Tactic_BVDecide_Frontend_DiagnosisM_diagnose_spec__5 + 416 4 libleanshared.dylib 0x000000011027cbe0 l_Lean_Elab_Tactic_BVDecide_Frontend_DiagnosisM_diagnose + 88 5 libleanshared.dylib 0x000000011497e2f0 lean_apply_7 + 184 6 libleanshared.dylib 0x000000011026e9bc l_Lean_Elab_Tactic_BVDecide_Frontend_DiagnosisM_run___lam__0 + 116 7 libleanshared.dylib 0x000000011497bf34 lean_apply_5 + 1020 8 libleanshared.dylib 0x000000011026e810 l_Lean_MVarId_withContext___at___00Lean_Elab_Tactic_BVDecide_Frontend_DiagnosisM_run_spec__0___redArg + 44 9 libleanshared.dylib 0x000000011027e104 l_Lean_Elab_Tactic_BVDecide_Frontend_explainCounterExampleQuality + 220 --- Fp/Division.lean | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/Fp/Division.lean b/Fp/Division.lean index 73d1267..6481fc9 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -422,7 +422,7 @@ If our exponent is smaller than -12, then we cannot fit it. -- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) -- choose e = 0, s = 1 then see that after normalization, we e. -- If our exponent is -def roundRNEFast (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := +def roundRNEFastUF (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := let inUf := inUf.normalize -- Great, we have a number of the form <1.xxxxx> * 2^(ex) or, 0 * 2^(ex). if inUf.sig = 0 @@ -486,16 +486,6 @@ def roundRNEFast (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := eufOut - -- I need to ensure that exponent is in - -- [minNormalExp-1, maxNormalExp] - -- normal plus subnormal case. - -- drop the number of bits that can fit into the exponent. - - -- let shouldRoundAway := shouldRoundAway .RNE - -- | normal case. - -- let round_result := roundSignificandAndStickyBit up.sign up.sig s RoundingMode.RTZ - -- EUnpackedFloat.mkNormal (up.ex) round_result - @[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) @@ -517,6 +507,14 @@ def div_on_unpackedFloat (a b : UnpackedFloat (exponentWidth e s) (s + 1)) (mode 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 + let ufIn : UnpackedFloat _ _ := { + sign + ex := BitVec.subExtending expNumerator expDenominator + sig := quot_with_sticky + } + roundRNEFastUF ufIn |>.pack + + /- -- | TODO: For the rounding, we still expand out into fixed point. -- We should instead use the cleverer rounder. if expNumerator.sge expDenominator then @@ -554,6 +552,7 @@ def div_on_unpackedFloat (a b : UnpackedFloat (exponentWidth e s) (s + 1)) (mode } } round _ _ mode quot_rshift + -/ /-- Division of two floating-point numbers, rounded to a floating point number From b38961a1d9c692638842d950334cb5f4ef58123c Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 11:06:41 +0100 Subject: [PATCH 8/9] chore: add bv_normalize calls --- Fp/Division.lean | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Fp/Division.lean b/Fp/Division.lean index 6481fc9..c60e7e8 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -328,6 +328,7 @@ def BitVec.shrExtending (bv : BitVec w) (shAmt : BitVec v) : BitVec v := (bv.zeroExtend v) >>> shAmt +@[bv_normalize] def EUnpackedFloat.incrSignificand {e s : Nat} (eu : EUnpackedFloat e s) : EUnpackedFloat e s := match eu.state with | .NaN => eu @@ -422,6 +423,7 @@ If our exponent is smaller than -12, then we cannot fit it. -- the number is sig.toNat * 2 ^(-sigPrec) * 2 ^ (ex.toInt) -- choose e = 0, s = 1 then see that after normalization, we e. -- If our exponent is +@[bv_normalize] def roundRNEFastUF (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := let inUf := inUf.normalize -- Great, we have a number of the form <1.xxxxx> * 2^(ex) or, 0 * 2^(ex). From 8e4cadccc921dc20928f838837c82aa887c0efdc Mon Sep 17 00:00:00 2001 From: Siddharth Bhat Date: Mon, 5 Jan 2026 11:45:07 +0100 Subject: [PATCH 9/9] chore: having widths that are not-constant causes fuckups at type level --- Fp/Division.lean | 23 +++++++++++++++-------- Fp/Tactics/Simps.lean | 6 ++++++ 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/Fp/Division.lean b/Fp/Division.lean index c60e7e8..4ac247d 100644 --- a/Fp/Division.lean +++ b/Fp/Division.lean @@ -255,21 +255,25 @@ def BitVec.monus (a : BitVec w) (b : BitVec w) : BitVec w := if a ≤ b then 0#w else a - b +@[bv_normalize] def BitVec.addExtending (a : BitVec v) (b : BitVec w) : BitVec (max v w + 1) := let a' := a.signExtend (max v w + 1) let b' := b.signExtend (max v w + 1) a' + b' +@[bv_normalize] def BitVec.subExtending (a : BitVec v) (b : BitVec w) : BitVec (max v w + 1) := let a' := a.signExtend (max v w + 1) let b' := b.signExtend (max v w + 1) a' - b' +@[bv_normalize] def BitVec.eqExtending (a : BitVec v) (b : BitVec w) : Prop := let a' := a.signExtend (max v w) let b' := b.signExtend (max v w) a' = b' +@[bv_normalize] def BitVec.sltExtending (a : BitVec v) (b : BitVec w) : Prop := let a' := a.signExtend (max v w) let b' := b.signExtend (max v w) @@ -278,12 +282,12 @@ def BitVec.sltExtending (a : BitVec v) (b : BitVec w) : Prop := instance : Decidable (BitVec.eqExtending a b) := by unfold BitVec.eqExtending - simp + simp infer_instance instance : Decidable (BitVec.sltExtending a b) := by unfold BitVec.sltExtending - simp + simp infer_instance /-- x ≥ y ↔ y ≤ x-/ @@ -293,6 +297,7 @@ def BitVec.sge (x y : BitVec w) : Bool := y.sle x #check round +@[bv_normalize] def BitVec.dropLsbs (bv : BitVec w) (n : Nat) : BitVec (w - n) := (bv >>> n).setWidth _ @@ -300,7 +305,7 @@ def BitVec.dropLsbs (bv : BitVec w) (n : Nat) : BitVec (w - n) := theorem getMsbD_dropLsbs {w n : Nat} (h : n < w) (bv : BitVec w) : bv.getMsbD i = (bv.dropLsbs n).getMsbD i := by simp [BitVec.dropLsbs, BitVec.getMsbD] - by_cases hi : i < w + by_cases hi : i < w · simp [hi] sorry · simp [hi] @@ -357,7 +362,7 @@ written in scientific notation. So, consider fixed point numbers with three digits, dot right after the first digit (as in normal floating point). This gives us: -→ 0.00 0.01 0.10 0.11 +→ 0.00 0.01 0.10 0.11 → 1.00 1.01 1.10 1.11 → 0.00 = _ | 0.01 = 1.00 * 10^(-2) | 0.10 = 1.00 * 10^(-1) | 0.11 = 1.10 * 10^(-1) @@ -369,7 +374,7 @@ Next, to encode the `1.00`, we write this as: → 1.00 = (100 * 10^(-2)) * 10^(0) | 1.01 = (101 * 10^(-2)) * 10^(0) | 1.10 = (110*10^(-2)) * 10^(0) | 1.11 = (111 * 10^(-2)) * 10^(0) In this way, we can see our number as a composition of scientific notation with fixed-width exponent. -In the UnpackedFloat, we only store the exponent, not the constant fixed-point shift (10^(-2)), +In the UnpackedFloat, we only store the exponent, not the constant fixed-point shift (10^(-2)), since this is common to all numbers, and also, morally, it is absorbed into the significand as a fixed-point. -/ @@ -385,7 +390,7 @@ Suppose we want to round to 0 bits of precision, using round to nearest even. Th + 10.5001 -> 11 * If it is non-zero, we *always* round up to 11, since it's strictly closer to 11 than 10. - 10.9 -> 11 -- 11.5????? +- 11.5????? + 11.5000 -> 12 * If it is zero, then we round to the nearest even, which is 12. + 11.5001 -> 12 @@ -405,7 +410,7 @@ Suppose we have 3 digits of precision, and the exponent can go from -9 to 7. Then, suppose we have an unnormalized number: ##### Smallest Exponent -- 0.001 * 10^7 +- 0.001 * 10^7 This upon normalization becomes: 1.000 * 10^(7-3) = 1.000 * 10^4 - 0.001 * 10^(-9) @@ -484,7 +489,7 @@ def roundRNEFastUF (inUf : UnpackedFloat e s) : EUnpackedFloat e' s' := if shouldRoundAway then -- we should round up. eufOut.incrSignificand - else + else eufOut @@ -611,4 +616,6 @@ theorem div_one_is_id_on_numerator_ex_smaller (a : PackedFloat 5 2) (h : ¬ a.un theorem div_self_is_one' (a : PackedFloat 5 2) (h : ¬a.isNaN ∧ ¬a.isInfinite ∧ ¬a.isZero) : (div a a .RTZ) = oneE5M2 := by + bv_normalize + simp at *; bv_decide diff --git a/Fp/Tactics/Simps.lean b/Fp/Tactics/Simps.lean index e3bc9c3..66d8467 100644 --- a/Fp/Tactics/Simps.lean +++ b/Fp/Tactics/Simps.lean @@ -24,6 +24,12 @@ open Lean Meta Simp in dsimproc [seval, simp, bv_normalize] reduceLog2 (Nat.log2 _) := Nat.reduceUnary ``Nat.log2 1 Nat.log2 +dsimproc [seval, simp, bv_normalize] reduceMax (Nat.max _ _) := + Nat.reduceBin ``Nat.max 2 Nat.max + +dsimproc [seval, simp, bv_normalize] reduceAdd (Nat.add _ _) := + Nat.reduceBin ``Nat.add 2 Nat.add + open Lean Meta Simp in simproc ↓ [bv_normalize] ite_eq_cond_proc (@ite _ _ _ _ _) := fun e => do let mkApp5 (.const ``ite [u]) α c hc t e := e | return .continue