Cramer’s formula: using ring with mathematical components

In an article with David Pichardie that dates back to 2001, we consider determinants to compute the oriented surface of triangles, of the following form

Latex formula

In some place, we need a formula that is called Cramer’s formula, which can be summarized as follows:
Latex formula

In more recent work, I wish to reuse this formula in a proof written in the context of the mathematical components library, where the notions of matrices and determinants are a nice description. It is natural to use the ring tactic for this purpose, but this turns out to be a little difficult. This articles shows how to get through.

As a first step, let’s see how to describe the determinants. We start by choosing a number type
in which to work. In the current case, we later want to work with orders, so that we choose to work
in a type
R that has to be a numDomainType (a ring with an order attached to it).

Section Cramer_formula.

Variable R : realDomainType.

We do this by first writing the matrices
and then the determinant of these matrices. To make the matrices appear in nice table form, we will
write them as sequences of sequences (the table is sequence of rows, and each row is given as a sequence).

Definition seq_to_matrix n m lrows : 'M[R]_(n,m) :=
  \matrix_(i < n, j < m)
     nth 0 (nth [::] lrows i) j.

the determinant
Latex formula
will thus be described by the following text in Coq.

\det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_c; y_c]])

Now, we would like to have this determinant expressed directly as a formula using only multiplications, additions, and subtractions. This is done thanks to the following lemma.

Lemma expand_det_surface x1 y1 x2 y2 x3 y3 :
  \det (seq_to_matrix 3 3 [:: [::1; x1 ; y1]; [::1; x2; y2]; [::1; x3; y3]]) =
  x2 * y3 - y2 * x3 - x1 * y3 + x1 * y2 + y1 * x3 - y1 * x2.
Proof.
by repeat rewrite (expand_det_row _ ord0) /cofactor /row' /col'
 !big_ord_recr /= big_ord0 add0r /= 
 !(mxE, expr1, expr0, expr2, det_mx00, mul1r, mulNr, mulrN, opprK, mulr1,
   mulrDr, opprD, addrA). 
Qed.

I actually expect the tactic that is used here to expand the determinant to work for any determinant of matrices of arbitrary sizes (but the size should not be to big, otherwise it will take a lot of time).

Stating the main lemma for Cramer’s formula

The equality we want to prove is simply expressed as follows, when compared to the above
mathematical formula, we needed to add information about the size of matrices involved. This
size is not guessed from the size of the tables given as argument. If lists that are too
short are given in the table, the missing data is filled with zeros.

Lemma Cramer_formula x_a y_a x_b y_b x_c y_c x_d y_d x_e y_e :
 \det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_d; y_d]]) * \det (seq_to_matrix 3 3
                                        [:: [:: 1; x_a; y_a];
                                            [:: 1; x_c; y_c];
                                            [:: 1; x_e; y_e]]) =
 \det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_e; y_e]]) * \det (seq_to_matrix 3 3
                                        [:: [:: 1; x_a; y_a];
                                            [:: 1; x_c; y_c];
                                            [:: 1; x_d; y_d]]) +
 \det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_c; y_c]]) * \det (seq_to_matrix 3 3
                                        [:: [:: 1; x_a; y_a];
                                            [:: 1; x_d; y_d];
                                            [:: 1; x_e; y_e]]).
Proof.

Using ring

To use the ring tactic, we need a ring structure to be declared for the type handled in our equality.
The first naive approach to declaring this ring structure is to write the following text, but

(* WARNING : FAILING EXAMPLE, DO NOT COPY FOR YOUR OWN USAGE. GOOD EXAMPLE COMES LATER. *) 
Definition R_theory := 
      mk_rt 0 1 +%R *%R (fun x y => x - y) (@GRing.opp _) (@eq R) 
        (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R) (@mulrA R) (@mulrDl R)
        (fun x y : R => erefl (x - y)) (@addrN R).

Add Ring R_Ring : R_theory.

At the time of filling up this description of the ring theory, we discover a small difficulty: the
traditional approach to a ring theory expects the simultaneous existence of a subtraction function and
a unary oposite function. In the mathematical component library, the custom is to define only the unary opposite function and to make the binary subtraction appear simply as a notation on the pattern (x + – y).
For this reason, we need to provide an anonymous function

(fun x y => x - y)

where mk_theory is expecting a name for the binary function, and the proof that binary subtraction and unary opposite are consistent can be a simple invocation to reflexivity of equality:

fun x y : R => erefl (x - y)

This appears as last but one argument to the mk_rt theory builder.

Coq accepts this definition without complaining, so we can proceed with our proof, hoping that everything works well. Here is a candidate script:

rewrite expand_det_surface.
ring.

Alas, after the ring command runs, we only get the following error message.

