feat(Falcon): decode FPR to ℝ, and prove all six of its error bounds - #514
Open
alik-eth wants to merge 37 commits into
Open
feat(Falcon): decode FPR to ℝ, and prove all six of its error bounds#514alik-eth wants to merge 37 commits into
alik-eth wants to merge 37 commits into
Conversation
The `FPR` denotation went through `Float.ofBits` and `Float.toRat0`, whose chain bottoms out on runtime primitives the kernel cannot evaluate, so nothing about `toReal` was provable. Decode the IEEE-754 binary64 fields from the word instead, in two layers: `FPR.decode` splits the sign, biased exponent and mantissa using `Nat` bit operations, and `FPR.Bits.toReal` denotes those fields, giving subnormals no implicit leading bit and sending non-finite patterns to zero. With the denotation reducible, the structural facts follow: `toReal_zero`, `toReal_one`, and `toReal_neg`, the last because the interpretation is odd in the sign bit in every branch. The per-operation rounding bounds remain open, and `HasRealSemantics FPR` still awaits them.
Reusable algebra on the IEEE-754 field decomposition, independent of any particular FPR operation: field bounds, a uniform reconstruction of the denotation as significand times a power of two that erases the subnormal/normal case split, the predicate carving out exactly representable reals, the spacing between adjacent representable values, and the bound relating that spacing to the magnitude of the value. The last of these is what any relative error bound ultimately rests on.
The per-operation error bounds cannot hold over all bit patterns: the emulation supports neither subnormal operands nor results outside the normal range, and a non-finite result denotes zero, so overflow breaks any relative bound. Name the domain on which the bounds are intended to hold: `FPR.IsNormal` for operands and `FPR.ExactInNormalRange` for the exact result. Also establish that magnitude order agrees with bit-pattern order (`abs_toReal_lt_iff_magKey_lt`), together with the low-63-bit mask bridge and the subtract-and-test-sign-bit comparison, which is what the addition algorithm's conditional swap rests on.
The four FPR operations share a small set of bit-manipulation idioms. This establishes each of them once, against the literal terms the kernels compute: - the sticky fold, which records in a single bit whether a right shift discarded anything, so alignment and renormalisation lose no information; - the leading-zero count behind normalisation, and that shifting by it and compensating the exponent preserves the denoted value; - the three-bit lookup shared by every operation's final assembly, shown to be round-to-nearest with ties to even, together with its half-ulp bound; - the assembly itself, whose mantissa carry propagates into the exponent field, and the round trip back through the field decomposition. The per-operation error bounds remain open.
Stated over all bit patterns, the four bounds are false. A finite result that overflows lands on the non-finite exponent, which denotes zero, so the error is the whole magnitude; and the emulation supports neither subnormal operands nor results that underflow out of the normal range. Square root needs no range condition on its result: the square root of a normal magnitude is hundreds of binades clear of both ends. Restrict the operands with `FPR.IsNormal` and, where the result can leave the representable range, restrict it with the magnitude bracket. Rename that bracket to `FPR.InNormalMagnitudeRange`, since it asserts nothing about exact representability. Every conclusion is unchanged. Witness that the hypotheses admit ordinary values, and that the result condition is not implied by the operand conditions: two normal operands can cancel into the subnormal range.
…ote the zero significand Remove the theory of `norm64`, which no code calls: the addition and scaling kernels inline its pattern, and lemmas already stated on that inlined form cover it. Restate the one fact unique to it about the term that occurs. Decide the comparison the addition kernel actually performs. It tests the sign of `za ||| ((za - 1) &&& x)`, not of `za`; off a tie these agree, and on a tie `za - 1` is all ones, so the test reads the sign of `x`. That settles the tie-break the module docstring recorded as open. Also drop a hypothesis its own siblings contradict, and denote `make_z` at a zero significand, the branch its flush-to-zero mask exists to serve.
The class asked every operation for a relative error bound on all inputs. No finite floating-point format can supply that: a result that leaves the representable range is not within any relative epsilon of the exact value, so the class had no instances and its only intended one stayed commented out. Carry the domain instead. A format now says which of its values the laws govern, and which exact results it can round; each operation bounds its error on that domain and, so that operations compose, states that its result stays in it. Negation keeps its unconditional law, being a sign flip. The compound bounds are unchanged; they take the domain conditions of the steps they chain.
🤖 PR Summary
Mathematical Formalization
Decode kernel. Proof Completion (
|
| component | units |
|---|---|
| fixed-point pipeline | 10 |
| 15 Chebyshev certificates | 3574 |
| Taylor truncation, deg. 19 | 80 |
| total | 3664 |
Fixed-point half: mulHi64 = high half of exact 128-bit product (toNat_mulHi64); mtwop63 = exactly ⌊2^63·toReal x⌋ (toNat_mtwop63) with no rounding. Twelve Horner floors contract because each step scales error by <0.694. UInt64 subtraction cannot wrap since mulHi64 a b ≤ b for any a. expm_p63_sub_trueArg_le bounds the total fixed-point deviation.
Certification layer. facctCoeffs is minimax, not Taylor: P−T has O(1) coefficients while being O(2^-51). Triangle-inequality bounds fail (coefficient-wise misses by ~38000×; midpoint+Lipschitz still misses by 730× at 256 subintervals). Subdivision restores locality: 15 Chebyshev certificates over dyadic breakpoints land at 3574. Dyadic avoids decimal log 2 numerals (which produce 337-digit text). certQ = Horner polynomial minus degree-18 Taylor; abs_le_of_chebCert + one concrete certificate cert0 for [0,1/16]. The certification layer costs ~25s build time; docs/agents/expm-certification.md documents the measurements and rejected approaches.
Protocols / Soundness
berExpReduce (new in SamplerZ.lean) reduces x ≥ 0 modulo log 2, returning (si, r) with floor rounding; berExp refactored to call it. Docstrings added for berExpReduce, berExp, samplerZLoop. Four tests added to runFalconProtocolTests: (1) k=1..199 berExpReduce yields r ∈ [0, log2Bound]; (2) berExpReduce zero = (0, zero); (3)/(4) xAt sign check across isigmaHi/isigmaLo boundary (σ₀).
Import alignment. FPRBridge.lean moves from Batteries.Data.Rat.Float to LatticeCrypto.Falcon.Concrete.FPR (with all modifier), ensuring consistent project-internal FPR definitions.
Infrastructure / CI
Full CI green; CI-lint clean; extern and interop isolation OK. Every new constant axiom-clean. No bv_decide/native_decide.
Documentation
docs/agents/expm-certification.md — working notes for the expm_p63 bound: obligation |expm_p63 x ccs / 2^63 − ccs·exp(−x)| ≤ 2^−51, the 4096-unit budget, the minimax fit, four ineffective bounding approaches, and the viable piecewise Chebyshev strategy. Records cost measurements (certificates over ℚ vs ℤ) and the three-part proof shape.
Architecture / Refactoring
File size note. FPRBridge.lean is now 6.8k lines (double the previous largest). A proposed six-way split is described (Decode / Kernel / Add / Mul / Loop / DivSqrt + re-export), deliberately deferred to a follow-up PR for cleaner diff attribution.
Original statements were false. Each of the five bounds was refuted by a kernel-checked ¬ (∀ …) against its original text (not included in this branch — they quote deleted statements — but available as evidence the new hypotheses are necessary).
Statistics
| Metric | Count |
|---|---|
| 📝 Files Changed | 7 |
| ✅ Lines Added | 8862 |
| ❌ Lines Removed | 94 |
Lean Declarations
✏️ Removed: 1 declaration(s)
Extern/Falcon/FPRBridge.lean (1)
theorem expm_p63_error (x ccs : FPR)
✏️ Added: 670 declaration(s)
Extern/Falcon/ExpmBridge.lean (138)
example : 0 ≤ toReal FPR.zero ∧ toReal FPR.zero ≤ 694 / 1000 ∧private def certN0 : ℕ → ℤprivate def certN1 : ℕ → ℤprivate def certN10 : ℕ → ℤprivate def certN11 : ℕ → ℤprivate def certN12 : ℕ → ℤprivate def certN13 : ℕ → ℤprivate def certN14 : ℕ → ℤprivate def certN2 : ℕ → ℤprivate def certN3 : ℕ → ℤprivate def certN4 : ℕ → ℤprivate def certN5 : ℕ → ℤprivate def certN6 : ℕ → ℤprivate def certN7 : ℕ → ℤprivate def certN8 : ℕ → ℤprivate def certN9 : ℕ → ℤprivate def half : FPRprivate def hornerMachine (z : UInt64) : ℕ → UInt64private noncomputable def certQ (x : ℝ) : ℝprivate noncomputable def hornerExact (ζ : ℝ) : ℕ → ℝprivate noncomputable def scaledArg (v : FPR) : ℝprivate noncomputable def taylorExpNeg (n : ℕ) (x : ℝ) : ℝprivate theorem abs_certQ_le (t : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 89 / 128) :private theorem abs_le_of_chebCert {p q s D B : ℝ} (hq : 0 < q) (hs : 0 < s) (hD : 0 < D)private theorem abs_taylorExpNeg_sub_exp_le (t : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 1) :private theorem and_63_of_lt {n : ℕ} (h : n < 64) : n &&& 63 = nprivate theorem and_63_or_65535 (n : ℕ) : (n ||| 65535) &&& 63 = 63private theorem cert0 (x : ℝ) (hx0 : (1 - 1 : ℝ) / 32 ≤ x)private theorem cert0_id (y : ℝ) :private theorem cert1 (x : ℝ) (hx0 : (9 - 1 : ℝ) / 128 ≤ x)private theorem cert10 (x : ℝ) (hx0 : (15 - 1 : ℝ) / 32 ≤ x)private theorem cert10_id (y : ℝ) :private theorem cert11 (x : ℝ) (hx0 : (17 - 1 : ℝ) / 32 ≤ x)private theorem cert11_id (y : ℝ) :private theorem cert12 (x : ℝ) (hx0 : (19 - 1 : ℝ) / 32 ≤ x)private theorem cert12_id (y : ℝ) :private theorem cert13 (x : ℝ) (hx0 : (21 - 1 : ℝ) / 32 ≤ x)private theorem cert13_id (y : ℝ) :private theorem cert14 (x : ℝ) (hx0 : (177 - 1 : ℝ) / 256 ≤ x)private theorem cert14_id (y : ℝ) :private theorem cert1_id (y : ℝ) :private theorem cert2 (x : ℝ) (hx0 : (11 - 1 : ℝ) / 128 ≤ x)private theorem cert2_id (y : ℝ) :private theorem cert3 (x : ℝ) (hx0 : (13 - 1 : ℝ) / 128 ≤ x)private theorem cert3_id (y : ℝ) :private theorem cert4 (x : ℝ) (hx0 : (15 - 1 : ℝ) / 128 ≤ x)private theorem cert4_id (y : ℝ) :private theorem cert5 (x : ℝ) (hx0 : (5 - 1 : ℝ) / 32 ≤ x)private theorem cert5_id (y : ℝ) :private theorem cert6 (x : ℝ) (hx0 : (7 - 1 : ℝ) / 32 ≤ x)private theorem cert6_id (y : ℝ) :private theorem cert7 (x : ℝ) (hx0 : (9 - 1 : ℝ) / 32 ≤ x)private theorem cert7_id (y : ℝ) :private theorem cert8 (x : ℝ) (hx0 : (11 - 1 : ℝ) / 32 ≤ x)private theorem cert8_id (y : ℝ) :private theorem cert9 (x : ℝ) (hx0 : (13 - 1 : ℝ) / 32 ≤ x)private theorem cert9_id (y : ℝ) :private theorem certQ_expand (x : ℝ) :private theorem certQ_shift (s u v : ℝ) (hv : s * u = v) :private theorem chebEval0 (y : ℝ) : (Chebyshev.T ℝ (0 : ℤ)).eval y = 1private theorem chebEval1 (y : ℝ) : (Chebyshev.T ℝ (1 : ℤ)).eval y = yprivate theorem chebEval10 (y : ℝ) : (Chebyshev.T ℝ (10 : ℤ)).eval y =private theorem chebEval11 (y : ℝ) : (Chebyshev.T ℝ (11 : ℤ)).eval y =private theorem chebEval12 (y : ℝ) : (Chebyshev.T ℝ (12 : ℤ)).eval y =private theorem chebEval13 (y : ℝ) : (Chebyshev.T ℝ (13 : ℤ)).eval y =private theorem chebEval14 (y : ℝ) : (Chebyshev.T ℝ (14 : ℤ)).eval y =private theorem chebEval15 (y : ℝ) : (Chebyshev.T ℝ (15 : ℤ)).eval y =private theorem chebEval16 (y : ℝ) : (Chebyshev.T ℝ (16 : ℤ)).eval y =private theorem chebEval17 (y : ℝ) : (Chebyshev.T ℝ (17 : ℤ)).eval y =private theorem chebEval18 (y : ℝ) : (Chebyshev.T ℝ (18 : ℤ)).eval y =private theorem chebEval2 (y : ℝ) : (Chebyshev.T ℝ (2 : ℤ)).eval y =private theorem chebEval3 (y : ℝ) : (Chebyshev.T ℝ (3 : ℤ)).eval y =private theorem chebEval4 (y : ℝ) : (Chebyshev.T ℝ (4 : ℤ)).eval y =private theorem chebEval5 (y : ℝ) : (Chebyshev.T ℝ (5 : ℤ)).eval y =private theorem chebEval6 (y : ℝ) : (Chebyshev.T ℝ (6 : ℤ)).eval y =private theorem chebEval7 (y : ℝ) : (Chebyshev.T ℝ (7 : ℤ)).eval y =private theorem chebEval8 (y : ℝ) : (Chebyshev.T ℝ (8 : ℤ)).eval y =private theorem chebEval9 (y : ℝ) : (Chebyshev.T ℝ (9 : ℤ)).eval y =private theorem decode_half : FPR.decode half = ⟨false, 1022, 0⟩private theorem expm_p63_eq (x ccs : FPR) :private theorem expm_p63_sub_exact_le (x ccs : FPR)private theorem expm_p63_sub_trueArg_le (x ccs : FPR)private theorem exponent_le_of_toReal_lt_one (x : FPR) (hn : FPR.IsNormal x)private theorem facctCoeffs_le {i : ℕ} (hi : i ≤ 12) : (facctCoeffs[i]!).toNat ≤ 2 ^ 63private theorem facctCoeffs_mono {i : ℕ} (hi : i < 12) :private theorem facctVal0 : (facctCoeffs[0]!).toNat = 19127174051private theorem facctVal1 : (facctCoeffs[1]!).toNat = 233346759686private theorem facctVal10 : (facctCoeffs[10]!).toNat = 4611686018427565056private theorem facctVal11 : (facctCoeffs[11]!).toNat = 9223372036854728704private theorem facctVal12 : (facctCoeffs[12]!).toNat = 9223372036854775808private theorem facctVal2 : (facctCoeffs[2]!).toNat = 2542029181962private theorem facctVal3 : (facctCoeffs[3]!).toNat = 25415798087749private theorem facctVal4 : (facctCoeffs[4]!).toNat = 228754078003076private theorem facctVal5 : (facctCoeffs[5]!).toNat = 1830034511206115private theorem facctVal6 : (facctCoeffs[6]!).toNat = 12810238987800554private theorem facctVal7 : (facctCoeffs[7]!).toNat = 76861433589428176private theorem facctVal8 : (facctCoeffs[8]!).toNat = 384307168197152512private theorem facctVal9 : (facctCoeffs[9]!).toNat = 1537228672812056320private theorem floor_two_pow_63_mul_toReal (x : FPR) (hn : FPR.IsNormal x)private theorem hornerExact_abs_le (t : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 694 / 1000) (i : ℕ)private theorem hornerExact_lipschitz (s t : ℝ) (hs0 : 0 ≤ s) (ht0 : 0 ≤ t)private theorem hornerExact_mem (t : ℝ) (ht0 : 0 ≤ t) (ht1 : t ≤ 694 / 1000) :private theorem hornerMachine_error {z : UInt64} {ζ : ℝ}private theorem hornerMachine_le_coeff (x : FPR) (h0 : 0 ≤ toReal x)private theorem hornerMachine_sub_hornerExact_le (x : FPR)private theorem horner_step {z y c : UInt64} {Y ζ : ℝ}private theorem mtwop63_neg (x : FPR) : mtwop63 (FPR.neg x) = mtwop63 xprivate theorem mulHi64_bracket (a b : UInt64) :private theorem mulHi_limbs (aHi aLo bHi bLo : ℕ) :private theorem or_two_pow_mod_two_pow_succ (A k : ℕ) :private theorem scaledArg_bracket (v : FPR)private theorem scaledArg_eq (v : FPR) (h0 : 0 ≤ toReal v)private theorem scaledArg_le_694 (x : FPR) (h0 : 0 ≤ toReal x)private theorem scaledArg_le_one (v : FPR) : scaledArg v ≤ 1private theorem scaledArg_le_toReal (x : FPR) (h0 : 0 ≤ toReal x)private theorem scaledArg_nonneg (v : FPR) : 0 ≤ scaledArg vprivate theorem toNat_add3_of_lt {x y z : UInt64}private theorem toNat_add4_of_lt {x y z w : UInt64}private theorem toNat_and_low32 (v : UInt64) : (v &&& 0xFFFFFFFF).toNat = v.toNat % 2 ^ 32private theorem toNat_e_of (x : FPR) (hexle : (FPR.decode x).exponent ≤ 1022) :private theorem toNat_finalColumn {X Y Z u v w : UInt64} {A B : ℕ}private theorem toNat_m_of (x : FPR) :private theorem toNat_mtwop63 (x : FPR)private theorem toNat_mtwop63_aux (x : FPR) (hexle : (FPR.decode x).exponent ≤ 1022) :private theorem toNat_mtwop63_of_exponent_eq_2047 (x : FPR)private theorem toNat_mtwop63_of_exponent_eq_zero (x : FPR)private theorem toNat_mulHi64 (a b : UInt64) :private theorem toNat_mulHi64_le (a b : UInt64) : (mulHi64 a b).toNat ≤ b.toNatprivate theorem toNat_shiftRight_32_uint64 (v : UInt64) : (v >>> 32).toNat = v.toNat / 2 ^ 32private theorem toNat_shiftRight_ue {m : UInt64} {e : UInt32} {k : ℕ} (he : e.toNat = k)private theorem toNat_ue_of (e : UInt32) (he : e.toNat < 2 ^ 16) :private theorem toReal_eq_significand_of_nonneg (x : FPR) (hn : FPR.IsNormal x)private theorem toReal_half : toReal half = 0.5private theorem toReal_lt_of_exponent_eq_zero (x : FPR)private theorem two_zpow_neg_63 : (2 : ℝ) ^ (-63 : ℤ) = 1 / 2 ^ (63 : ℕ)theorem expm_p63_error (x ccs : FPR)theorem expm_p63_error_abs (x ccs : FPR)theorem expm_p63_neg (x ccs : FPR) : FPR.expm_p63 (FPR.neg x) ccs = FPR.expm_p63 x ccs
Extern/Falcon/FPRBridge.lean (531)
def FPR.Bits.IsNormal (b : FPR.Bits) : Propdef FPR.Bits.magKey (b : FPR.Bits) : ℕdef FPR.Bits.significand (b : FPR.Bits) : ℕdef FPR.Bits.workExp (b : FPR.Bits) : ℕdef FPR.InNormalMagnitudeRange (r : ℝ) : Propdef FPR.IsNormal (x : FPR) : Propdef FPR.decode (x : FPR) : FPR.Bits wheredef FPR.maxFiniteReal : ℝdef FPR.minNormalReal : ℝdef FPR.ulpOfExponent (e : ℕ) : ℝdef IsFPRRepresentable (r : ℝ) : Propdef roundQuarterTiesEven (n : ℕ) : ℕ
…and 520 more not listed.
✏️ Affected: 12 declaration(s) (line number changed)
theorem add_result_bounds (a b : F) (ha : self.Valid a) (hb : self.Valid b)inExtern/Falcon/ApproxArith.leanmoved from L91 to L128theorem butterfly_add_error (a b w : F)inExtern/Falcon/ApproxArith.leanmoved from L166 to L217theorem butterfly_sub_error (a b w : F)inExtern/Falcon/ApproxArith.leanmoved from L190 to L245theorem compound_add_mul_error (a b c d : F)inExtern/Falcon/ApproxArith.leanmoved from L108 to L148theorem horner_step_error (a x b : F)inExtern/Falcon/ApproxArith.leanmoved from L139 to L186theorem mul_result_bounds (a b : F) (ha : self.Valid a) (hb : self.Valid b)inExtern/Falcon/ApproxArith.leanmoved from L98 to L137def toReal (x : FPR) : ℝinExtern/Falcon/FPRBridge.leanmoved from L63 to L103def verifyPrimitives (p : Falcon.Params) (hn : p.n = 2 ^ p.logn) : Falcon.Primitives p whereinExtern/Falcon/FPRBridge.leanmoved from L70 to L3815theorem add_error (a b : FPR) (ha : FPR.IsNormal a) (hb : FPR.IsNormal b)inExtern/Falcon/FPRBridge.leanmoved from L84 to L3835theorem div_error (a b : FPR) (hb : toReal b ≠ 0) (ha : FPR.IsNormal a) (hb' : FPR.IsNormal b)inExtern/Falcon/FPRBridge.leanmoved from L96 to L5767theorem mul_error (a b : FPR) (ha : FPR.IsNormal a) (hb : FPR.IsNormal b)inExtern/Falcon/FPRBridge.leanmoved from L90 to L4899theorem sqrt_error (a : FPR) (ha' : FPR.IsNormal a) (ha : 0 ≤ toReal a) :inExtern/Falcon/FPRBridge.leanmoved from L102 to L6555
sorry Tracking
✅ Removed: 5 `sorry`(s)
Extern/Falcon/FPRBridge.lean (5)
theorem add_error (a b : FPR) :(L87)theorem div_error (a b : FPR) (hb : toReal b ≠ 0) :(L99)theorem expm_p63_error (x ccs : FPR)(L115)theorem mul_error (a b : FPR) :(L93)theorem sqrt_error (a : FPR) (ha : 0 ≤ toReal a) :(L105)
Coverage Notes
- AI file summarization partially analyzed 2 file(s) because their individual diffs exceeded the per-file size budget. Statistics and Lean signal tracking still cover the full PR.
Partially Analyzed Files
Extern/Falcon/ExpmBridge.lean(+1896/-0)Extern/Falcon/FPRBridge.lean(+6657/-39)
📄 **Per-File Summaries**
- Extern.lean: In
Extern.lean, a newpublic import Extern.Falcon.ExpmBridgeline has been added. This makes theExpmBridgemodule, which likely provides the exponential function bridge for the Falcon project, publicly accessible to downstream modules. - Extern/Falcon/ApproxArith.lean: The
HasRealSemanticsclass inExtern/Falcon/ApproxArith.leanis extended with two new predicates —Valid : F → PropandInRange : ℝ → Prop— and six new closure fields (add_valid,mul_valid,div_valid,sqrt_valid,neg_valid,sub_valid) that guarantee operands stay valid after each operation when the exact result is in range. All existing error fields (add_error,mul_error,div_error,sqrt_error,sub_error) now require the caller to supplyValidandInRangehypotheses for the operands and exact result. The derived theoremsadd_result_bounds,mul_result_bounds,compound_add_mul_error,horner_step_error,butterfly_add_error, andbutterfly_sub_errorare updated to accept these same hypotheses, and their proofs use the new*_validfields to derive validity of intermediate results. The commented‑outFPRinstance is rewritten: it records thatinterpis now a pure bit‑field decoding (sointerp_zero,interp_one,neg_exactare fully proved), that the error fields are available fromFPRBridge, and that the remaining open obligations are the_validclosure fields andsub_error(which should reduce compositionally). The instance remains commented out pending those proofs. - Extern/Falcon/ExpmBridge.lean: This new file introduces
Extern/Falcon/ExpmBridge.lean, the fixed-point verification layer for Falcon'sexpm_p63function. It provides a complete formal proof that theUInt64Horner machinehornerMachineand the finalmulHi64scaling inexpm_p63track the exact real Horner recurrencehornerExactto within an error bound of10units of2^(-63)when the argumenttoReal xis in[0, 694/1000]and the scaletoReal ccslies in[0,1). Key steps include: the theoremtoNat_mulHi64establishing thatmulHi64computes the high half of the 128-bit product; the theoremtoNat_mtwop63proving thatmtwop63yields the exact floor⌊2^63 * toReal x⌋forxin[0,1); and the cumulative error theoremexpm_p63_sub_trueArg_lewhich bounds the total deviation. The file also provides the structural skeleton for a Chebyshev-based numerical certification of the coefficient polynomial againstReal.exp, specifically the definitioncertQ(the difference between the Horner polynomial and a degree‑18 Taylor expansion) and an abstract certificate lemmaabs_le_of_chebCert, together with one concrete certificate (cert0) for the subinterval[0, 1/16]that bounds|certQ| ≤ 3045. Nosorryoradmitstatements appear in the diff. - Extern/Falcon/FPRBridge.lean: The import
Batteries.Data.Rat.Floatis removed, and the importLatticeCrypto.Falcon.Concrete.FPRis added (with theallmodifier, replacing the earlierpublic importof that same module). This aligns the file's dependency on floating-point representations with the project's ownConcrete.FPRmodule rather than the externalBatterieslibrary, ensuring that all definitions and theorems related to FPR (floating-point rationals) are consistently drawn from the project-internal source. - Extern/Falcon/SamplerZ.lean: The diff adds the new function
berExpReduce(which reducesx ≥ 0modulolog 2returning(si, r)with floor rounding), refactorsberExpto call it instead of inlining the reduction, and adds docstrings forberExpReduce,berExp, andsamplerZLoopthat explain their semantics, thex ≥ 0precondition, and the correctness reasoning behind the saturation logic and the rejection probability. - LatticeCryptoTest/Falcon/Main.lean: Added four new test checks to
runFalconProtocolTests. The first verifies that forkfrom 1 to 199,berExpReduce (scaled (Int64.ofNat k) (-5))returns anrthat is non-negative and ≤log2Bound(the representable upper bound oflog 2), reporting the count of failures and the first failingk. The second checksberExpReduce (F := FPR) zeroreturns(0, zero). The third and fourth tests examine the sign ofxcomputed by the localxAtfunction (which models the exponent argument insamplerZ), using two adjacent representable values of1/σ(isigmaHiandisigmaLo) straddling the semantic1/σ₀boundary: the third asserts zero negative values amongz0in[0:26]forisigmaHi(σ ≤ σ₀), and the fourth asserts at least one negative value forisigmaLo(σ > σ₀ by one ulp). Nosorryoradmitwere added. - docs/agents/expm-certification.md: This documentation file adds working notes for resolving the final
sorryinExtern/Falcon/FPRBridge.lean, which involves certifying the error bound forexpm_p63onx ∈ [0, log 2)andccs ∈ [0, 1). It specifies the obligation as|expm_p63 x ccs / 2^63 - ccs * exp(-x)| ≤ 2^(-51), defines the per-unit budget of 4096 units, and explains the tightness due toFPR.facctCoeffsbeing a minimax fit (not a Taylor truncation) with uniform error of 3562 units. The document surveys four ineffective bounding approaches (Taylor coefficient-wise, Chebyshev whole-interval, midpoint+Lipschitz, and coefficient-wise Lipschitz), then presents viable piecewise strategies (Taylor form and Chebyshev) with their bounds and margins, concluding that Chebyshev at 32 pieces is the preferred basis for certificate cost reasons. It also discusses certificate generation cost overℚvs.ℤ, reporting that clearing denominators yields 2.54 s per certificate (compared to 11.5 s overℚ), estimates total certification build times (~88 s for Chebyshev/32, ~332 s for Taylor/128), and outlines the three-part proof shape: fixed-point layer, certification layer, and assembly.
Last updated: 2026-08-14 09:29 UTC.
…dition bound Pin the addition kernel to a named record of its intermediates, proved equal to the kernel by reduction, so later steps can name what each stage computes. On that footing: the significand handed to the final assembly lies in the normalized range, the conditional swap preserves the operand pair and orders it by magnitude, and the extracted fields denote the operands' significands, exponents and signs. The addition bound itself remains open.
The bound now follows from the pipeline record: the conditional swap preserves the sum and orders the magnitudes, the aligned significand is the sticky shift of the smaller operand and so within one unit in the last place of it, renormalisation preserves the value, and the final assembly rounds to nearest. The alignment allowance is not the naive half of the budget. Matching signs cost a quarter of an ulp; differing signs with a gap of at most three are exact, the significands being scaled by eight; and beyond that cancellation is bounded away. The assembly's own range condition also needed widening by one binade, since a sum of two largest-magnitude normals still lies inside the representable range. The range hypothesis earns itself three times over: it bounds the exponent from below, rules out the carry at the top binade, and its zero case is exact cancellation.
alik-eth
changed the title
Subtraction is addition against a negated operand, and negation flips the sign bit, leaving the exponent field and so normality alone. The bound therefore transfers from addition with no further rounding analysis.
alik-eth
changed the title
The bound constrained x but placed no condition on ccs. The routine's output is a UInt64 read at scale 2^63, so the left side is always below 2, while toReal ccs * exp (-toReal x) is unbounded; and at ccs = 1 the fixed-point conversion wraps to zero, collapsing the product for every x in range. Require 0 <= toReal ccs and toReal ccs < 1, and add a witness showing the four side conditions are jointly satisfiable. The 2^-51 conclusion is unchanged: checked against 300-bit reference arithmetic over 460k operand pairs, the worst error is 3562 units of 2^-63 against a budget of 4096.
Models FPR.mul as a MulPipeline record pinned to the kernel term by rfl, then proves the 2^-52 relative bound on normal operands whose exact product is representable. The carry chain is proved twice over: the 25-bit limb identity as pure arithmetic on N, away from machine words, and separately a no-overflow bound for every intermediate (the widest, z2, reaches 2^29). Together they give zu = xu*yu / 2^50 with the two residues carrying exactly the discarded low bits. The discarded bits are folded twice, at bit 50 and again at bit 1 when the product carries into bit 55. Sticky shifts compose, so the pair is a single sticky shift and truncation costs one ulp rather than two: 2^-53 for rounding plus 2^-54 for truncation, a quarter of the budget spare. Composing the two brackets separately would have spent 7/8 of it. The flush-to-zero guard is inactive on normal operands, and the wrapping UInt32 exponent lands on the two's-complement encoding of the intended value. As for addition, the result-range hypothesis is what pins the exponent into the window the rounding step needs. Sorry count in the module drops from four to three.
alik-eth
changed the title
FPR.div and FPR.sqrt are the two kernels that are not straight-line: each runs a fixed-length `for _ in [0:n]` loop over a mutable state, so neither could be pinned to a model by rfl the way add and mul were. Neither loop body reads the index, so the fold is iteration of a single step function. Core's forIn_eq_forIn_range' reduces the Range fold to a List.range' fold; induction on that gives forIn = f^[n], and from there an invariant of the loop state is ordinary induction on the step count. Adds DivPipeline and SqrtPipeline in the style of addPipeline and mulPipeline, each naming the loop result as a field pinned by rfl, plus divPipeline_loop_induction and sqrtPipeline_loop_induction as the entry points. divPipeline_quotient_lt carries the first invariant through the 55 iterations, bounding the quotient handed to the renormalisation step. Shared prerequisite for div_error and sqrt_error; no sorry is closed by this commit.
Carries the restoring-division invariant through all 55 iterations: yu * q + r = 2^55 * xu with the running remainder held below 2 * yu. That pins the quotient handed to the sticky and renormalisation steps, and is the mathematical core of div_error. Route: divStep_mask reads the loop's comparison as a plain comparison (the bounds rule out the wraparound the sign-bit test would otherwise pick up), divStep_toNat reads one step as conditional subtraction plus a doubling, and divLoop_invariant inducts over the iterate. The step count has to be bounded for the invariant to hold: past 62 iterations the quotient's doubling would overflow UInt64, so the lemma is stated for n <= 60 and instantiated at 55. Adds xu / yu fields to DivPipeline to mirror MulPipeline, and factors the extended-significand computation into significand_pack_toNat, now shared by the mul, div and sqrt pipelines. No sorry closed; div_error still needs its rounding and assembly half.
Closes the gap between the division loop's output and the exact quotient: the word handed to the renormalisation step satisfies |q0 * yu - 2^55 * xu| < yu, one unit in the last place either way. The sticky bit is what buys the sharpness. The loop leaves an even quotient and a remainder below 2 * yu, which alone only pins the quotient to two units; folding in `r ||| (0 - r)` — nonzero exactly when the division was inexact — recovers the odd values and halves the bracket. Measured against exact rational arithmetic over 150k in-domain operand pairs, div's end-to-end error peaks at 0.498 of the 2^-52 budget, consistent with one ulp before rounding plus 2^-53 of rounding. Also proves the closing sticky test itself: r ||| (0 - r) has its top bit set exactly when r is nonzero, since one of a word and its two's complement negation always does.
Carries the one-ulp quotient bracket through the renormalisation step: the significand handed to make brackets the exact 2^55 * xu / yu to one unit in its own last place, and lands in [2^54, 2^55) as make's rounding analysis requires. The sharpness needs a lemma the generic sticky theory does not give. stickyShift_mul_lt / lt_stickyShift_mul_add allow two units of movement, which composed with the quotient's own unit would overrun the budget; stickyShift_one_bracket pins a one-bit sticky shift to a single unit, which is what it actually costs. Also clears the flush-to-zero guard: on a normal numerator the exponent field is at least 1, so tbmask is zero and neither the exponent nor the significand is rewritten. No sorry closed; div_error still needs its exponent window and assembly.
Closes div_error: the 2^-52 relative bound on normal operands whose exact quotient is representable. Sorry count in the module drops from three to two. The pieces: the restoring-division loop's invariant and the one-ulp quotient bracket (already landed), the exponent chain through the wrapping UInt32 arithmetic, the magnitude and sign of the exact quotient, the exponent window forced by the result-range hypothesis, and the no-carry case at the top of that window. mul_error_combine assembles it unchanged, instantiated with the renormalised quotient, 2^es as the scale between its units and the exact value's, and 2^(ex-ey-55) as the common exponent. Writing that lemma abstractly rather than inlining it into mul_error is what made division reuse it directly.
alik-eth
changed the title
The 2^-51 in that bound is very nearly saturated, and by the approximation rather than the arithmetic: facctCoeffs is a minimax fit whose uniform error against exp(-x) over [0, log 2) is already about 2^-51.2, roughly 87% of the bound. The fixed-point conversion and the twelve Horner steps get the remaining eighth. Perturbing the table or widening either operand range is liable to make the statement false, not just harder to prove. That also constrains how it can be proved. Comparing the table to the Taylor coefficients term by term and adding a Lagrange remainder lands around 2^-35.8, four orders of magnitude short, because a minimax fit earns its accuracy from cancellation across the interval and a per-coefficient triangle inequality discards exactly that. The bound needs a rigorous enclosure of |P x - exp (-x)| over the interval instead, tight to within the eighth of budget left over.
Carries the square-root invariant through all 54 iterations: writing q for the partial root, s for its double, r for the bit weight and xp for the scaled remainder, the loop maintains q^2 * 2^k + 2^53 * xp = 2^(54+k) * xu and xp < 4q + 2^(55-k) alongside s = 2q and r = 2^(53-k). The invariant and the bit weight are mutually dependent: the quadratic term telescopes precisely because r * 2^(k+1) = 2^54, so both must be carried together. The comparison mask is generalised out of the division proof at the same time. sub_mask_of_lt reads (a - c) >>> 63 - 1 as a plain comparison for any operands below 2^63, which is what both loops need; divStep_mask now delegates to it, dropping 25 lines of duplicated wraparound reasoning. No sorry closed; sqrt_error still needs its integer-root to Real.sqrt bridge, its exponent analysis and its assembly.
The bit weight reaches zero on the final iteration, so the r = 2^(53-k) conjunct of sqrtLoop_invariant stops at k = 53. The remaining conjuncts survive one more step, and that step is what turns the invariant into an exact integer square root: q^2 <= 2^54 * xu < (q+1)^2 at k = 54. nat_sqrt_bracket then carries that to the reals, bracketing Real.sqrt (2^54 * xu) between q and q + 1. That bridge was the piece with no analogue elsewhere in the file, and it turns out to be short once the loop is known to compute a floor rather than an approximation. The closing step is proved by ring rather than omega. Stating it over 2^107 and 2^108 directly exhausts the recursion limit even at eight thousand; naming the remainder difference and rewriting through 2^108 * xu = 2 * (2^107 * xu) makes it a polynomial identity instead, and needs no limit raised. No sorry closed; sqrt_error still needs its exponent analysis and assembly.
Closes sqrt_error. The module's only remaining sorry is expm_p63_error. Both exponent parities collapse to a single formula: the operand is exactly xu * 2^(2e - 52), with xu the possibly-doubled significand and e the halved exponent, because the odd case's doubling of xu compensates the floor in e = arsh32(e_, 1) precisely. That removes what would otherwise be a duplicated case split running through the whole proof. The exponent needs no result-range hypothesis, which is what the statement always claimed: a normal operand gives e in [-511, 511], so e' = e - 54 lands in [-565, 457], well inside make_z's window, and the carry case at the top of that window cannot arise. Assembly is therefore strictly simpler than div_error's. mul_error_combine applies once more, here with K = 1 and the exact root 2 * sqrt (2^54 * xu) as the target value; the bracket is one unit because the loop computes a floor and the sticky bit recovers the odd case.
alik-eth
changed the title
Records what a proof of the last sorry has to beat, with measurements rather than estimates. The bound is 87% consumed by the minimax polynomial alone, leaving about 534 units of 2^-63 for everything else. Every bound that applies a triangle inequality to coefficients dies to the same cancellation: the fit earns its accuracy across the interval, while P - T_n has O(1) coefficients and is O(2^-51). Taylor coefficient-wise misses by ~38000x, and midpoint-plus-Lipschitz still misses by 730x at 256 subintervals, because bounding sup|Q'| coefficient-wise has the identical pathology. Subdivision works. Chebyshev on 32 pieces lands at 3585 units, matching Taylor form on 128 - a 4x reduction in certificates for the price of needing |T_j| <= 1. Taylor degree turns out irrelevant to the bound; 18 suffices. The binding constraint is build time, not mathematics: certificates cost about 15s each, and reducing the degree from 26 to 18 made it slower, not faster, since the cost tracks the rational arithmetic. Two untested levers are recorded in priority order, along with the observation that the fixed-point half can be built and landed on its own.
Working the certificates over Z instead of Q takes the marginal cost from about 11.5s to 2.54s each, measured by putting 1, 4 and 8 certificates in one file and taking the slope (9.44s / 16.90s / 27.04s, on a 6.8s import). The digit counts were denominator-dominated, so one common scaling removes most of the arithmetic. That puts Chebyshev on 32 pieces at about 88 seconds, which is unremarkable for this repo, and settles the basis question: Taylor form would need 128 pieces and about 332 seconds, so the |T_j| <= 1 lemma Chebyshev requires is worth paying for. Build time is no longer the binding constraint, and the second lever - proving the binomial shift once and instantiating with norm_num - is not needed. Also records that reducing the Taylor degree made things slower, not faster, since the cost tracks coefficient arithmetic and not degree.
Adds Extern/Falcon/ExpmBridge.lean, a new file rather than more weight on FPRBridge, which is already an outlier at 6.8k lines. mulHi_limbs is proved: the high half of a 64-bit product equals the top partial product plus the carries out of the two cross terms and the middle column. That is the whole mathematical content of mulHi64, stated as arithmetic on N with every quotient and remainder named so the identity is genuinely polynomial - ring cannot see through div and mod atoms, which is what defeated the first attempt. toNat_mulHi64, the UInt64 wrapper, is left as a documented sorry: what remains is pushing toNat through the unfolded expression tree, one lemma per node with a no-overflow side condition, all of which are already available in the proof context. WIP branch; nothing here is on the PR path.
Adds every fact the wrapper needs: the four limb bounds, the four partial products, the high and low halves of each, the middle column's no-overflow bound, and bounds on all four summands of the final column. The remaining step threads toNat through four nested UInt64 additions and closes with mulHi_limbs. Three syntax traps cost most of the time and are worth recording: >>> binds tighter than *, so (a >>> 32) * x >>> 32 parses as (a >>> 32) * (x >>> 32); mask literals need an explicit : UInt64 ascription or they elaborate as N and HAnd fails to synthesize; and the pretty-printer drops exactly the parentheses that disambiguate both, so the displayed goal is misleading.
mulHi64 is now proved to be the high half of the exact product, sorry-free and axiom-clean. Corrects the previous commit's diagnosis, which was wrong. It claimed the outermost no-overflow condition failed because omega lacked b11. Supplying b11 could not have helped: the four summand bounds are jointly too weak, since (2^32-1)^2 + 3*(2^32-1) exceeds 2^64 by 2^32 - 2. No bound on the summands can discharge that condition, because the sum genuinely can be made to look too large by them alone. The fact has to come from the column's value rather than its parts: the sum equals a*b/2^64, which is below 2^64 because a and b are. So the assembly establishes the value first via mulHi_limbs and reads the bound off it. b11, b10, b01 and bmid are consequently dead and removed. Adds toNat_add3_of_lt / toNat_add4_of_lt for register-fitting sums, and states the assembly over short variables so the masked and shifted terms are matched by unification and never retyped - which sidesteps the operator-precedence and literal-ascription traps entirely rather than navigating them. Also wraps a pre-existing 101-character line.
Proves mtwop63 x = floor(2^63 * toReal x) for a normal operand in [0, 1): the conversion introduces no rounding of its own, so all of expm_p63's fixed-point error lives in the Horner steps. The three machine facts underneath: the shift-and-or leaves m = 2^10 * significand, since only the exponent's low bit survives at position 62 and the implicit-bit OR overwrites it; e = 1022 - ex without wraparound, because toReal x < 1 forces ex <= 1022; and the (63 - e) >>> 16 idiom saturates the shift at 63, proved as an exact Nat.testBit fact rather than a bound. The saturating branch needs no exclusion: when it bites the operand is below 2^-64 and both sides are zero. Adds Mathlib.Algebra.Order.Floor.Semifield to the imports for Nat.floor_div_natCast, which is what turns the exact real quotient into Nat division and is not in FPRBridge's closure.
The fixed-point half was estimated at ~15 units of 2^-63. Measured, it is 1.55: mtwop63 is exact, so the only truncation is the twelve mulHi64 floors plus the final one, and each step's floor costs at most 1 against a previous error scaled by zeta < 0.694. Geometric series bounds the whole pipeline below about 4. That figure was the one the entire margin analysis rested on, and it was a guess. It is now measured, and the certification layer has about 530 units of slack rather than 534 minus a 15-unit reserve.
Completes the fixed-point half. The pipeline tracks the exact real Horner recurrence at the same quantised argument to within 4 at every iterate, and to within 5 including the final scaling multiply - against a budget of 4096 units, so the fixed-point work is essentially free and the whole margin belongs to the polynomial certification. The UInt64 subtraction in the loop body cannot wrap, and the reason is simpler than expected: mulHi64 a b <= b for every a, since a < 2^64 makes a*b/2^64 < b. So mulHi64 z y <= y <= facctCoeffs[i], needing only that the coefficient table is nondecreasing - no error bound, no coefficient magnitudes, no dependence on zeta. Carried as a conjunct of the same induction, since that is what chains, but the bound is free. The loop reads its index, so forIn_range_eq_iterate does not apply; the twelve iterations are unrolled by a single simp only through forIn_eq_forIn_range' and List.range'_succ. rfl alone exhausts the recursion limit, but that route does not need it raised. Adds Mathlib.Analysis.Complex.ExponentialBounds for Real.log_two_lt_d9. Only scaledArg_le_694 uses it; every lemma below takes the numeric bound directly, so the import can be dropped by restating the hypothesis as toReal x <= 694/1000 if that is preferred.
Completes the fixed-point half. The pipeline tracks
toReal ccs * hornerExact (toReal x) 12 to within 10, against a budget of
4096 units, so all of the margin belongs to the polynomial
certification.
The constant is 10 rather than the 8 estimated, because the sketch asked
for a Lipschitz constant the induction cannot deliver. The step is
d_i <= t*d_{i-1} + 2^63*|s-t|, whose fixed point needs at least
2^63/(1 - 0.694) = 3.27 * 2^63; 2^63 was the true constant but not a
provable one. Using 2^65 charges 4 for the argument gap instead of 1.
Sharpening it needs per-index magnitude bounds on hornerExact, which
duplicates certification-layer machinery for margin that is not needed.
The bridge is usually exact rather than merely small: 2^63 * toReal v is
an integer whenever the exponent field is at least 1012, i.e. whenever
toReal v >= 2^-11, so the floor does nothing and scaledArg v = toReal v.
Only operands below 2^-11 pay anything, and there the argument is
negligible. Realized end-to-end error on samples is under 1 unit against
the proved 10.
Builds the per-subinterval bound - if Q composed with the affine map to [-1,1] equals a Chebyshev sum, then |Q| is bounded by the sum of the coefficients' magnitudes - on Mathlib's abs_eval_T_real_le_one, plus one worked certificate over [0, 1/16]. The breakpoints must be dyadic. Anchoring the split at a decimal upper bound for log 2 gives coefficients with 337-digit numerators, and the longLine linter fires at 100 characters on a numeral that cannot be wrapped, so that design is unimplementable rather than merely awkward. With endpoints k/2^m the clearing scale is 18! * 2^(18(m+1)) and the numerals drop to about 63 digits. Dyadic breakpoints also need fewer pieces, since narrower intervals near the ends buy more than uniformity: 15 pieces reach 3573 units against a 4096 budget, where 32 uniform pieces reached 3585. Measured marginal cost is 1.26s per certificate on 7.8s of machinery, so the whole layer is about 25 seconds rather than the 88 estimated. certQ_expand and certQ_shift clear 18! and s^18 once each, so the per-certificate ring is integral and division-free. Sampled sup on the first piece is 2893 against the certified 3045 - tight, not vacuous.
Closes the last of the six bounds. The certification layer lands: 15 Chebyshev certificates over dyadic breakpoints, plus the fixed-point pipeline and the Taylor remainder. fixed-point pipeline 10 15 Chebyshev certificates 3574 Taylor truncation, n = 19 80 total 3664 <= 4096 Margin 10.5%. The Taylor term is 80 rather than negligible only because |-t|^19 is bounded by 1 instead of 0.694^19; the true value is 0.075, so tightening is available if the margin is ever wanted. Normality had to come out of the fixed-point chain. The target carries no such hypothesis, and non-normal operands are inside its domain in two distinct ways: exponent 0, where toReal < 2^-1022, and exponent 2047, where toReal is 0 by the Inf/NaN convention. The second wraps 1022 - 2047 and so escapes the general shift-amount reasoning entirely; it needs its own computation showing the shift still saturates. Rather than shadow the chain with primed variants, hn is removed from toNat_mtwop63 itself and mechanically from the eight lemmas that only threaded it. That is a strict strengthening. The file's own non-vacuity witness uses x = FPR.zero, which is exactly a non-normal operand, so a version carrying hn would not have applied to the example the file already contains. Proved as expm_p63_error' here; FPRBridge.lean is untouched and still carries the sorry. Relocating the statement is a decision for the file owner, alongside the split that file already needs.
The statement lived in FPRBridge.lean, which ExpmBridge.lean imports, so the proof could not close it in place. Moving it the other way lets the sorry go, and FPRBridge.lean now carries none: all six per-operation bounds are proved. The docstring keeps the domain argument, which is intrinsic to the statement, and replaces the note on how a proof would have to be shaped with how the bound actually decomposes - 10 units for the fixed-point pipeline, 3574 for the Chebyshev certificates, 80 for the Taylor truncation, against 4096. The half constant and its two lemmas move with the non-vacuity witness that uses them, which now also records that FPR.zero is not a normal operand, so the witness exercises the zero-exponent branch rather than only the ordinary one.
alik-eth
changed the title
berExp reduced x modulo log 2 with rint (round to nearest) where the reference uses fpr_trunc; on x >= 0 truncation is floor, so the remainder r belongs in [0, log 2). Rounding to nearest centres r on zero instead, putting it below zero for about half of all inputs -- 100 of the 199 points x = k/32, k < 200, and 49.8% over a uniform sweep of the range Falcon actually reaches. That is a correctness bug rather than a rounding nicety, because expm_p63 reaches its argument only through mtwop63, which reads the significand and the exponent field and never the sign bit. It therefore computes exp(-|r|), so a negative r yields exp(r) where the algorithm needs exp(-r) -- an acceptance weight wrong by a factor of exp(2r), up to 2x. The Bernoulli step is what shapes samplerZ's output distribution, so the error lands on the discrete Gaussian itself. The Lean sampler was never cross-checked against the reference: the signing tests all sign via FFI and only verify in Lean, and the one berExp test uses x = 0, where the two roundings agree. - SamplerZ: split the reduction out as berExpReduce, so the test and any future proof share one subject, and round with floor_. - Main: pin berExpReduce's output to [0, log 2] over x = k/32, k < 200. The check fails at 100 of those points under round-to-nearest. - ExpmBridge: state expm_p63_error's domain as x <= 694/1000 rather than x < log 2. 0.694 is the contraction factor the Horner error induction already runs on, and the Chebyshev certificates reach 89/128, so this costs nothing. It is needed: berExpReduce computes r by rounding a floating-point quotient, and r lands a few ulps above log 2 when x is near a multiple of it, so no statement closed at log 2 could apply. - ExpmBridge: prove the sign-blindness that makes the bug a bug -- mtwop63_neg, expm_p63_neg, and expm_p63_error_abs, the error bound over the symmetric domain against exp(-|x|). All new results are axiom-clean.
Collaborator
|
Reviewed by building the branch (3080 jobs, clean) and testing the claims rather than reading them. The central one holds up, and it is a good deal stronger than "the bounds needed tightening". VerifiedThe decoder is right. All seven bounds are axiom-clean. The refutation reproduces exactly. For So the hypothesis Three things before merge1. New prose in
That text is added by this PR and is false as of this PR. Anyone reading the file after merge concludes the bounds are unproven. The rest of that block (the 2. Worth stating the scope in the description. 3. A question about the The substantive question: the saturation logic is written for So Nothing above blocks the PR. (1) is a paragraph, (2) is description text, (3) is a question I would like answered before the sampler change lands. Agreed on the file split being its own PR — "pure code motion, checkable by comparing |
…>= 0 Review follow-ups on the FPR decoder branch. ApproxArith: the "What remains open" block claimed the four error bounds were still `sorry`. They are proved on this branch, so the claim was false as of the commit that introduced it. What actually remains is the `_valid` half of the class, which the block now says, together with the consequence worth stating plainly: with the instance still commented out, the bounds do not yet reach any downstream Falcon theorem. SamplerZ: berExp's `x >= 0` was asserted without saying where it comes from or what happens without it. Both now recorded. The saturation `s ||| (63 - s) >>> 26` detects `si >= 64` through the borrow bit and has no negative branch, so a negative `si` truncates to `63 + si` and drives the acceptance weight to zero, where the correct weight for x < 0 is above one. The precondition itself comes from |z - r| >= z0 together with sigma <= sigma0, the latter a key-generation invariant rather than anything the sampler checks. That invariant is tight in floating point, which the new tests pin from both sides: at the representable 1/sigma just above 1/sigma0, r = 0 and b = 0 give x = 0 exactly for every z0 in [0, 25]; one ulp down, sigma passes sigma0 and 25 of those 26 go negative. Testing only the positive side would not have distinguished the invariant from a coincidence. Also pins berExpReduce(0), which the k/32 sweep skipped.
Contributor
Author
|
Thanks — all three land, and (3) turned up something I had not checked. 1. 2. Scope — agreed, description updated. On 3. The And it is worse than you put it: on
What rules out
So the precondition holds for every On adding an explicit negative branch: I would rather not without your call. The reference has the identical structure and the identical unproved precondition — One thing I should flag rather than leave implicit: the end-to-end statement "the reduced argument is in range" is not proved, and it is harder than it looks. Sterbenz makes the subtraction exact, so |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The
FPRdenotation went throughFloat.ofBitsandFloat.toRat0, whose chain bottoms out on runtime primitives the kernel cannot evaluate. Nothing abouttoRealwas provable, and the five error bounds below it had satsorrysince.This decodes the IEEE-754 fields from the word instead.
Floatno longer appears in the file, andtoRealreduces in the kernel.That made the bounds refutable — and all five were false
add_error,mul_error— two largest-finite normals overflow to the non-finite exponent, which denotes zero, so the error is the entire magnitude.div_error— overflow wraps to a negative, small, finite value.sqrt_error— needs no overflow at all: at+∞the hypothesis0 ≤ toReal aholds and the right-hand side is0, whilesqrtreturns2^512.expm_p63_error— constrainedxbut said nothing aboutccs. Atccs = 1the fixed-point conversion wraps and the routine returns0against a true value ofexp(-x); above1the claim fails outright.Each is now restricted to the domain where it holds. Every conclusion is unchanged — same
2^-52, same2^-51, same shape. In-file witnesses show the hypotheses admit ordinary values, and that the result-range condition is not implied by the operand conditions.All six are now proved
add_error,sub_error,mul_error,div_error,sqrt_error,expm_p63_error— every one axiom-clean[propext, Classical.choice, Quot.sound], nobv_decide, nonative_decide.Extern/Falcon/FPRBridge.leancarries nosorry.Each kernel is modelled as a pipeline record pinned to the kernel term by
rfl, so the proofs are about the code that runs rather than a restatement of it.addandmulare straight-line.divandsqrtare not: each runs a fixed-lengthfor _ in [0:n]loop, and neither body reads the loop index, so the fold is iteration of one step function — core'sforIn_eq_forIn_range'reduces it to aList.range'fold, induction givesforIn = f^[n], and a loop invariant becomes induction on the step count.divis exact restoring division:yu·q + r = 2^55·xu, remainder below2·yu, over 55 iterations. The step count has to be bounded — past 62 the quotient's doubling overflowsUInt64.sqrtis a digit recurrence computing an exact integer square root:q² ≤ 2^54·xu < (q+1)². Its bit weight hits zero on the last iteration, so the invariant holds for 53 steps and is extended by one; that final step is what turns it into a floor. Both exponent parities then collapse to a single formula, because the odd case's doubling of the significand exactly compensates the floor in the halved exponent.mul_error_combineassemblesmul,divandsqrtalike. Writing it abstractly rather than inlining it intomul_erroris what let each later operation reuse it unchanged.Two places where the obvious bound was true but too loose
Both would have overrun the budget once composed with the
2^-53rounding step, and both needed a sharper fact rather than more care:mul's two folds cost one ulp, not two.Numbers first, lemmas second: measured against exact rational arithmetic,
div's end-to-end error peaks at 0.498 of the2^-52budget over 150k in-domain pairs. That is what said a sharper fact had to exist.The sixth: certifying a minimax polynomial
expm_p63_erroris a different kind of obligation, and it lands in a new file,Extern/Falcon/ExpmBridge.lean. Its2^-51is very nearly saturated, so the budget has to be spent carefully. In units of2^-63, against 4096:The fixed-point half is ordinary:
mulHi64is the high half of the exact product,mtwop63is exactly⌊2^63 · toReal x⌋with no rounding of its own, and the twelve Horner floors contract because each step scales the previous error byζ < 0.694. The loop'sUInt64subtraction cannot wrap for a pleasingly cheap reason:mulHi64 a b ≤ bfor everya.The certification half is where the difficulty sits.
facctCoeffsis a minimax fit, not a Taylor truncation, soP − TcarriesO(1)coefficients while beingO(2^-51)— and any bound applying a triangle inequality to those coefficients discards exactly the cancellation the fit relies on. Measured: coefficient-wise against Taylor misses by ~38000×, and midpoint-plus-Lipschitz still misses by 730× at 256 subintervals, because boundingsup |Q′|coefficient-wise has the identical pathology. Subdividing restores locality; 15 Chebyshev certificates over dyadic breakpoints land at 3574. Dyadic matters: anchoring the split at a decimal bound forlog 2produces 337-digit numerals, and the line-length linter cannot wrap a numeral.The whole certification layer costs ~25s of build time, and
docs/agents/expm-certification.mdrecords the measurements behind every choice — including the three estimates that were wrong until checked.A sampler bug the branch also fixes
berExpreducedxmodulolog 2withrint(round to nearest) where the reference usesfpr_trunc—c-fn-dsa/sign_sampler.c:874, "We can use fpr_trunc() because x >= 0". Onx ≥ 0truncation is floor, so the remainder belongs in[0, log 2); rounding to nearest centres it on zero and puts it below zero for about half of all inputs — 100 of the 199 pointsx = k/32, k < 200, and 49.8% over a uniform sweep of the range Falcon reaches.That is a correctness bug rather than a rounding nicety, because
expm_p63reaches its argument only throughmtwop63, which reads the significand and the exponent field and never bit 63. It therefore computesexp(-|r|), so a negativeryieldsexp(r)where the algorithm needsexp(-r)— an acceptance weight wrong byexp(2r), up to 2×, in the step that shapessamplerZ's output distribution. This branch proves that sign-blindness rather than asserting it:mtwop63_neg,expm_p63_neg, andexpm_p63_error_abs, the bound over the symmetric domain againstexp(-|x|).It survived because the Lean sampler is never cross-checked against the reference — every signing test signs via FFI and only verifies in Lean — and the one
berExptest usedx = 0, where the two roundings agree.The change is right precisely on
x ≥ 0, and strictly worse without it: atx = -0.01,rintlands within 2% of the correct weight by accident, whilefloor_givessShift = 63and rejects always.x ≥ 0follows from|z − r| ≥ z0together withσ ≤ σ₀, the latter a key-generation invariant rather than anything the sampler checks; in floating point its margin is exactly zero at the boundary, so the tests pin both sides —x = 0at the representable1/σjust above1/σ₀, andx < 0one ulp past it.expm_p63_error's domain moved fromx < log 2tox ≤ 694/1000for this reason: a computed reduction lands a few ulps abovelog 2whenxis near a multiple of it, so nothing stated atlog 2could ever be applied to one.0.694is the contraction factor the Horner induction already ran on, and the certificates reach89/128, so the widening cost no constant churn. The end-to-end0 ≤ ris not proved: Sterbenz makes the subtraction exact, reducing it tofl(si·L') ≤ x, butinvLog2Const · log2Const = 1 − 0.43·uwhile the naive chain must absorb(1+u)², and the fact is tight — worst margin 0 ulp, equality atsi = 1, x = fl(log 2).Scope: what this does not yet reach
HasRealSemanticsstill has no instances and no users outsideApproxArith.lean; the FPR instance stays commented out, blocked on the five_validclosure fields, which are aboutIsNormalpreservation rather than about error. So the seven bounds proved here do not yet reach any downstream Falcon theorem. Closing that, not the file split, is the more important half of what remains.Relatedly,
Extern/Falcon/Instance.lean:140still denotes a raw FPR word as(Float.ofBits x).toRat0in the idealℝmodel, whichExtern/Falcon/FFT.leanreads the GM twiddle table through — so two denotations of the same word coexist. Repointing it at the decoder should be semantically inert (they agree on every finite pattern, and both send non-finite to0), but it changes what the ideal FFT model is built on, so it belongs in its own commit.One thing to know when reviewing
The statements are textually identical to
main, buttoRealnow means the decoder. Agreement with Lean'sFloatrests on a 30-pattern differential check, not a bridge lemma: the gap is closed by replacement.Follow-up: this file is now too big
At 6.8k lines
FPRBridge.leanis nearly double the repo's previous largest (SimulateQ.lean, 3.5k);ExpmBridge.leanadds a further 1.9k. It wants splitting along the seams it already has — the layering is acyclic, so the cut is mechanical:Decode(decode,toReal, domain predicates) ·Kernel(sticky folds, leading-zero count, rounding table,make/make_z) ·Add·Mul·Loop(theforInprinciple and the two loop pipelines) ·DivSqrt· a thin re-export retainingexpm_p63_errorand the verifier bridge.Deliberately not done here: git would render a 1→6 split as thousands of deleted plus thousands of added lines and bury the content worth reviewing. Better as its own PR whose whole claim is "pure code motion" — checkable by comparing
#print axiomsfor the five bounds before and after, rather than by reading. Attribution on the new files is the file owner's call, so I would rather ask than guess.Verification
Full CI set green; CI-lint clean; extern and interop isolation OK. Every new declaration is axiom-clean — a sweep of all constants in the module found nothing else. No
bv_decideand nonative_decideanywhere: both emit custom trust axioms.Each original statement was refuted by a kernel-checked
¬ (∀ …)against its verbatim text, axiom-clean to the same standard. Those proofs are not in this branch — they quote statements it deletes — but they are the evidence that the added hypotheses are necessary rather than convenient, and I am happy to land them as a companion file if that is the more useful record.