Grace/Cumulus/LinearAlgebra/Elimination.g17

Open in Workbench Download source

Imports

Matrix · Rat.Computable

Declarations

Source

1-- Executable exact Gaussian elimination on finite rational arrays.2-- A cell (a,b,d) denotes (a-b)/d; inputs use d > 0.3-- Part of Cumulus, the Grace17 mathematics library.4-- Copyright (c) 2026 Grant Passmore.5-- Released under the Apache License 2.0. See LICENSE in the repository root.6-- Authors: Grant Passmore78module Grace.Cumulus.LinearAlgebra.Elimination9import Matrix10import Rat.Computable11-- Nested executable containers use Set until their carrier sets are constructed.12-- Constructor equations and mathematical correctness certificates remain checked.1314namespace RationalMatrixOps1516-- A missing coordinate reads as zero; inputs for elimination are rectangular.17def head (row : Set) → Set18  | List.nil => RationalCellOps.zero19  | List.cons(x, xs) => x2021def drop (row : Set) → Set22  | List.nil => List.nil23  | List.cons(x, xs) => xs2425def entry (row : Set, column : Nat) → Set26  | row, zero => head(row)27  | row, succ(column) => entry(drop(row), column)2829def scale (row : Set, factor : RationalCell(ℕ)) → Set30  | List.nil, factor => List.nil31  | List.cons(x, xs), factor =>32    List.cons(RationalCellOps.cellMul(factor, x), scale(xs, factor))3334def addScaled (target : Set, source : Set, factor : RationalCell(ℕ)) → Set35  | List.nil, source, factor => List.nil36  | List.cons(x, xs), source, factor =>37    List.cons(RationalCellOps.cellAdd(x, RationalCellOps.cellMul(factor, head(source))),38      addScaled(xs, drop(source), factor))3940def firstRow (rows : Set) → Set41  | List.nil => List.nil42  | List.cons(row, rest) => row4344def restRows (rows : Set) → Set45  | List.nil => List.nil46  | List.cons(row, rest) => rest4748def insertSecond (row : Set, result : Set) → Set49  | row, Option.none => Option.none50  | row, Option.some(picked) =>51    Option.some(List.cons(firstRow(picked), List.cons(row, restRows(picked))))5253def selectPivot (zero : Bool, row : Set, rest : Set, candidate : Set) → Set54  | Bool.false, row, rest, candidate => Option.some(List.cons(row, rest))55  | Bool.true, row, rest, candidate => insertSecond(row, candidate)5657-- Search is structural recursion on rows, preserving their relative order.58def promote (rows : Set, column : Nat) → Set59  | List.nil, column => Option.none60  | List.cons(row, rest), column =>61    selectPivot(RationalCellOps.isZero(entry(row, column)), row, rest, promote(rest, column))6263def clearBelow (rows : Set, pivot : Set, column : Nat) → Set64  | List.nil, pivot, column => List.nil65  | List.cons(row, rest), pivot, column =>66    List.cons(addScaled(row, pivot, RationalCellOps.cellNeg(entry(row, column))),67      clearBelow(rest, pivot, column))6869def normalizedPivot (candidate : Set, column : Nat) → Set70  | Option.none, column => Option.none71  | Option.some(picked), column =>72    Option.some(scale(firstRow(picked), RationalCellOps.inverse(entry(firstRow(picked), column))))7374def remaining (candidate : Set, original : Set, column : Nat) → Set75  | Option.none, original, column => original76  | Option.some(picked), original, column =>77    clearBelow(restRows(picked),78      scale(firstRow(picked), RationalCellOps.inverse(entry(firstRow(picked), column))), column)7980def attach (pivot : Set, rest : Set) → Set81  | Option.none, rest => rest82  | Option.some(row), rest => List.cons(row, rest)8384-- Structural recursion on remaining columns. A pivot consumes one row;85-- an empty pivot column advances without changing the row list.86def forward (columns column : Nat, rows : Set) → Set87  | zero, column, rows => rows88  | succ(columns), column, rows =>89    attach(normalizedPivot(promote(rows, column), column),90      forward(columns, succ(column), remaining(promote(rows, column), rows, column)))9192-- Back substitution uses each normalized lower pivot row to clear its93-- column in the rows above. The search ignores the augmented right side.94def choosePivotIndex (zero : Bool, current : Nat, later : Option(ℕ)) → Option(ℕ)95  | Bool.false, current, later => Option.some(current)96  | Bool.true, current, later => later9798def pivotIndex (row : Set, columns current : Nat) → Set99  | row, zero, current => Option.none100  | row, succ(columns), current =>101    choosePivotIndex(RationalCellOps.isZero(entry(row, current)), current,102      pivotIndex(row, columns, succ(current)))103104def clearAt (pivot : Option(ℕ), target source : Set) → Set105  | Option.none, target, source => target106  | Option.some(column), target, source =>107    addScaled(target, source, RationalCellOps.cellNeg(entry(target, column)))108109def clearUsing (pivots : Set, target : Set, columns : Nat) → Set110  | List.nil, target, columns => target111  | List.cons(pivot, rest), target, columns =>112    clearUsing(rest, clearAt(pivotIndex(pivot, columns, 0), target, pivot), columns)113114def prependReduced (row : Set, tail : Set, columns : Nat) : Set :=115  List.cons(clearUsing(tail, row, columns), tail)116117def backward (rows : Set, columns : Nat) → Set118  | List.nil, columns => List.nil119  | List.cons(row, rest), columns => prependReduced(row, backward(rest, columns), columns)120121def reducedEchelon (columns : Nat, rows : Set) : Set :=122  backward(forward(columns, 0, rows), columns)123124-- Reject malformed fraction data before executing a pivot. This check is125-- computational; the proof side uses the corresponding denominator premises.126def validRow (row : Set) → Set127  | List.nil => Bool.true128  | List.cons(x, xs) => BoolOps.and(RationalCellOps.valid(x), validRow(xs))129130def validRows (rows : Set) → Set131  | List.nil => Bool.true132  | List.cons(row, rest) => BoolOps.and(validRow(row), validRows(rest))133134def guardForward (valid : Bool, columns column : Nat, rows : Set) → Set135  | Bool.false, columns, column, rows => Option.none136  | Bool.true, columns, column, rows => Option.some(forward(columns, column, rows))137138def forwardChecked (columns column : Nat, rows : Set) : Set :=139  guardForward(validRows(rows), columns, column, rows)140141def guardReduced (valid : Bool, columns : Nat, rows : Set) → Set142  | Bool.false, columns, rows => Option.none143  | Bool.true, columns, rows => Option.some(reducedEchelon(columns, rows))144145def reducedEchelonChecked (columns : Nat, rows : Set) : Set :=146  guardReduced(validRows(rows), columns, rows)147148end RationalMatrixOps149150namespace RationalEliminationExamples151152def rowOne : Set :=153  List.cons(RationalCell.fraction(1, 0, 1),154    List.cons(RationalCell.fraction(1, 0, 1),155      List.cons(RationalCell.fraction(3, 0, 1), List.nil)))156def rowTwo : Set :=157  List.cons(RationalCell.fraction(1, 0, 1),158    List.cons(RationalCell.fraction(0, 1, 1),159      List.cons(RationalCell.fraction(1, 0, 1), List.nil)))160def system : Set := List.cons(rowOne, List.cons(rowTwo, List.nil))161162def reduced : Set := RationalMatrixOps.forward(2, 0, system)163164theorem pivot_exists165  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalMatrixOps.firstRow(reduced), 0)) = Bool.false166proof167  compute168169theorem eliminated_below170  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalMatrixOps.firstRow(RationalMatrixOps.restRows(reduced)), 0)) = Bool.true171proof172  compute173174theorem second_pivot_exists175  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalMatrixOps.firstRow(RationalMatrixOps.restRows(reduced)), 1)) = Bool.false176proof177  compute178179-- A missing leading pivot is swapped in; a negative leading coefficient is180-- inverted with the correct sign. These are executable regression proofs.181def swapped : Set := List.cons(182  List.cons(RationalCell.fraction(0, 0, 1), List.cons(RationalCell.fraction(1, 0, 1), List.nil)),183  List.cons(List.cons(RationalCell.fraction(0, 1, 1), List.cons(RationalCell.fraction(1, 0, 1), List.nil)), List.nil))184185theorem swapped_pivot_is_one186  shows RationalMatrixOps.entry(RationalMatrixOps.firstRow(RationalMatrixOps.forward(1, 0, swapped)), 0) =187    RationalCell.fraction(1, 0, 1)188proof189  compute190191theorem skipped_column192  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalMatrixOps.firstRow(193    RationalMatrixOps.forward(2, 0, List.cons(194      List.cons(RationalCellOps.zero, List.cons(RationalCellOps.one, List.nil)), List.nil))), 1)) = Bool.false195proof196  compute197198theorem malformed_denominator_rejected199  shows RationalMatrixOps.forwardChecked(1, 0,200    List.cons(List.cons(RationalCell.fraction(1, 0, 0), List.nil), List.nil)) = Option.none201proof202  compute203204theorem example_is_accepted205  shows RationalMatrixOps.forwardChecked(2, 0, system) = Option.some(reduced)206proof207  compute208209def completelyReduced : Set :=210  RationalMatrixOps.reducedEchelon(2, system)211212theorem cleared_above_second_pivot213  shows RationalCellOps.isZero(RationalMatrixOps.entry(214    RationalMatrixOps.firstRow(completelyReduced), 1)) = Bool.true215proof216  compute217218theorem left_pivot_remains219  shows RationalCellOps.isZero(RationalMatrixOps.entry(220    RationalMatrixOps.firstRow(completelyReduced), 0)) = Bool.false221proof222  compute223224theorem reduced_system_is_accepted225  shows RationalMatrixOps.reducedEchelonChecked(2, system) = Option.some(completelyReduced)226proof227  compute228229theorem solution_x_is_two230  shows RationalCellOps.sameValue(RationalMatrixOps.entry(231    RationalMatrixOps.firstRow(completelyReduced), 2),232      RationalCell.fraction(2, 0, 1)) = Bool.true233proof234  compute235236theorem solved_second_coordinate237  shows RationalMatrixOps.entry(RationalMatrixOps.firstRow(RationalMatrixOps.restRows(reduced)), 2) =238    RationalCell.fraction(3, 1, 2)239proof240  compute241242eval RationalCellOps.isZero(RationalMatrixOps.entry(RationalMatrixOps.firstRow(RationalMatrixOps.restRows(reduced)), 0))243244end RationalEliminationExamples245246eval RationalEliminationExamples.reduced247248eval RationalEliminationExamples.completelyReduced249250-- Finite row-operation certificates. Their constructors are computational251-- data; replay is structural recursion over the certificate list.252datatype RationalRowOperation(Index : Set)253  | swap(first second : Index)254  | scale(row : Index, factor : Set)255  | addMultiple(target source : Index, factor : Set)256257namespace RationalRowCertificates258259def replaceHead (rows : Set, value : Set) → Set260  | List.nil, value => List.nil261  | List.cons(row, rest), value => List.cons(value, rest)262263def replaceTail (rows rest : Set) → Set264  | List.nil, rest => List.nil265  | List.cons(row, tail), rest => List.cons(row, rest)266267def rowAt (rows : Set, index : Nat) → Set268  | rows, zero => RationalMatrixOps.firstRow(rows)269  | rows, succ(index) => rowAt(RationalMatrixOps.restRows(rows), index)270271def replaceRow (rows : Set, index : Nat, value : Set) → Set272  | rows, zero, value => replaceHead(rows, value)273  | rows, succ(index), value =>274    replaceTail(rows, replaceRow(RationalMatrixOps.restRows(rows), index, value))275276def swapRows (rows : Set, first second : Nat) : Set :=277  replaceRow(replaceRow(rows, first, rowAt(rows, second)), second, rowAt(rows, first))278279def scaleRow (rows : Set, row : Nat, factor : RationalCell(ℕ)) : Set :=280  replaceRow(rows, row, RationalMatrixOps.scale(rowAt(rows, row), factor))281282def addRow (rows : Set, target source : Nat, factor : RationalCell(ℕ)) : Set :=283  replaceRow(rows, target,284    RationalMatrixOps.addScaled(rowAt(rows, target), rowAt(rows, source), factor))285286def execute (rows : Set, operation : RationalRowOperation(ℕ)) → Set287  | rows, RationalRowOperation.swap(first, second) => swapRows(rows, first, second)288  | rows, RationalRowOperation.scale(row, factor) => scaleRow(rows, row, factor)289  | rows, RationalRowOperation.addMultiple(target, source, factor) =>290    addRow(rows, target, source, factor)291292def replay (operations : Set, rows : Set) → Set293  | List.nil, rows => rows294  | List.cons(step, rest), rows => replay(rest, execute(rows, step))295296-- The factor fields are sets so the datatype can store a mix of swaps and297-- rational operations. Validate their rational representation before replay.298def validStep (operation : RationalRowOperation(ℕ)) → Set299  | RationalRowOperation.swap(first, second) => Bool.true300  | RationalRowOperation.scale(row, factor) =>301    BoolOps.and(RationalCellOps.valid(factor), BoolOps.not(RationalCellOps.isZero(factor)))302  | RationalRowOperation.addMultiple(target, source, factor) => RationalCellOps.valid(factor)303304def validTrace (operations : Set) → Set305  | List.nil => Bool.true306  | List.cons(step, rest) => BoolOps.and(validStep(step), validTrace(rest))307308def inBounds (operation : RationalRowOperation(ℕ), count : Nat) → Bool309  | RationalRowOperation.swap(first, second), count =>310    BoolOps.and(RationalCellOps.less(first, count), RationalCellOps.less(second, count))311  | RationalRowOperation.scale(row, factor), count => RationalCellOps.less(row, count)312  | RationalRowOperation.addMultiple(target, source, factor), count =>313    BoolOps.and(RationalCellOps.less(target, count), RationalCellOps.less(source, count))314315def boundedTrace (operations : Set, count : Nat) → Set316  | List.nil, count => Bool.true317  | List.cons(step, rest), count =>318    BoolOps.and(inBounds(step, count), boundedTrace(rest, count))319320def rowWidth (row : Set) → Set321  | List.nil => 0322  | List.cons(cell, rest) => succ(rowWidth(rest))323324def rowCount (rows : Set) → Set325  | List.nil => 0326  | List.cons(row, rest) => succ(rowCount(rest))327328def rectangular (rows : Set, width : Nat) → Set329  | List.nil, width => Bool.true330  | List.cons(row, rest), width =>331    BoolOps.and(equalNat(rowWidth(row), width), rectangular(rest, width))332333def guardReplay (accepted : Bool, operations : Set, rows : Set) → Set334  | Bool.false, operations, rows => Option.none335  | Bool.true, operations, rows => Option.some(replay(operations, rows))336337def replayChecked (operations : Set, rows : Set) : Set :=338  guardReplay(BoolOps.and(RationalMatrixOps.validRows(rows),339    BoolOps.and(rectangular(rows, rowWidth(RationalMatrixOps.firstRow(rows))),340      BoolOps.and(validTrace(operations), boundedTrace(operations, rowCount(rows))))),341    operations, rows)342343end RationalRowCertificates344345namespace RationalRowSolutions346347-- Adding a multiple of an established equation is reversible: subtracting348-- that same multiple recovers the original equation in quotient ℚ.349theorem addMultiple.equation_iff (p q r s t : Rational)350  requires first : p = q351  shows (r = s) ↔ (r + t * p = s + t * q)352proof353  constructor354  . intro second355    rw [first, second]356    rfl357  . intro changed358    rw first at changed359    have leftCancel := rationalAddValue.add_sub_cancel_right at r, t * q360    have rightCancel := rationalAddValue.add_sub_cancel_right at s, t * q361    have cancel : (r + t * q) - t * q = (s + t * q) - t * q362      rw changed363      rfl364    rwa [leftCancel, rightCancel] at cancel365366theorem addMultiple.two_variable_solutions (a b c d e f t x y : Rational)367  requires first : a * x + b * y = c368  shows (d * x + e * y = f) ↔369    ((d + t * a) * x + (e + t * b) * y = f + t * c)370proof371  have shape : (d + t * a) * x + (e + t * b) * y =372    (d * x + e * y) + t * (a * x + b * y)373    ring374  rw shape375  exact addMultiple.equation_iff with first376377end RationalRowSolutions378379namespace GeneralRowSolutions380381-- A row is a real sequence. The finite dot product is already defined by382-- recursion on its width; its add and scale laws were proved by induction.383def addMultipleRow (target source : Set, C : Real) : Set :=384  addSequence(target, mulSequence(constantSequence(C), source))385386theorem addMultipleRow.sequence (target source : Set, C : Real)387  requires targetSequence : RealSequence(target)388  requires sourceSequence : RealSequence(source)389  shows RealSequence(addMultipleRow(target, source, C))390proof391  have constant := constantSequence.sequence at C392  have scaled := mulSequence.sequence with constant, sourceSequence393  unfold addMultipleRow394  exact addSequence.sequence with targetSequence, scaled395396theorem addMultipleRow.dot (target source vector : Set, C : Real, n : Nat)397  requires targetSequence : RealSequence(target)398  requires sourceSequence : RealSequence(source)399  requires vectorSequence : RealSequence(vector)400  shows dotProduct(addMultipleRow(target, source, C), vector, n) =401    (dotProduct(target, vector, n):R) + C * (dotProduct(source, vector, n):R)402proof403  have constant := constantSequence.sequence at C404  have scaled := mulSequence.sequence with constant, sourceSequence405  have added := dotProduct.add_left with targetSequence, scaled, vectorSequence406  have product := dotProduct.scale_left with sourceSequence, vectorSequence407  unfold addMultipleRow408  simp [added, product]409410theorem addMultipleRow.solution_iff (target source vector : Set, C sourceValue targetValue : Real, n : Nat)411  requires targetSequence : RealSequence(target)412  requires sourceSequence : RealSequence(source)413  requires vectorSequence : RealSequence(vector)414  requires first : dotProduct(source, vector, n) = sourceValue415  shows (dotProduct(target, vector, n) = targetValue) ↔416    (dotProduct(addMultipleRow(target, source, C), vector, n) =417      targetValue + C * sourceValue)418proof419  have shape := addMultipleRow.dot at target, source, vector, C, n with420    targetSequence, sourceSequence, vectorSequence421  have targetReal := dotProduct.real with targetSequence, vectorSequence422  have productReal := realMulValue.real at C, sourceValue with423    *{2}424  constructor425  . intro second426    rw [shape, first, second]427    rfl428  . intro changed429    rw [shape, first] at changed430    have leftCancel := realSubValue.add_cancel_right with targetReal, productReal431    have rightCancel := realSubValue.add_cancel_right at432      targetValue, C * sourceValue with assumption, productReal433    have cancel : ((dotProduct(target, vector, n):R) + C * sourceValue) - C * sourceValue =434      (targetValue + C * sourceValue) - C * sourceValue435      rw changed436      rfl437    rwa [leftCancel, rightCancel] at cancel438439predicate SolvesMatrix (M b u : Set, m n : Nat) :=440  ∀ i ∈ ℕ, i < m → dotProduct(matrixRow(M, i), u, n) = b[i]441442-- This theorem is independent of the number of rows and columns. A matrix443-- update is specified by its changed row and its unchanged entries; the444-- executable update can discharge those equations separately.445theorem addMultiple.matrix_solutions (M N b c u : Set, m n source target : Nat, C : Real)446  requires matrixM : RealMatrix(M)447  requires matrixN : RealMatrix(N)448  requires vector : RealSequence(u)449  requires rhsB : RealSequence(b)450  requires rhsC : RealSequence(c)451  requires sourceBound : source < m452  requires targetBound : target < m453  requires different : source ≠ target454  requires changedRow : matrixRow(N, target) =455    addMultipleRow(matrixRow(M, target), matrixRow(M, source), C)456  requires changedValue : c[target] = (b[target]:R) + C * (b[source]:R)457  requires otherRows : ∀ i ∈ ℕ, i < m → i ≠ target →458    matrixRow(N, i) = matrixRow(M, i) ∧ c[i] = b[i]459  shows SolvesMatrix(M, b, u, m, n) ↔ SolvesMatrix(N, c, u, m, n)460proof461  have sourceSequence := matrixRow.sequence at M, source with matrixM462  have targetSequence := matrixRow.sequence at M, target with matrixM463  have sourceNat := nat.omega.forward at source with assumption464  have targetNat := nat.omega.forward at target with assumption465  have differentSource : source ≠ target by exact different466  obtain ⟨sourceRow, sourceValue⟩ := otherRows with source, sourceNat, sourceBound, differentSource467  have preserve := addMultipleRow.solution_iff at468    matrixRow(M, target), matrixRow(M, source), u, C,469    b[source], b[target], n with470    targetSequence, sourceSequence, vector471  constructor472  . intro old473    intros i naturalI within474    have alternatives : i = target ∨ i ≠ target by tauto475    cases alternatives with equal | apart476    . rw equal477      have sourceEquation := old with source, sourceNat, sourceBound478      have targetEquation := old with target, targetNat, targetBound479      have transformed := preserve with sourceEquation480      have result := transformed.forward with targetEquation481      rw [changedRow, changedValue]482      exact result483    . have unchanged := otherRows with i, naturalI, within, apart484      rcases unchanged with ⟨rowEqual, valueEqual⟩485      rw [rowEqual, valueEqual]486      exact old with i, naturalI, within487  . intro new488    intros i naturalI within489    have alternatives : i = target ∨ i ≠ target by tauto490    cases alternatives with equal | apart491    . rw equal492      have sourceEquation := new with source, sourceNat, sourceBound493      rw [sourceRow, sourceValue] at sourceEquation494      have targetEquation := new with target, targetNat, targetBound495      rw [changedRow, changedValue] at targetEquation496      have transformed := preserve with sourceEquation497      exact transformed.backward with targetEquation498    . have unchanged := otherRows with i, naturalI, within, apart499      rcases unchanged with ⟨rowEqual, valueEqual⟩500      have equation := new with i, naturalI, within501      rw [rowEqual, valueEqual] at equation502      exact equation503504end GeneralRowSolutions505506namespace RationalRowTraceExamples507508def cell (n : Nat) : RationalCell(ℕ) := RationalCell.fraction(n, 0, 1)509def minusOne : RationalCell(ℕ) := RationalCell.fraction(0, 1, 1)510511-- y+z=5, x+y+z=6, x+2y+3z=14. The first pivot needs a swap.512def threeRows : Set :=513  List.cons(List.cons(cell(0), List.cons(cell(1), List.cons(cell(1), List.cons(cell(5), List.nil)))),514  List.cons(List.cons(cell(1), List.cons(cell(1), List.cons(cell(1), List.cons(cell(6), List.nil)))),515  List.cons(List.cons(cell(1), List.cons(cell(2), List.cons(cell(3), List.cons(cell(14), List.nil)))), List.nil)))516517def threeTrace : Set :=518  List.cons(RationalRowOperation.swap(0, 1),519  List.cons(RationalRowOperation.addMultiple(2, 0, minusOne),520  List.cons(RationalRowOperation.addMultiple(2, 1, minusOne),521  List.cons(RationalRowOperation.addMultiple(1, 2, minusOne),522  List.cons(RationalRowOperation.addMultiple(0, 2, minusOne),523  List.cons(RationalRowOperation.addMultiple(0, 1, minusOne), List.nil))))))524525def threeResult : Set := RationalRowCertificates.replay(threeTrace, threeRows)526527def hasFourValues (row : Set, a b c d : RationalCell(ℕ)) : Set :=528  BoolOps.and(RationalCellOps.sameValue(RationalMatrixOps.entry(row, 0), a),529    BoolOps.and(RationalCellOps.sameValue(RationalMatrixOps.entry(row, 1), b),530      BoolOps.and(RationalCellOps.sameValue(RationalMatrixOps.entry(row, 2), c),531        RationalCellOps.sameValue(RationalMatrixOps.entry(row, 3), d))))532533theorem threeTrace.checked534  shows RationalRowCertificates.replayChecked(threeTrace, threeRows) = Option.some(threeResult)535proof536  compute537538theorem threeTrace.reducedIdentity539  shows BoolOps.and(hasFourValues(RationalRowCertificates.rowAt(threeResult, 0), cell(1), cell(0), cell(0), cell(1)),540    BoolOps.and(hasFourValues(RationalRowCertificates.rowAt(threeResult, 1), cell(0), cell(1), cell(0), cell(2)),541      hasFourValues(RationalRowCertificates.rowAt(threeResult, 2), cell(0), cell(0), cell(1), cell(3)))) = Bool.true542proof543  compute544545theorem threeTrace.x546  shows RationalCellOps.sameValue(RationalMatrixOps.entry(RationalRowCertificates.rowAt(threeResult, 0), 3), cell(1)) = Bool.true547proof548  compute549550theorem threeTrace.y551  shows RationalCellOps.sameValue(RationalMatrixOps.entry(RationalRowCertificates.rowAt(threeResult, 1), 3), cell(2)) = Bool.true552proof553  compute554555theorem threeTrace.z556  shows RationalCellOps.sameValue(RationalMatrixOps.entry(RationalRowCertificates.rowAt(threeResult, 2), 3), cell(3)) = Bool.true557proof558  compute559560-- The reduced right-hand side is an equality in the set-theoretic rationals;561-- the executable fraction comparison is lifted to that meaning.562theorem threeTrace.quotient_x563  shows RationalCellOps.meaning(RationalMatrixOps.entry(564    RationalRowCertificates.rowAt(threeResult, 0), 3)) = rationalClass(1, 0, 1)565proof566  have evaluated : RationalMatrixOps.entry(RationalRowCertificates.rowAt(threeResult, 0), 3) =567    RationalCell.fraction(31, 30, 1)568    compute569  rw evaluated570  unfold RationalCellOps.meaning571  compute572573theorem threeTrace.quotient_y574  shows RationalCellOps.meaning(RationalMatrixOps.entry(575    RationalRowCertificates.rowAt(threeResult, 1), 3)) = rationalClass(2, 0, 1)576proof577  have evaluated : RationalMatrixOps.entry(RationalRowCertificates.rowAt(threeResult, 1), 3) =578    RationalCell.fraction(16, 14, 1)579    compute580  rw evaluated581  unfold RationalCellOps.meaning582  compute583584theorem threeTrace.quotient_z585  shows RationalCellOps.meaning(RationalMatrixOps.entry(586    RationalRowCertificates.rowAt(threeResult, 2), 3)) = rationalClass(3, 0, 1)587proof588  have evaluated : RationalMatrixOps.entry(RationalRowCertificates.rowAt(threeResult, 2), 3) =589    RationalCell.fraction(14, 11, 1)590    compute591  rw evaluated592  unfold RationalCellOps.meaning593  compute594595-- Equal left sides with unequal constants expose 0 = 1.596def inconsistentRows : Set :=597  List.cons(List.cons(cell(1), List.cons(cell(1), List.cons(cell(2), List.nil))),598    List.cons(List.cons(cell(1), List.cons(cell(1), List.cons(cell(3), List.nil))), List.nil))599def inconsistentTrace : Set :=600  List.cons(RationalRowOperation.addMultiple(1, 0, minusOne), List.nil)601def inconsistentResult : Set :=602  RationalRowCertificates.replay(inconsistentTrace, inconsistentRows)603604theorem inconsistentTrace.checked605  shows RationalRowCertificates.replayChecked(inconsistentTrace, inconsistentRows) = Option.some(inconsistentResult)606proof607  compute608609theorem inconsistentTrace.zeroLeft610  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalRowCertificates.rowAt(inconsistentResult, 1), 0)) = Bool.true ∧611    RationalCellOps.isZero(RationalMatrixOps.entry(RationalRowCertificates.rowAt(inconsistentResult, 1), 1)) = Bool.true612proof613  constructor614  . compute615  . compute616617theorem inconsistentTrace.nonzeroRight618  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalRowCertificates.rowAt(inconsistentResult, 1), 2)) = Bool.false619proof620  compute621622-- A dependent equation reduces to a zero row, leaving one free variable.623def dependentRows : Set :=624  List.cons(List.cons(cell(1), List.cons(cell(1), List.cons(cell(2), List.nil))),625    List.cons(List.cons(cell(2), List.cons(cell(2), List.cons(cell(4), List.nil))), List.nil))626def dependentTrace : Set :=627  List.cons(RationalRowOperation.addMultiple(1, 0, RationalCell.fraction(0, 2, 1)), List.nil)628def dependentResult : Set :=629  RationalRowCertificates.replay(dependentTrace, dependentRows)630631theorem dependentTrace.checked632  shows RationalRowCertificates.replayChecked(dependentTrace, dependentRows) = Option.some(dependentResult)633proof634  compute635636theorem dependentTrace.zeroRow637  shows RationalCellOps.isZero(RationalMatrixOps.entry(RationalRowCertificates.rowAt(dependentResult, 1), 0)) = Bool.true ∧638    RationalCellOps.isZero(RationalMatrixOps.entry(RationalRowCertificates.rowAt(dependentResult, 1), 1)) = Bool.true ∧639    RationalCellOps.isZero(RationalMatrixOps.entry(RationalRowCertificates.rowAt(dependentResult, 1), 2)) = Bool.true640proof641  constructor642  . compute643  . constructor644    . compute645    . compute646647end RationalRowTraceExamples648649namespace RationalRowCertificates650651-- The executable row operations commute with quotient interpretation at an652-- entry. The tails are arbitrary, so these facts apply inside longer rows.653theorem scale.tail (x xs factor : Set)654  shows RationalMatrixOps.drop(RationalMatrixOps.scale(List.cons(x, xs), factor)) =655    RationalMatrixOps.scale(xs, factor)656proof657  simp [RationalMatrixOps.drop, RationalMatrixOps.scale]658659theorem addScaled.tail (x xs source factor : Set)660  shows RationalMatrixOps.drop(RationalMatrixOps.addScaled(List.cons(x, xs), source, factor)) =661    RationalMatrixOps.addScaled(xs, RationalMatrixOps.drop(source), factor)662proof663  simp [RationalMatrixOps.drop, RationalMatrixOps.addScaled]664665theorem scale.first_meaning (a b m c d n : Nat, rest : Set)666  requires factorNonzero : m ≠ 0667  requires entryNonzero : n ≠ 0668  shows RationalCellOps.meaning(RationalMatrixOps.entry(669    RationalMatrixOps.scale(List.cons(RationalCell.fraction(c, d, n), rest),670      RationalCell.fraction(a, b, m)), 0)) =671    rationalMulValue(rationalClass(a, b, m), rationalClass(c, d, n))672proof673  have sound := RationalCellOps.mul.sound at a, b, m, c, d, n with674    factorNonzero, entryNonzero675  simpa [RationalMatrixOps.entry, RationalMatrixOps.head,676    RationalMatrixOps.scale] using sound677678theorem addScaled.first_meaning (a b m c d n e f k : Nat, targetTail sourceTail : Set)679  requires targetNonzero : m ≠ 0680  requires factorNonzero : n ≠ 0681  requires sourceNonzero : k ≠ 0682  shows RationalCellOps.meaning(RationalMatrixOps.entry(683    RationalMatrixOps.addScaled(List.cons(RationalCell.fraction(a, b, m), targetTail),684      List.cons(RationalCell.fraction(e, f, k), sourceTail),685      RationalCell.fraction(c, d, n)), 0)) =686    rationalAddValue(rationalClass(a, b, m),687      rationalMulValue(rationalClass(c, d, n), rationalClass(e, f, k)))688proof689  have product := RationalCellOps.mul.sound at c, d, n, e, f, k with690    factorNonzero, sourceNonzero691  have denominator : n * k ≠ 0 by exact mul.nonzero at n, k with factorNonzero, sourceNonzero692  have sum := RationalCellOps.add.sound at a, b, m,693    c * e + d * f, c * f + d * e, n * k with targetNonzero, denominator694  have intermediate : RationalCellOps.meaning(RationalMatrixOps.entry(695    RationalMatrixOps.addScaled(List.cons(RationalCell.fraction(a, b, m), targetTail),696      List.cons(RationalCell.fraction(e, f, k), sourceTail),697      RationalCell.fraction(c, d, n)), 0)) =698    rationalAddValue(rationalClass(a, b, m),699      RationalCellOps.meaning(RationalCellOps.cellMul(RationalCell.fraction(c, d, n),700        RationalCell.fraction(e, f, k))))701    simpa [RationalMatrixOps.entry, RationalMatrixOps.head,702      RationalMatrixOps.addScaled, RationalCellOps.cellMul,703      RationalCellOps.multiplyBy, RationalCellOps.meaning] using sum704  rwa product at intermediate705706theorem replay.empty (rows : Set)707  shows replay(List.nil, rows) = rows708proof709  simp [replay]710711theorem replay.cons (step rest rows : Set)712  shows replay(List.cons(step, rest), rows) = replay(rest, execute(rows, step))713proof714  simp [replay]715716theorem swapRows.two (a b : Set)717  shows swapRows(List.cons(a, List.cons(b, List.nil)), 0, 1) =718    List.cons(b, List.cons(a, List.nil))719proof720  settle with swapRows, replaceRow, replaceHead, replaceTail, rowAt,721    RationalMatrixOps.firstRow, RationalMatrixOps.restRows722723theorem swapRows.two_involutive (a b : Set)724  shows swapRows(swapRows(List.cons(a, List.cons(b, List.nil)), 0, 1), 0, 1) =725    List.cons(a, List.cons(b, List.nil))726proof727  have first := swapRows.two with a, b728  rw first729  exact swapRows.two with b, a730731-- A replayable witness for the two-equation augmented system above.732-- r₂ ← r₂ − r₁; r₂ ← −r₂/2; r₁ ← r₁ − r₂.733def twoEquationTrace : Set :=734  List.cons(RationalRowOperation.addMultiple(1, 0, RationalCell.fraction(0, 1, 1)),735    List.cons(RationalRowOperation.scale(1, RationalCell.fraction(0, 1, 2)),736      List.cons(RationalRowOperation.addMultiple(0, 1, RationalCell.fraction(0, 1, 1)), List.nil)))737738def twoEquationResult : Set :=739  replay(twoEquationTrace, RationalEliminationExamples.system)740741theorem twoEquationTrace.valid742  shows validTrace(twoEquationTrace) = Bool.true743proof744  compute745746theorem twoEquationTrace.checked747  shows replayChecked(twoEquationTrace, RationalEliminationExamples.system) =748    Option.some(twoEquationResult)749proof750  compute751752theorem invalidScale.rejected753  shows replayChecked(List.cons(RationalRowOperation.scale(0, RationalCellOps.zero), List.nil),754    RationalEliminationExamples.system) = Option.none755proof756  compute757758theorem invalidIndex.rejected759  shows replayChecked(List.cons(RationalRowOperation.swap(0, 2), List.nil),760    RationalEliminationExamples.system) = Option.none761proof762  compute763764theorem malformedFactor.rejected765  shows replayChecked(List.cons(RationalRowOperation.addMultiple(0, 1,766    RationalCell.fraction(1, 0, 0)), List.nil), RationalEliminationExamples.system) = Option.none767proof768  compute769770theorem raggedRows.rejected771  shows replayChecked(List.nil, List.cons(RationalEliminationExamples.rowOne,772    List.cons(List.cons(RationalCellOps.one, List.nil), List.nil))) = Option.none773proof774  compute775776theorem malformedInput.rejected777  shows replayChecked(List.nil,778    List.cons(List.cons(RationalCell.fraction(1, 0, 0), List.nil), List.nil)) = Option.none779proof780  compute781782theorem twoEquationTrace.firstPivot783  shows RationalCellOps.sameValue(RationalMatrixOps.entry(rowAt(twoEquationResult, 0), 0),784    RationalCellOps.one) = Bool.true785proof786  compute787788theorem twoEquationTrace.clearedAbove789  shows RationalCellOps.isZero(RationalMatrixOps.entry(rowAt(twoEquationResult, 0), 1)) = Bool.true790proof791  compute792793theorem twoEquationTrace.clearedBelow794  shows RationalCellOps.isZero(RationalMatrixOps.entry(rowAt(twoEquationResult, 1), 0)) = Bool.true795proof796  compute797798theorem twoEquationTrace.secondPivot799  shows RationalCellOps.sameValue(RationalMatrixOps.entry(rowAt(twoEquationResult, 1), 1),800    RationalCellOps.one) = Bool.true801proof802  compute803804theorem twoEquationTrace.solutionX805  shows RationalCellOps.sameValue(RationalMatrixOps.entry(rowAt(twoEquationResult, 0), 2),806    RationalCell.fraction(2, 0, 1)) = Bool.true807proof808  compute809810theorem twoEquationTrace.solutionY811  shows RationalCellOps.sameValue(RationalMatrixOps.entry(rowAt(twoEquationResult, 1), 2),812    RationalCellOps.one) = Bool.true813proof814  compute815816end RationalRowCertificates817