In nested Ltac calls to "ring" and
"ring_lookup (tactic0) [ (constr_list) ] (ne_constr_list)", last call failed.
In nested Ltac calls to "ring" and
"ring_lookup (tactic0) [ (constr_list) ] (ne_constr_list)", last call failed.
cannot find a declared ring structure over "(GRing.Ring.sort R)"

This message may look quite cryptic, but it contains some valuable information: the type for which
the ring structure should exist is

(GRing.Ring.sort R)

We never wrote such a formula, for us the two members of the equality are in type R, but the Coq
system has a different view. Actually, R was declared as a realDomainType, but the system thinks
that the formula live in some ring structure. Let’s see what is the type that was stored in
our description of the ring theory.

Check R_theory.

The result is:

R_theory
     : ring_theory 0 1 +%R  *%R (fun x y : R => x - y) -%R eq

This tells us that the R_theory is an object of the ring_theory type family. Let’s investigate
ring_theory.

Print ring_theory.
Record ring_theory (R : Type) (rO rI : R) (radd rmul rsub : R -> R -> R)
(ropp : R -> R) (req : R -> R -> Prop) : Prop := mk_rt
  { Radd_0_l : forall x : R, req (radd rO x) x;
    Radd_comm : forall x y : R, req (radd x y) (radd y x);
    Radd_assoc : forall x y z : R,
                 req (radd x (radd y z)) (radd (radd x y) z);
    Rmul_1_l : forall x : R, req (rmul rI x) x;
    Rmul_comm : forall x y : R, req (rmul x y) (rmul y x);
    Rmul_assoc : forall x y z : R,
                 req (rmul x (rmul y z)) (rmul (rmul x y) z);
    Rdistr_l : forall x y z : R,
               req (rmul (radd x y) z) (radd (rmul x z) (rmul y z));
    Rsub_def : forall x y : R, req (rsub x y) (radd x (ropp y));
    Ropp_def : forall x : R, req (radd x (ropp x)) rO }

For ring_theory: Argument R is implicit
For mk_rt: Argument R is implicit
For ring_theory: Argument scopes are [type_scope _ _ function_scope
                   function_scope function_scope function_scope
                   function_scope]
For mk_rt: Argument scopes are [type_scope _ _ function_scope function_scope
             function_scope function_scope function_scope function_scope
             function_scope function_scope function_scope function_scope
             function_scope function_scope function_scope function_scope]

So this tells us that the type argument R was implicit when calling mk_rt. To see the type that
is stored ring’s table, we need to have the implicit argument of ring_theory to be displayed when
displaying the type of R_theory.

Set Printing Implicit.
Set Printing Coercions.
Check R_theory.

The result looks as follows:

R_theory
     : @ring_theory
         (GRing.Zmodule.sort
            (GRing.Ring.zmodType (Num.RealDomain.ringType R))) 0 1 
         +%R  *%R
         (fun
            x
             y : GRing.Zmodule.sort
                   (GRing.Ring.zmodType (Num.RealDomain.ringType R)) => 
          x - y) -%R (@eq (GRing.Ring.sort (Num.RealDomain.ringType R)))

So this explains that the type for which the ring theory has been added as the form GRing.Zmodule.sort …
but the type that appears in the current equality is different. We also see that the equality
for which ring will work will be for a type of the form GRing.Ring.sort.

Definition R_theory :=
  @mk_rt (GRing.Ring.sort R) 0 1 +%R *%R (fun x y => x - y) (@GRing.opp _)
   (@eq (GRing.Ring.sort (Num.RealDomain.ringType R)))
   (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)
     (@mulrA R) (@mulrDl R) (fun x y : R => erefl (x - y)) (@addrN R).

Add Ring R_ring : R_theory.

We can now finish our proof with the ring tactic configured in this fashion.

The whole script is given here:

From mathcomp Require Import all_ssreflect all_algebra.

Set Implicit Arguments.
Unset Strict Implicit.
Unset Printing Implicit Defensive.

Import GRing.Theory Num.Theory.

Open Scope ring_scope.

Section Cramer_formula.

Variable R : realDomainType.

Definition opp (y : R) := - y.

Definition R_theory :=
  @mk_rt (GRing.Ring.sort R) 0 1 +%R *%R (fun x y => x - y) (@GRing.opp _)
   (@eq (GRing.Ring.sort (Num.RealDomain.ringType R)))
   (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)
     (@mulrA R) (@mulrDl R) (fun x y : R => erefl (x - y)) (@addrN R).

Add Ring R_Ring : R_theory.

Definition seq_to_matrix n m lrows : 'M[R]_(n,m) :=
  \matrix_(i < n, j < m)
     nth 0 (nth [::] lrows i) j.

