{"id":188,"date":"2018-02-23T17:22:34","date_gmt":"2018-02-23T16:22:34","guid":{"rendered":"http:\/\/project.inria.fr\/coqexchange\/?p=188"},"modified":"2019-05-07T16:51:37","modified_gmt":"2019-05-07T14:51:37","slug":"cramers-formula-using-ring-with-mathematical-components","status":"publish","type":"post","link":"https:\/\/project.inria.fr\/coqexchange\/cramers-formula-using-ring-with-mathematical-components\/","title":{"rendered":"Cramer&#8217;s formula: using ring with mathematical components"},"content":{"rendered":"<p>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<\/p>\n\\(<br \/>\n\\left|\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_b&amp;y_b\\\\<br \/>\n1&amp;x_c&amp;y_c<br \/>\n\\end{array}\\right|<br \/>\n\\)\n<p>In some place, we need a formula that is called Cramer&#8217;s formula, which can be summarized as follows:<br \/>\n\\(<br \/>\n\\left|\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_b&amp;y_b\\\\<br \/>\n1&amp;x_d&amp;y_d<br \/>\n\\end{array}\\right|<br \/>\n\\left|\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_c&amp;y_c\\\\<br \/>\n1&amp;x_e&amp;y_e<br \/>\n\\end{array}\\right|<br \/>\n=<br \/>\n\\left|<br \/>\n\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_b&amp;y_b\\\\<br \/>\n1&amp;x_c&amp;y_c<br \/>\n\\end{array}\\right|<br \/>\n\\left|\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_d&amp;y_d\\\\<br \/>\n1&amp;x_e&amp;y_e<br \/>\n\\end{array}\\right|+<br \/>\n\\left|<br \/>\n\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_b&amp;y_b\\\\<br \/>\n1&amp;x_e&amp;y_e<br \/>\n\\end{array}\\right|<br \/>\n\\left|\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_c&amp;y_c\\\\<br \/>\n1&amp;x_e&amp;y_e<br \/>\n\\end{array}\\right|<br \/>\n\\)<\/p>\n<p>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.<\/p>\n<p>As a first step, let&#8217;s see how to describe the determinants.  We start by choosing a number type<br \/>\nin which to work.  In the current case, we later want to work with orders, so that we choose to work<br \/>\n in a type<br \/>\nR that has to be a numDomainType (a ring with an order attached to it).<\/p>\n<pre>\r\nSection Cramer_formula.\r\n\r\nVariable R : realDomainType.\r\n<\/pre>\n<p>We do this by first writing the matrices<br \/>\nand then the determinant of these matrices. To make the matrices appear in nice table form, we will<br \/>\nwrite them as sequences of sequences (the table is sequence of rows, and each row is given as a sequence).<\/p>\n<pre>Definition seq_to_matrix n m lrows : 'M[R]_(n,m) :=\r\n  \\matrix_(i &lt; n, j &lt; m)\r\n     nth 0 (nth [::] lrows i) j.\r\n<\/pre>\n<p>the determinant<br \/>\n\\(<br \/>\n\\left|\\begin{array}{ccc}1&amp;x_a &amp;y_a\\\\<br \/>\n1&amp;x_b&amp;y_b\\\\<br \/>\n1&amp;x_c&amp;y_c<br \/>\n\\end{array}\\right|<br \/>\n\\)<br \/>\nwill thus be described by the following text in Coq.<\/p>\n<pre>\\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_c; y_c]])\r\n<\/pre>\n<p>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.<\/p>\n<pre>Lemma expand_det_surface x1 y1 x2 y2 x3 y3 :\r\n  \\det (seq_to_matrix 3 3 [:: [::1; x1 ; y1]; [::1; x2; y2]; [::1; x3; y3]]) =\r\n  x2 * y3 - y2 * x3 - x1 * y3 + x1 * y2 + y1 * x3 - y1 * x2.\r\nProof.\r\nby repeat rewrite (expand_det_row _ ord0) \/cofactor \/row' \/col'\r\n !big_ord_recr \/= big_ord0 add0r \/= \r\n !(mxE, expr1, expr0, expr2, det_mx00, mul1r, mulNr, mulrN, opprK, mulr1,\r\n   mulrDr, opprD, addrA). \r\nQed.\r\n<\/pre>\n<p>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).<\/p>\n<h2>Stating the main lemma for Cramer&#8217;s formula<\/h2>\n<p>The equality we want to prove is simply expressed as follows, when compared to the above<br \/>\nmathematical formula, we needed to add information about the size of matrices involved. This<br \/>\nsize is not guessed from the size of the tables given as argument. If lists that are too<br \/>\nshort are given in the table, the missing data is filled with zeros.<\/p>\n<pre>Lemma Cramer_formula x_a y_a x_b y_b x_c y_c x_d y_d x_e y_e :\r\n \\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_d; y_d]]) * \\det (seq_to_matrix 3 3\r\n                                        [:: [:: 1; x_a; y_a];\r\n                                            [:: 1; x_c; y_c];\r\n                                            [:: 1; x_e; y_e]]) =\r\n \\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_e; y_e]]) * \\det (seq_to_matrix 3 3\r\n                                        [:: [:: 1; x_a; y_a];\r\n                                            [:: 1; x_c; y_c];\r\n                                            [:: 1; x_d; y_d]]) +\r\n \\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_c; y_c]]) * \\det (seq_to_matrix 3 3\r\n                                        [:: [:: 1; x_a; y_a];\r\n                                            [:: 1; x_d; y_d];\r\n                                            [:: 1; x_e; y_e]]).\r\nProof.\r\n<\/pre>\n<h2>Using ring<\/h2>\n<p>To use the ring tactic, we need a ring structure to be declared for the type handled in our equality.<br \/>\nThe first naive approach to declaring this ring structure is to write the following text, but<br \/>\n<div class=\"alert alert-warning\" role=\"alert\"><p class=\"printonly\"><strong>Warning!<\/strong><\/p>beware that this does not work.<\/div><\/p>\n<pre>\r\n(* WARNING : FAILING EXAMPLE, DO NOT COPY FOR YOUR OWN USAGE. GOOD EXAMPLE COMES LATER. *) \r\nDefinition R_theory := \r\n      mk_rt 0 1 +%R *%R (fun x y =&gt; x - y) (@GRing.opp _) (@eq R) \r\n        (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R) (@mulrA R) (@mulrDl R)\r\n        (fun x y : R =&gt; erefl (x - y)) (@addrN R).\r\n\r\nAdd Ring R_Ring : R_theory.\r\n<\/pre>\n<p>At the time of filling up this description of the ring theory, we discover a small difficulty: the<br \/>\ntraditional approach to a ring theory expects the simultaneous existence of a subtraction function and<br \/>\na 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 + &#8211; y).<br \/>\nFor this reason, we need to provide an anonymous function<\/p>\n<pre>(fun x y =&gt; x - y)\r\n<\/pre>\n<p>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:<\/p>\n<pre>\r\nfun x y : R =&gt; erefl (x - y)\r\n<\/pre>\n<p>This appears as last but one argument to the mk_rt theory builder.<\/p>\n<p>Coq accepts this definition without complaining, so we can proceed with our proof, hoping that everything works well. Here is a candidate script:<\/p>\n<pre>rewrite expand_det_surface.\r\nring.\r\n<\/pre>\n<p>Alas, after the ring command runs, we only get the following error message.<\/p>\n<pre>In nested Ltac calls to \"ring\" and\r\n\"ring_lookup (tactic0) [ (constr_list) ] (ne_constr_list)\", last call failed.\r\nIn nested Ltac calls to \"ring\" and\r\n\"ring_lookup (tactic0) [ (constr_list) ] (ne_constr_list)\", last call failed.\r\ncannot find a declared ring structure over \"(GRing.Ring.sort R)\"\r\n<\/pre>\n<p>This message may look quite cryptic, but it contains some valuable information: the type for which<br \/>\nthe ring structure should exist is<\/p>\n<pre>(GRing.Ring.sort R)\r\n<\/pre>\n<p>We never wrote such a formula, for us the two members of the equality are in type R, but the Coq<br \/>\nsystem has a different view.  Actually, R was declared as a realDomainType, but the system thinks<br \/>\nthat the formula live in some ring structure.  Let&#8217;s see what is the type that was stored in<br \/>\nour description of the ring theory.<\/p>\n<pre>\r\nCheck R_theory.\r\n<\/pre>\n<p>The result is:<\/p>\n<pre>\r\nR_theory\r\n     : ring_theory 0 1 +%R  *%R (fun x y : R => x - y) -%R eq\r\n<\/pre>\n<p>This tells us that the R_theory is an object of the ring_theory type family.  Let&#8217;s investigate<br \/>\nring_theory.<\/p>\n<pre>\r\nPrint ring_theory.\r\n<\/pre>\n<pre>\r\nRecord ring_theory (R : Type) (rO rI : R) (radd rmul rsub : R -> R -> R)\r\n(ropp : R -> R) (req : R -> R -> Prop) : Prop := mk_rt\r\n  { Radd_0_l : forall x : R, req (radd rO x) x;\r\n    Radd_comm : forall x y : R, req (radd x y) (radd y x);\r\n    Radd_assoc : forall x y z : R,\r\n                 req (radd x (radd y z)) (radd (radd x y) z);\r\n    Rmul_1_l : forall x : R, req (rmul rI x) x;\r\n    Rmul_comm : forall x y : R, req (rmul x y) (rmul y x);\r\n    Rmul_assoc : forall x y z : R,\r\n                 req (rmul x (rmul y z)) (rmul (rmul x y) z);\r\n    Rdistr_l : forall x y z : R,\r\n               req (rmul (radd x y) z) (radd (rmul x z) (rmul y z));\r\n    Rsub_def : forall x y : R, req (rsub x y) (radd x (ropp y));\r\n    Ropp_def : forall x : R, req (radd x (ropp x)) rO }\r\n\r\nFor ring_theory: Argument R is implicit\r\nFor mk_rt: Argument R is implicit\r\nFor ring_theory: Argument scopes are [type_scope _ _ function_scope\r\n                   function_scope function_scope function_scope\r\n                   function_scope]\r\nFor mk_rt: Argument scopes are [type_scope _ _ function_scope function_scope\r\n             function_scope function_scope function_scope function_scope\r\n             function_scope function_scope function_scope function_scope\r\n             function_scope function_scope function_scope function_scope]\r\n<\/pre>\n<p>So this tells us that the type argument R was implicit when calling mk_rt.  To see the type that<br \/>\nis stored ring&#8217;s table, we need to have the implicit argument of ring_theory to be displayed when<br \/>\ndisplaying the type of R_theory.<\/p>\n<pre>\r\nSet Printing Implicit.\r\nSet Printing Coercions.\r\nCheck R_theory.\r\n<\/pre>\n<p>The result looks as follows:<\/p>\n<pre>\r\nR_theory\r\n     : @ring_theory\r\n         (GRing.Zmodule.sort\r\n            (GRing.Ring.zmodType (Num.RealDomain.ringType R))) 0 1 \r\n         +%R  *%R\r\n         (fun\r\n            x\r\n             y : GRing.Zmodule.sort\r\n                   (GRing.Ring.zmodType (Num.RealDomain.ringType R)) => \r\n          x - y) -%R (@eq (GRing.Ring.sort (Num.RealDomain.ringType R)))\r\n<\/pre>\n<p>So this explains that the type for which the ring theory has been added as the form GRing.Zmodule.sort &#8230;<br \/>\nbut the type that appears in the current equality is different.  We also see that the equality<br \/>\nfor which ring will work will be for a type of the form GRing.Ring.sort.<\/p>\n<pre>\r\nDefinition R_theory :=\r\n  @mk_rt (GRing.Ring.sort R) 0 1 +%R *%R (fun x y => x - y) (@GRing.opp _)\r\n   (@eq (GRing.Ring.sort (Num.RealDomain.ringType R)))\r\n   (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)\r\n     (@mulrA R) (@mulrDl R) (fun x y : R => erefl (x - y)) (@addrN R).\r\n\r\nAdd Ring R_ring : R_theory.\r\n<\/pre>\n<p>We can now finish our proof with the ring tactic configured in this fashion.<\/p>\n<p>The whole script is given here:<\/p>\n<pre>\r\nFrom mathcomp Require Import all_ssreflect all_algebra.\r\n\r\nSet Implicit Arguments.\r\nUnset Strict Implicit.\r\nUnset Printing Implicit Defensive.\r\n\r\nImport GRing.Theory Num.Theory.\r\n\r\nOpen Scope ring_scope.\r\n\r\nSection Cramer_formula.\r\n\r\nVariable R : realDomainType.\r\n\r\nDefinition opp (y : R) := - y.\r\n\r\nDefinition R_theory :=\r\n  @mk_rt (GRing.Ring.sort R) 0 1 +%R *%R (fun x y => x - y) (@GRing.opp _)\r\n   (@eq (GRing.Ring.sort (Num.RealDomain.ringType R)))\r\n   (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)\r\n     (@mulrA R) (@mulrDl R) (fun x y : R => erefl (x - y)) (@addrN R).\r\n\r\nAdd Ring R_Ring : R_theory.\r\n\r\nDefinition seq_to_matrix n m lrows : 'M[R]_(n,m) :=\r\n  \\matrix_(i < n, j < m)\r\n     nth 0 (nth [::] lrows i) j.\r\n\r\nLemma expand_det_surface x1 y1 x2 y2 x3 y3 :\r\n  \\det (seq_to_matrix 3 3 [:: [::1; x1 ; y1]; [::1; x2; y2]; [::1; x3; y3]]) =\r\n  x2 * y3 - y2 * x3 - x1 * y3 + x1 * y2 + y1 * x3 - y1 * x2 :> GRing.Ring.sort R.\r\nProof.\r\nrepeat (rewrite det_mx00 || rewrite (expand_det_row _ ord0) \/cofactor \/row' \/col'\r\n  !big_ord_recr big_ord0 \/=).\r\nrewrite !mxE \/= !(expr0, exprS).\r\nby rewrite !(mul1r, add0r, mulr1, mulNr, mulrN, opprK, mulrDr, mul0r, mulr0,\r\n         opprD, oppr0, addrA, addr0) .\r\nQed.\r\n\r\nLemma Cramer_formula x_a y_a x_b y_b x_c y_c x_d y_d x_e y_e :\r\n \\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_d; y_d]]) * \\det (seq_to_matrix 3 3\r\n                                        [:: [:: 1; x_a; y_a];\r\n                                            [:: 1; x_c; y_c];\r\n                                            [:: 1; x_e; y_e]]) =\r\n \\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_e; y_e]]) * \\det (seq_to_matrix 3 3\r\n                                        [:: [:: 1; x_a; y_a];\r\n                                            [:: 1; x_c; y_c];\r\n                                            [:: 1; x_d; y_d]]) +\r\n \\det (seq_to_matrix 3 3\r\n        [:: [:: 1; x_a; y_a];\r\n            [:: 1; x_b; y_b];\r\n            [:: 1; x_c; y_c]]) * \\det (seq_to_matrix 3 3\r\n                                        [:: [:: 1; x_a; y_a];\r\n                                            [:: 1; x_d; y_d];\r\n                                            [:: 1; x_e; y_e]]).\r\nProof.\r\nrewrite !expand_det_surface.\r\nring.\r\nQed.\r\n\r\nEnd Cramer_formula.\r\n<\/pre>\n<p>This script has been tested with Coq 8.7 and mathcomp 1.6.2 installed using opam on February 23rd, 2018.<\/p>\n<p>Unfortunately, there are other contexts where this solution still does not work.  In particular,<br \/>\nI wish I could prove Lemma expand_det_surface using ring, too, but for this lemma the proof fails with<br \/>\na different error being announced.<\/p>\n<h2>A more invasive solution<\/h2>\n<p>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.<\/p>\n<pre>\r\nDefinition R' := (R : Type).\r\n\r\nDefinition mul : R' -> R' -> R' := @GRing.mul _.\r\nDefinition add : R' -> R' -> R' := @GRing.add _.\r\nDefinition sub : R' -> R' -> R' := (fun x y => x - y).\r\nDefinition opp : R' -> R' := @GRing.opp _.\r\nDefinition zero : R' := 0.\r\nDefinition one : R' := 1.\r\n\r\nDefinition R2_theory :=\r\n  @mk_rt R' zero one add mul sub opp\r\n   (@eq R')\r\n   (@add0r R) (@addrC R) (@addrA R) (@mul1r R) (@mulrC R)\r\n     (@mulrA R) (@mulrDl R) (fun x y : R => erefl (x - y)) (@addrN R).\r\n\r\nAdd Ring R2_Ring : R2_theory.\r\n\r\nLtac mc_ring :=\r\nrewrite ?mxE \/= !(expr0, exprS, mulrS, mulr0n) -?[@GRing.add _]\/add\r\n   -?[@GRing.mul _]\/mul\r\n   -?[@GRing.opp _]\/opp -?[1]\/one -?[0]\/zero;\r\nmatch goal with |- @eq ?X _ _ => change X with R' end;\r\nring.\r\n\r\nLemma expand_det_surface' x1 y1 x2 y2 x3 y3 :\r\n  \\det (seq_to_matrix 3 3 [:: [::1; x1 ; y1]; [::1; x2; y2]; [::1; x3; y3]]) =\r\n  x2 * y3 - y2 * x3 - x1 * y3 + x1 * y2 + y1 * x3 - y1 * x2.\r\nProof.\r\nrepeat (rewrite det_mx00 || rewrite (expand_det_row _ ord0) \/cofactor \/row' \/col'\r\n  !big_ord_recr big_ord0 \/=).\r\nmc_ring.\r\nQed.\r\n<\/pre>\n<p>So now, the tactic mc_ring forces all instances of + to be seen as instances of the specific add, and it<br \/>\nalso forces the equality we are considering to be viewed as an equality on type R&#8217;.  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.<\/p>\n<p>Here is another example, where injections of natural numbers into the ring are also present:<\/p>\n<pre>\r\nTheorem expand_cube (x : R) :\r\n  (x - 1) ^+ 3 = \r\n  x ^+ 3 - 3%:R * x ^+ 2 + x  + x + x - 1%:R.\r\nProof.\r\nby mc_ring.\r\nQed.\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>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 In some place, we need a formula that is called Cramer&#8217;s formula, which can be summarized as follows: In more recent work, I wish to\u2026<\/p>\n<p> <a class=\"continue-reading-link\" href=\"https:\/\/project.inria.fr\/coqexchange\/cramers-formula-using-ring-with-mathematical-components\/\"><span>Continue reading<\/span><i class=\"crycon-right-dir\"><\/i><\/a> <\/p>\n","protected":false},"author":1013,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-188","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/posts\/188","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/users\/1013"}],"replies":[{"embeddable":true,"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/comments?post=188"}],"version-history":[{"count":19,"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/posts\/188\/revisions"}],"predecessor-version":[{"id":236,"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/posts\/188\/revisions\/236"}],"wp:attachment":[{"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/media?parent=188"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/categories?post=188"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/project.inria.fr\/coqexchange\/wp-json\/wp\/v2\/tags?post=188"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}