3.2. Meta-level implementation of Bird's algorithm
3.2.1. A certified meta-level Bird determinant
This note contains an implementation of a certified Bird determinant algorithm.
open Lean Meta Qq Mathlib.Tactic Mathlib.Tactic.Ring Determinant.Cert
namespace Determinant.BirdCertified
variable
{u : Level}
{α : Q(Type u)}
{sα : Q(CommSemiring $α)}
abbrev CertMatrix
{u : Level}
{α : Q(Type u)}
(sα : Q(CommSemiring $α)) :=
Array (CertExpr sα)
private def certIdx (n i j : Nat) : Nat :=
Determinant.BirdSpec.idx n i j
/-- Read a certified matrix cell. -/
def certGet (A : CertMatrix sα) (n i j : Nat) : AtomM (CertExpr sα) := do
let some r := A[certIdx n i j]?
| throwError "certGet: index {i} {j} out of bounds"
pure r
/-- Write a certified matrix cell. -/
def certSet (A : CertMatrix sα) (n i j : Nat) (x : CertExpr sα) : CertMatrix sα :=
A.set! (certIdx n i j) x
def certEntry (e : Q($α)) : CertM sα (CertExpr sα) := do
certAtom e
def certZero : CertM sα (CertExpr sα) := certAtom q((0 : $α))
def certOne : CertM sα (CertExpr sα) := certAtom q((1 : $α))
/-- Certified sign factor `(-1) ^ (n - 1)` -/
def certSign (n : Nat) : CertM sα (CertExpr sα) := do
let env ← read
let some _ := env.cache.rα
| throwError m!"Missing CommRing instance"
let k : Q(Nat) := mkNatLit (n - 1)
certAtom q(((-1 : $α) ^ $k : $α))
/-- Sum certified terms using `certAdd` starting from certified zero -/
def certSumRange (xs : Array (CertExpr sα)) : CertM sα (CertExpr sα) := do
let mut acc ← certZero
for x in xs do
acc ← certAdd acc x
pure acc
/-- Certified matrix multiplication over `CertMatrix` cells. -/
def certMatrixMul (n : Nat) (A B : CertMatrix sα) : CertM sα (CertMatrix sα) := do
let mut out := Array.replicate (n * n) (← certZero)
for i in [0:n] do
for j in [0:n] do
let mut terms := #[]
for k in [0:n] do
terms := terms.push (← certMul (← certGet A n i k) (← certGet B n k j))
out := certSet out n i j (← certSumRange terms)
pure out
/-- Certified version of the Bird determinant μ function -/
def certMu (n : Nat) (A : CertMatrix sα) : CertM sα (CertMatrix sα) := do
let mut out := Array.replicate (n * n) (← certZero)
for i in [0:n] do
for j in [0:n] do
let x ←
if j < i then certZero
else if i = j then
let mut terms := #[]
for k in [0:n] do
terms := terms.push (← if i < k then certGet A n k k else certZero)
certNeg (← certSumRange terms)
else certGet A n i j
out := certSet out n i j x
pure out
/-- Certified iteration of Bird's recurrence. -/
partial def certIterF
(n : Nat)
(A : CertMatrix sα) :
Nat → CertMatrix sα → CertM sα (CertMatrix sα)
| 0, X => pure X
| k + 1, X => do
let X' ← certIterF n A k X
certMatrixMul n (← certMu n X') A
/-- Run Bird's determinant recurrence over certified matrix entries. -/
def birdDetCert
(n : Nat)
(A : CertMatrix sα) :
CertM sα (CertExpr sα) := do
let absDet ←
match n with
| 0 => certOne
| k + 1 => do
let B ← certIterF n A k A
certGet B n 0 0
match n with
| 0 => pure absDet
| _ + 1 => certMul (← certSign n) absDet
end Determinant.BirdCertified
3.2.2. Certified ring operations
Certified ring operations and functions to combine certificates. These are wrappers for Mathlib.Tactic.Ring.Common.Result.
open Lean Meta Qq
open Mathlib.Tactic
open Mathlib.Tactic.Ring
open Ring.Common
namespace Determinant.Cert
variable
{u : Level}
{α : Q(Type u)}
{sα : Q(CommSemiring $α)}
/--
A certified ring expression.
`raw` is the expression the caller wants to normalize. `result` is the certificate
proving that `raw` is equal to the normalized expression stored in the result.
-/
structure CertExpr
{u : Level}
{α : Q(Type u)}
(sα : Q(CommSemiring $α)) where
raw : Q($α)
result : Result (Common.ExSum RatCoeff sα) raw
-- TODO: Rename to Context
/-- An environment for the certificate run. -/
structure Env
{u : Level}
{α : Q(Type u)}
(sα : Q(CommSemiring $α)) where
cache: Cache sα
rc: RingCompute RatCoeff sα
def mkEnv : MetaM (Env sα) := do
let cache ← mkCache sα
pure {
cache
rc := ringCompute cache
}
/-- Construct an `Env` with a prepopulated Cache` -/
def mkEnvCache (cache : Cache sα) : MetaM (Env sα) := do
pure {
cache
rc := ringCompute cache
}
/-- The monad used in certificate producing functions. -/
abbrev CertM
{u : Level}
{α : Q(Type u)}
(sα : Q(CommSemiring $α))
:= ReaderT (Env sα) AtomM
/-- Run the certificate computation with a given `Env` and then run the
underlying `AtomM`. -/
def CertM.runEnv {β : Type} (cert : CertM sα β) (env : Env sα) : MetaM β :=
AtomM.run .reducible (cert.run env)
namespace CertExpr
/-- The normalized expression represented by a `CertExpr`. -/
def norm (x : CertExpr sα) : Q($α) :=
x.result.expr
/-- The `ExBase` normal-form stored in a `CertExpr`. -/
def nf (x : CertExpr sα) : Ring.ExSum sα x.norm :=
x.result.val
/-- The proof that the raw expression equals the normalized expression. -/
def proof (x : CertExpr sα) : Q($x.raw = $x.norm) :=
x.result.proof
end CertExpr
/-- Normalize an expression and package the normal-form certificate as a
`CertExpr` -/
def normalizeExpr (e : Q($α)) : CertM sα (CertExpr sα) := do
let env ← read
let result ← eval rcℕ env.rc env.cache e
pure {
raw := e
result
}
/-- Normalize an expression -/
def certAtom (e : Q($α)) : CertM sα (CertExpr sα) :=
normalizeExpr e
/-- Combine two certificates into a certificate for the sum of their raw
expressions. -/
def certAdd (x y : CertExpr sα) : CertM sα (CertExpr sα) := do
let env ← read
let r ← evalAdd env.rc rcℕ x.nf y.nf
pure {
raw := q($x.raw + $y.raw)
result := {
expr := r.expr
val := r.val
proof := q(add_congr $x.proof $y.proof $r.proof)
}
}
/-- Combine two certificates into a certificate for the product of their raw
expressions. -/
def certMul (x y : CertExpr sα) : CertM sα (CertExpr sα) := do
let env ← read
let r ← evalMul env.rc rcℕ x.nf y.nf
pure {
raw := q($x.raw * $y.raw)
result := {
expr := r.expr
val := r.val
proof := q(mul_congr $x.proof $y.proof $r.proof)
}
}
/-- Combine two certificates into a certificate for the difference of their raw
expressions. -/
def certSub (x y : CertExpr sα) : CertM sα (CertExpr sα) := do
let env ← read
let some rα := env.cache.rα
| throwError m!"Missing CommRing instance"
let r ← evalSub env.rc rcℕ rα x.nf y.nf
pure {
raw := q($x.raw - $y.raw)
result := {
expr := r.expr
val := r.val
proof := q(sub_congr $x.proof $y.proof $r.proof)
}
}
/-- Produce a certificate for the negation of a raw expression -/
def certNeg (x : CertExpr sα) : CertM sα (CertExpr sα) := do
let env ← read
let some rα := env.cache.rα
| throwError m!"Missing CommRing instance"
let r ← evalNeg env.rc rα x.nf
pure {
raw := q(-$x.raw)
result := {
expr := r.expr
val := r.val
proof := q(neg_congr $x.proof $r.proof)
}
}
/-- Recursively build a certificate for a ring expression.
The parser recognizes addition, multiplication, negation and subtractions and
uses the corresponding certificate combinator. Anything else is treated as an
atom.
-/
partial def certExpr (e : Q($α)) : CertM sα (CertExpr sα) := do
let env ← read
let some _ := env.cache.rα | throwError "missing CommRing instance"
match e with
| ~q($a + $b) => certAdd (← certExpr a) (← certExpr b)
| ~q($a * $b) => certMul (← certExpr a) (← certExpr b)
| ~q(-$a) => certNeg (← certExpr a)
| ~q($a - $b) => certSub (← certExpr a) (← certExpr b)
| _ => certAtom e
/-- Turn two expressions with the same normal form into an equality of the
original expressions. -/
theorem ofSameNormalizedEq
{R : Type*}
{lhs rhs norm : R}
(hl : lhs = norm)
(hr : rhs = norm) : lhs = rhs := R:Type u_1lhs:Rrhs:Rnorm:Rhl:lhs = normhr:rhs = norm⊢ lhs = rhs
All goals completed! 🐙
/-- Prove equality of two ring expressions by normalizing both sides in a
shared AtomM context. -/
def proveEq (cache : Cache sα) (lhs rhs : Q($α)) : MetaM Q($lhs = $rhs) := do
let env ← mkEnvCache cache
let mkCert : CertM sα (CertExpr sα × CertExpr sα) := do
let lhsCert ← certExpr lhs
let rhsCert ← certExpr rhs
pure (lhsCert, rhsCert)
let (lhsCert, rhsCert) ← mkCert.runEnv env
unless lhsCert.nf.eq rcℕ env.rc rhsCert.nf do
throwError "ring expressions not equal"
let lhsProof : Q($lhs = $lhsCert.norm) := lhsCert.proof
let rhsProof : Q($rhs = $rhsCert.norm) := rhsCert.proof
have : $lhsCert.norm =Q $rhsCert.norm := ⟨⟩
pure q(ofSameNormalizedEq $lhsProof $rhsProof)
def proveCertEq (lhsCert rhsCert : CertExpr sα) : CertM sα Q($lhsCert.raw = $rhsCert.raw) := do
unless lhsCert.nf.eq rcℕ (← read).rc rhsCert.nf do
throwError m!"proveCertEq not equal lhs: {lhsCert.norm} rhs: {rhsCert.norm}"
let lhsProof : Q($lhsCert.raw = $lhsCert.norm) := lhsCert.proof
let rhsProof : Q($rhsCert.raw = $rhsCert.norm) := rhsCert.proof
have : $lhsCert.norm =Q $rhsCert.norm := ⟨⟩
pure q(ofSameNormalizedEq $lhsProof $rhsProof)
end Determinant.Cert
3.2.3. Bird implementation using a FlatMatrix
An implementation of the Bird determinant using a single Array, with entries in row-major order.
namespace Determinant.BirdSpec
-- set_option trace.compiler.ir.result true
/-- Dense Matrix representation.
Entries are stored in row-major order in an `Array`.
-/
abbrev FlatMatrix (R : Type*) :=
Array R
/-- Row-major index for the `(i, j)` entry of an `n` by `n` `FlatMatrix -/
def idx (n i j : Nat) : Nat :=
i * n + j
-- TODO: Try defining a GetElem instance
/-- `FlatMatrix` lookup. Returns `0` when the index is out of bounds -/
def get {R : Type*} [Zero R] (n : Nat) (A : FlatMatrix R) (i j : Nat) : R :=
A.getD (idx n i j) 0
/-- Sum `f 0 + ... + f (n - 1)`. -/
def sumRange {R : Type*} [Zero R] [Add R] (n : Nat) (f : Nat → R) : R :=
(List.range n).foldl (fun acc k => acc + f k) 0
/-- Build an `n` by `n` `FlatMatrix` from a function of row / column indices -/
def ofFn {R : Type*} (n : Nat) (f : Nat → Nat → R) : FlatMatrix R :=
((List.range (n * n)).map fun p => f (p / n) (p % n)).toArray
/-- Matrix mulitplication -/
def mul {R : Type*} [CommRing R]
(n : Nat) (A B : FlatMatrix R) : FlatMatrix R :=
ofFn n fun i j =>
sumRange n fun k => get n A i k * get n B k j
/-- The Bird algorithm μ function applied to a `FlatMatrix` -/
def mu {R : Type*} [CommRing R]
(n : Nat) (A : FlatMatrix R) : FlatMatrix R :=
ofFn n fun i j =>
if j < i then 0
else if i = j then
-sumRange n fun k => if i < k then get n A k k else 0
else
get n A i j
/-- The Bird algorithm F^{n - 1} function. -/
def iterF {R : Type*} [CommRing R]
(n : Nat) (A : FlatMatrix R) :
Nat → FlatMatrix R → FlatMatrix R
| 0, X => X
| k + 1, X => mul n (mu n (iterF n A k X)) A
/-- Computation of the determinant of a `FlatMatrix` using Bird's algorithm. -/
def birdDetFlat {R : Type*} [CommRing R]
(n : Nat) (A : FlatMatrix R) : R :=
match n with
| 0 => 1
| k + 1 => (-1 : R) ^ (n - 1) * get n (iterF n A k A) 0 0
end Determinant.BirdSpec
3.2.4. A Bird determinant tactic specialized for 2x2 matrices
Certified ring operations and functions to combine certificates. These are wrappers for Mathlib.Tactic.Ring.Common.Result.
open Lean Elab Tactic Meta Qq Mathlib.Tactic.Ring.Common Determinant.Cert Determinant.BirdCertified
namespace RingCert
def checkBird2x2Eq (aStx bStx cStx dStx rhsStx : Syntax) : TacticM Unit := do
withMainContext do
let aExpr ← elabTerm aStx none
let αExpr ← instantiateMVars (← inferType aExpr)
let bExpr ← elabTerm bStx (some αExpr)
let cExpr ← elabTerm cStx (some αExpr)
let dExpr ← elabTerm dStx (some αExpr)
let rhsExpr ← elabTerm rhsStx (some αExpr)
let ⟨_, α⟩ ← getLevelQ' αExpr
have a : Q($α) := aExpr
have b : Q($α) := bExpr
have c : Q($α) := cExpr
have d : Q($α) := dExpr
have rhs : Q($α) := rhsExpr
let sα : Q(CommSemiring $α) ← synthInstanceQ q(CommSemiring $α)
let cache ← mkCache sα
let env ← mkEnvCache cache
let action : CertM sα (CertExpr sα × CertExpr sα) := do
let aCert ← certAtom a
let bCert ← certAtom b
let cCert ← certAtom c
let dCert ← certAtom d
let A : CertMatrix sα := #[aCert, bCert, cCert, dCert]
let detCert ← birdDetCert 2 A
let rhsCert ← certExpr rhs
logInfo m!"{detCert.raw} = {rhsCert.raw}"
logInfo m!"{detCert.norm} = {rhsCert.norm}"
logInfo m!"{← inferType detCert.proof} = {← inferType rhsCert.proof}"
pure (detCert, rhsCert)
let check : CertM sα Unit := do
let (detCert, rhsCert) ← action
discard <| proveCertEq detCert rhsCert
check.runEnv env
syntax "check_bird_2x2 " term:arg term:arg term:arg term:arg " = " term : tactic
elab_rules : tactic
| `(tactic| check_bird_2x2 $a:term $b:term $c:term $d:term = $rhs:term) =>
checkBird2x2Eq a b c d rhs
3.2.5. Examples of using the checkbird2x2 assertion
set_option linter.unusedTactic false
set_option linter.unusedVariables false
theorem birdCert_test1 {R : Type*} [CommRing R] (a b c d : R) : True := R:Type u_1inst✝:CommRing Ra:Rb:Rc:Rd:R⊢ True
R:Type u_1inst✝:CommRing Ra:Rb:Rc:Rd:R⊢ True
All goals completed! 🐙
theorem birdCert_test2 {R : Type*} [CommRing R] (a b c d : R) : True := R:Type u_1inst✝:CommRing Ra:Rb:Rc:Rd:R⊢ True
R:Type u_1inst✝:CommRing Ra:Rb:Rc:Rd:R⊢ True
All goals completed! 🐙
theorem birdCert_test3 : True := ⊢ True
⊢ True
All goals completed! 🐙