Lemma expand_det_surface x1 y1 x2 y2 x3 y3 :
  \det (seq_to_matrix 3 3 [:: [::1; x1 ; y1]; [::1; x2; y2]; [::1; x3; y3]]) =
  x2 * y3 - y2 * x3 - x1 * y3 + x1 * y2 + y1 * x3 - y1 * x2 :> GRing.Ring.sort R.
Proof.
repeat (rewrite det_mx00 || rewrite (expand_det_row _ ord0) /cofactor /row' /col'
  !big_ord_recr big_ord0 /=).
rewrite !mxE /= !(expr0, exprS).
by rewrite !(mul1r, add0r, mulr1, mulNr, mulrN, opprK, mulrDr, mul0r, mulr0,
         opprD, oppr0, addrA, addr0) .
Qed.

Lemma Cramer_formula x_a y_a x_b y_b x_c y_c x_d y_d x_e y_e :
 \det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_d; y_d]]) * \det (seq_to_matrix 3 3
                                        [:: [:: 1; x_a; y_a];
                                            [:: 1; x_c; y_c];
                                            [:: 1; x_e; y_e]]) =
 \det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_e; y_e]]) * \det (seq_to_matrix 3 3
                                        [:: [:: 1; x_a; y_a];
                                            [:: 1; x_c; y_c];
                                            [:: 1; x_d; y_d]]) +
 \det (seq_to_matrix 3 3
        [:: [:: 1; x_a; y_a];
            [:: 1; x_b; y_b];
            [:: 1; x_c; y_c]]) * \det (seq_to_matrix 3 3
                                        [:: [:: 1; x_a; y_a];
                                            [:: 1; x_d; y_d];
                                            [:: 1; x_e; y_e]]).
Proof.
rewrite !expand_det_surface.
ring.
Qed.

End Cramer_formula.

This script has been tested with Coq 8.7 and mathcomp 1.6.2 installed using opam on February 23rd, 2018.

Unfortunately, there are other contexts where this solution still does not work. In particular,
I wish I could prove Lemma expand_det_surface using ring, too, but for this lemma the proof fails with
a different error being announced.

A more invasive solution

The problem with the proof of expand_det_surface is that a different type comes up as the type in the equality and as the parameter of the ring operations. The type is only different in syntactic form, it is convertible, but the ring tactic does not look this type up modulo conversion. A solution for us is to force a monomorphic shape for all operations before calling ring. This is done in the following excerpt, to be inserted in the Cramer_formula section.

Definition R' := (R : Type).

Definition mul : R' -> R' -> R' := @GRing.mul _.
Definition add : R' -> R' -> R' := @GRing.add _.
Definition sub : R' -> R' -> R' := (fun x y => x - y).
Definition opp : R' -> R' := @GRing.opp _.
Definition zero : R' := 0.
Definition one : R' := 1.

Definition R2_theory :=
  @mk_rt R' zero one add mul sub opp
   (@eq R')
   (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)
     (@mulrA R) (@mulrDl R) (fun x y : R => erefl (x - y)) (@addrN R).

Add Ring R2_Ring : R2_theory.

Ltac mc_ring :=
rewrite ?mxE /= !(expr0, exprS, mulrS, mulr0n) -?[@GRing.add _]/add
   -?[@GRing.mul _]/mul
   -?[@GRing.opp _]/opp -?[1]/one -?[0]/zero;
match goal with |- @eq ?X _ _ => change X with R' end;
ring.

Lemma expand_det_surface' x1 y1 x2 y2 x3 y3 :
  \det (seq_to_matrix 3 3 [:: [::1; x1 ; y1]; [::1; x2; y2]; [::1; x3; y3]]) =
  x2 * y3 - y2 * x3 - x1 * y3 + x1 * y2 + y1 * x3 - y1 * x2.
Proof.
repeat (rewrite det_mx00 || rewrite (expand_det_row _ ord0) /cofactor /row' /col'
  !big_ord_recr big_ord0 /=).
mc_ring.
Qed.

So now, the tactic mc_ring forces all instances of + to be seen as instances of the specific add, and it
also forces the equality we are considering to be viewed as an equality on type R’. As a result, the subsequent call to ring finds the right table to apply and is able to conclude, at least in the example of expand_det_surface, as can be seen here.

Here is another example, where injections of natural numbers into the ring are also present:

Theorem expand_cube (x : R) :
  (x - 1) ^+ 3 = 
  x ^+ 3 - 3%:R * x ^+ 2 + x  + x + x - 1%:R.
Proof.
by mc_ring.
Qed.

One Comment

  1. Pierre-Yves Strub developed a similar layer to integrate ring with mathematical components. Pointers can be found [Here](https://github.com/jasmin-lang/jasmin/blob/master/proofs/3rdparty/ssrring.v)

Leave a Reply

Your email address will not be published.