diff options
Diffstat (limited to 'plugins/fourier')
| -rw-r--r-- | plugins/fourier/Fourier.v | 21 | ||||
| -rw-r--r-- | plugins/fourier/Fourier_util.v | 222 | ||||
| -rw-r--r-- | plugins/fourier/fourier.ml | 205 | ||||
| -rw-r--r-- | plugins/fourier/fourierR.ml | 629 | ||||
| -rw-r--r-- | plugins/fourier/fourier_plugin.mllib | 4 | ||||
| -rw-r--r-- | plugins/fourier/g_fourier.ml4 | 17 | ||||
| -rw-r--r-- | plugins/fourier/vo.itarget | 2 |
7 files changed, 1100 insertions, 0 deletions
diff --git a/plugins/fourier/Fourier.v b/plugins/fourier/Fourier.v new file mode 100644 index 00000000..07b2973a --- /dev/null +++ b/plugins/fourier/Fourier.v @@ -0,0 +1,21 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(************************************************************************) + +(* $Id$ *) + +(* "Fourier's method to solve linear inequations/equations systems.".*) + +Require Export LegacyRing. +Require Export LegacyField. +Require Export DiscrR. +Require Export Fourier_util. +Declare ML Module "fourier_plugin". + +Ltac fourier := abstract (fourierz; field; discrR). + +Ltac fourier_eq := apply Rge_antisym; fourier. diff --git a/plugins/fourier/Fourier_util.v b/plugins/fourier/Fourier_util.v new file mode 100644 index 00000000..0fd92d60 --- /dev/null +++ b/plugins/fourier/Fourier_util.v @@ -0,0 +1,222 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(************************************************************************) + +(* $Id$ *) + +Require Export Rbase. +Comments "Lemmas used by the tactic Fourier". + +Open Scope R_scope. + +Lemma Rfourier_lt : forall x1 y1 a:R, x1 < y1 -> 0 < a -> a * x1 < a * y1. +intros; apply Rmult_lt_compat_l; assumption. +Qed. + +Lemma Rfourier_le : forall x1 y1 a:R, x1 <= y1 -> 0 < a -> a * x1 <= a * y1. +red in |- *. +intros. +case H; auto with real. +Qed. + +Lemma Rfourier_lt_lt : + forall x1 y1 x2 y2 a:R, + x1 < y1 -> x2 < y2 -> 0 < a -> x1 + a * x2 < y1 + a * y2. +intros x1 y1 x2 y2 a H H0 H1; try assumption. +apply Rplus_lt_compat. +try exact H. +apply Rfourier_lt. +try exact H0. +try exact H1. +Qed. + +Lemma Rfourier_lt_le : + forall x1 y1 x2 y2 a:R, + x1 < y1 -> x2 <= y2 -> 0 < a -> x1 + a * x2 < y1 + a * y2. +intros x1 y1 x2 y2 a H H0 H1; try assumption. +case H0; intros. +apply Rplus_lt_compat. +try exact H. +apply Rfourier_lt; auto with real. +rewrite H2. +rewrite (Rplus_comm y1 (a * y2)). +rewrite (Rplus_comm x1 (a * y2)). +apply Rplus_lt_compat_l. +try exact H. +Qed. + +Lemma Rfourier_le_lt : + forall x1 y1 x2 y2 a:R, + x1 <= y1 -> x2 < y2 -> 0 < a -> x1 + a * x2 < y1 + a * y2. +intros x1 y1 x2 y2 a H H0 H1; try assumption. +case H; intros. +apply Rfourier_lt_le; auto with real. +rewrite H2. +apply Rplus_lt_compat_l. +apply Rfourier_lt; auto with real. +Qed. + +Lemma Rfourier_le_le : + forall x1 y1 x2 y2 a:R, + x1 <= y1 -> x2 <= y2 -> 0 < a -> x1 + a * x2 <= y1 + a * y2. +intros x1 y1 x2 y2 a H H0 H1; try assumption. +case H0; intros. +red in |- *. +left; try assumption. +apply Rfourier_le_lt; auto with real. +rewrite H2. +case H; intros. +red in |- *. +left; try assumption. +rewrite (Rplus_comm x1 (a * y2)). +rewrite (Rplus_comm y1 (a * y2)). +apply Rplus_lt_compat_l. +try exact H3. +rewrite H3. +red in |- *. +right; try assumption. +auto with real. +Qed. + +Lemma Rlt_zero_pos_plus1 : forall x:R, 0 < x -> 0 < 1 + x. +intros x H; try assumption. +rewrite Rplus_comm. +apply Rle_lt_0_plus_1. +red in |- *; auto with real. +Qed. + +Lemma Rlt_mult_inv_pos : forall x y:R, 0 < x -> 0 < y -> 0 < x * / y. +intros x y H H0; try assumption. +replace 0 with (x * 0). +apply Rmult_lt_compat_l; auto with real. +ring. +Qed. + +Lemma Rlt_zero_1 : 0 < 1. +exact Rlt_0_1. +Qed. + +Lemma Rle_zero_pos_plus1 : forall x:R, 0 <= x -> 0 <= 1 + x. +intros x H; try assumption. +case H; intros. +red in |- *. +left; try assumption. +apply Rlt_zero_pos_plus1; auto with real. +rewrite <- H0. +replace (1 + 0) with 1. +red in |- *; left. +exact Rlt_zero_1. +ring. +Qed. + +Lemma Rle_mult_inv_pos : forall x y:R, 0 <= x -> 0 < y -> 0 <= x * / y. +intros x y H H0; try assumption. +case H; intros. +red in |- *; left. +apply Rlt_mult_inv_pos; auto with real. +rewrite <- H1. +red in |- *; right; ring. +Qed. + +Lemma Rle_zero_1 : 0 <= 1. +red in |- *; left. +exact Rlt_zero_1. +Qed. + +Lemma Rle_not_lt : forall n d:R, 0 <= n * / d -> ~ 0 < - n * / d. +intros n d H; red in |- *; intros H0; try exact H0. +generalize (Rgt_not_le 0 (n * / d)). +intros H1; elim H1; try assumption. +replace (n * / d) with (- - (n * / d)). +replace 0 with (- -0). +replace (- (n * / d)) with (- n * / d). +replace (-0) with 0. +red in |- *. +apply Ropp_gt_lt_contravar. +red in |- *. +exact H0. +ring. +ring. +ring. +ring. +Qed. + +Lemma Rnot_lt0 : forall x:R, ~ 0 < 0 * x. +intros x; try assumption. +replace (0 * x) with 0. +apply Rlt_irrefl. +ring. +Qed. + +Lemma Rlt_not_le_frac_opp : forall n d:R, 0 < n * / d -> ~ 0 <= - n * / d. +intros n d H; try assumption. +apply Rgt_not_le. +replace 0 with (-0). +replace (- n * / d) with (- (n * / d)). +apply Ropp_lt_gt_contravar. +try exact H. +ring. +ring. +Qed. + +Lemma Rnot_lt_lt : forall x y:R, ~ 0 < y - x -> ~ x < y. +unfold not in |- *; intros. +apply H. +apply Rplus_lt_reg_r with x. +replace (x + 0) with x. +replace (x + (y - x)) with y. +try exact H0. +ring. +ring. +Qed. + +Lemma Rnot_le_le : forall x y:R, ~ 0 <= y - x -> ~ x <= y. +unfold not in |- *; intros. +apply H. +case H0; intros. +left. +apply Rplus_lt_reg_r with x. +replace (x + 0) with x. +replace (x + (y - x)) with y. +try exact H1. +ring. +ring. +right. +rewrite H1; ring. +Qed. + +Lemma Rfourier_gt_to_lt : forall x y:R, y > x -> x < y. +unfold Rgt in |- *; intros; assumption. +Qed. + +Lemma Rfourier_ge_to_le : forall x y:R, y >= x -> x <= y. +intros x y; exact (Rge_le y x). +Qed. + +Lemma Rfourier_eqLR_to_le : forall x y:R, x = y -> x <= y. +exact Req_le. +Qed. + +Lemma Rfourier_eqRL_to_le : forall x y:R, y = x -> x <= y. +exact Req_le_sym. +Qed. + +Lemma Rfourier_not_ge_lt : forall x y:R, (x >= y -> False) -> x < y. +exact Rnot_ge_lt. +Qed. + +Lemma Rfourier_not_gt_le : forall x y:R, (x > y -> False) -> x <= y. +exact Rnot_gt_le. +Qed. + +Lemma Rfourier_not_le_gt : forall x y:R, (x <= y -> False) -> x > y. +exact Rnot_le_lt. +Qed. + +Lemma Rfourier_not_lt_ge : forall x y:R, (x < y -> False) -> x >= y. +exact Rnot_lt_ge. +Qed. diff --git a/plugins/fourier/fourier.ml b/plugins/fourier/fourier.ml new file mode 100644 index 00000000..73fb4929 --- /dev/null +++ b/plugins/fourier/fourier.ml @@ -0,0 +1,205 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(************************************************************************) + +(* $Id$ *) + +(* Méthode d'élimination de Fourier *) +(* Référence: +Auteur(s) : Fourier, Jean-Baptiste-Joseph + +Titre(s) : Oeuvres de Fourier [Document électronique]. Tome second. Mémoires publiés dans divers recueils / publ. par les soins de M. Gaston Darboux,... + +Publication : Numérisation BnF de l'édition de Paris : Gauthier-Villars, 1890 + +Pages: 326-327 + +http://gallica.bnf.fr/ +*) + +(* Un peu de calcul sur les rationnels... +Les opérations rendent des rationnels normalisés, +i.e. le numérateur et le dénominateur sont premiers entre eux. +*) +type rational = {num:int; + den:int} +;; +let print_rational x = + print_int x.num; + print_string "/"; + print_int x.den +;; + +let rec pgcd x y = if y = 0 then x else pgcd y (x mod y);; + + +let r0 = {num=0;den=1};; +let r1 = {num=1;den=1};; + +let rnorm x = let x = (if x.den<0 then {num=(-x.num);den=(-x.den)} else x) in + if x.num=0 then r0 + else (let d=pgcd x.num x.den in + let d= (if d<0 then -d else d) in + {num=(x.num)/d;den=(x.den)/d});; + +let rop x = rnorm {num=(-x.num);den=x.den};; + +let rplus x y = rnorm {num=x.num*y.den + y.num*x.den;den=x.den*y.den};; + +let rminus x y = rnorm {num=x.num*y.den - y.num*x.den;den=x.den*y.den};; + +let rmult x y = rnorm {num=x.num*y.num;den=x.den*y.den};; + +let rinv x = rnorm {num=x.den;den=x.num};; + +let rdiv x y = rnorm {num=x.num*y.den;den=x.den*y.num};; + +let rinf x y = x.num*y.den < y.num*x.den;; +let rinfeq x y = x.num*y.den <= y.num*x.den;; + +(* {coef;hist;strict}, où coef=[c1; ...; cn; d], représente l'inéquation +c1x1+...+cnxn < d si strict=true, <= sinon, +hist donnant les coefficients (positifs) d'une combinaison linéaire qui permet d'obtenir l'inéquation à partir de celles du départ. +*) + +type ineq = {coef:rational list; + hist:rational list; + strict:bool};; + +let pop x l = l:=x::(!l);; + +(* sépare la liste d'inéquations s selon que leur premier coefficient est +négatif, nul ou positif. *) +let partitionne s = + let lpos=ref [] in + let lneg=ref [] in + let lnul=ref [] in + List.iter (fun ie -> match ie.coef with + [] -> raise (Failure "empty ineq") + |(c::r) -> if rinf c r0 + then pop ie lneg + else if rinf r0 c then pop ie lpos + else pop ie lnul) + s; + [!lneg;!lnul;!lpos] +;; +(* initialise les histoires d'une liste d'inéquations données par leurs listes de coefficients et leurs strictitudes (!): +(add_hist [(equation 1, s1);...;(équation n, sn)]) += +[{équation 1, [1;0;...;0], s1}; + {équation 2, [0;1;...;0], s2}; + ... + {équation n, [0;0;...;1], sn}] +*) +let add_hist le = + let n = List.length le in + let i=ref 0 in + List.map (fun (ie,s) -> + let h =ref [] in + for k=1 to (n-(!i)-1) do pop r0 h; done; + pop r1 h; + for k=1 to !i do pop r0 h; done; + i:=!i+1; + {coef=ie;hist=(!h);strict=s}) + le +;; +(* additionne deux inéquations *) +let ie_add ie1 ie2 = {coef=List.map2 rplus ie1.coef ie2.coef; + hist=List.map2 rplus ie1.hist ie2.hist; + strict=ie1.strict || ie2.strict} +;; +(* multiplication d'une inéquation par un rationnel (positif) *) +let ie_emult a ie = {coef=List.map (fun x -> rmult a x) ie.coef; + hist=List.map (fun x -> rmult a x) ie.hist; + strict= ie.strict} +;; +(* on enlève le premier coefficient *) +let ie_tl ie = {coef=List.tl ie.coef;hist=ie.hist;strict=ie.strict} +;; +(* le premier coefficient: "tête" de l'inéquation *) +let hd_coef ie = List.hd ie.coef +;; + +(* calcule toutes les combinaisons entre inéquations de tête négative et inéquations de tête positive qui annulent le premier coefficient. +*) +let deduce_add lneg lpos = + let res=ref [] in + List.iter (fun i1 -> + List.iter (fun i2 -> + let a = rop (hd_coef i1) in + let b = hd_coef i2 in + pop (ie_tl (ie_add (ie_emult b i1) + (ie_emult a i2))) res) + lpos) + lneg; + !res +;; +(* élimination de la première variable à partir d'une liste d'inéquations: +opération qu'on itère dans l'algorithme de Fourier. +*) +let deduce1 s = + match (partitionne s) with + [lneg;lnul;lpos] -> + let lnew = deduce_add lneg lpos in + (List.map ie_tl lnul)@lnew + |_->assert false +;; +(* algorithme de Fourier: on élimine successivement toutes les variables. +*) +let deduce lie = + let n = List.length (fst (List.hd lie)) in + let lie=ref (add_hist lie) in + for i=1 to n-1 do + lie:= deduce1 !lie; + done; + !lie +;; + +(* donne [] si le système a des solutions, +sinon donne [c,s,lc] +où lc est la combinaison linéaire des inéquations de départ +qui donne 0 < c si s=true + ou 0 <= c sinon +cette inéquation étant absurde. +*) +let unsolvable lie = + let lr = deduce lie in + let res = ref [] in + (try (List.iter (fun e -> + match e with + {coef=[c];hist=lc;strict=s} -> + if (rinf c r0 && (not s)) || (rinfeq c r0 && s) + then (res := [c,s,lc]; + raise (Failure "contradiction found")) + |_->assert false) + lr) + with _ -> ()); + !res +;; + +(* Exemples: + +let test1=[[r1;r1;r0],true;[rop r1;r1;r1],false;[r0;rop r1;rop r1],false];; +deduce test1;; +unsolvable test1;; + +let test2=[ +[r1;r1;r0;r0;r0],false; +[r0;r1;r1;r0;r0],false; +[r0;r0;r1;r1;r0],false; +[r0;r0;r0;r1;r1],false; +[r1;r0;r0;r0;r1],false; +[rop r1;rop r1;r0;r0;r0],false; +[r0;rop r1;rop r1;r0;r0],false; +[r0;r0;rop r1;rop r1;r0],false; +[r0;r0;r0;rop r1;rop r1],false; +[rop r1;r0;r0;r0;rop r1],false +];; +deduce test2;; +unsolvable test2;; + +*) diff --git a/plugins/fourier/fourierR.ml b/plugins/fourier/fourierR.ml new file mode 100644 index 00000000..3f490bab --- /dev/null +++ b/plugins/fourier/fourierR.ml @@ -0,0 +1,629 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(************************************************************************) + +(* $Id$ *) + + + +(* La tactique Fourier ne fonctionne de manière sûre que si les coefficients +des inéquations et équations sont entiers. En attendant la tactique Field. +*) + +open Term +open Tactics +open Clenv +open Names +open Libnames +open Tacticals +open Tacmach +open Fourier +open Contradiction + +(****************************************************************************** +Opérations sur les combinaisons linéaires affines. +La partie homogène d'une combinaison linéaire est en fait une table de hash +qui donne le coefficient d'un terme du calcul des constructions, +qui est zéro si le terme n'y est pas. +*) + +type flin = {fhom:(constr , rational)Hashtbl.t; + fcste:rational};; + +let flin_zero () = {fhom=Hashtbl.create 50;fcste=r0};; + +let flin_coef f x = try (Hashtbl.find f.fhom x) with _-> r0;; + +let flin_add f x c = + let cx = flin_coef f x in + Hashtbl.remove f.fhom x; + Hashtbl.add f.fhom x (rplus cx c); + f +;; +let flin_add_cste f c = + {fhom=f.fhom; + fcste=rplus f.fcste c} +;; + +let flin_one () = flin_add_cste (flin_zero()) r1;; + +let flin_plus f1 f2 = + let f3 = flin_zero() in + Hashtbl.iter (fun x c -> let _=flin_add f3 x c in ()) f1.fhom; + Hashtbl.iter (fun x c -> let _=flin_add f3 x c in ()) f2.fhom; + flin_add_cste (flin_add_cste f3 f1.fcste) f2.fcste; +;; + +let flin_minus f1 f2 = + let f3 = flin_zero() in + Hashtbl.iter (fun x c -> let _=flin_add f3 x c in ()) f1.fhom; + Hashtbl.iter (fun x c -> let _=flin_add f3 x (rop c) in ()) f2.fhom; + flin_add_cste (flin_add_cste f3 f1.fcste) (rop f2.fcste); +;; +let flin_emult a f = + let f2 = flin_zero() in + Hashtbl.iter (fun x c -> let _=flin_add f2 x (rmult a c) in ()) f.fhom; + flin_add_cste f2 (rmult a f.fcste); +;; + +(*****************************************************************************) +open Vernacexpr + +type ineq = Rlt | Rle | Rgt | Rge + +let string_of_R_constant kn = + match Names.repr_con kn with + | MPfile dir, sec_dir, id when + sec_dir = empty_dirpath && + string_of_dirpath dir = "Coq.Reals.Rdefinitions" + -> string_of_label id + | _ -> "constant_not_of_R" + +let rec string_of_R_constr c = + match kind_of_term c with + Cast (c,_,_) -> string_of_R_constr c + |Const c -> string_of_R_constant c + | _ -> "not_of_constant" + +let rec rational_of_constr c = + match kind_of_term c with + | Cast (c,_,_) -> (rational_of_constr c) + | App (c,args) -> + (match (string_of_R_constr c) with + | "Ropp" -> + rop (rational_of_constr args.(0)) + | "Rinv" -> + rinv (rational_of_constr args.(0)) + | "Rmult" -> + rmult (rational_of_constr args.(0)) + (rational_of_constr args.(1)) + | "Rdiv" -> + rdiv (rational_of_constr args.(0)) + (rational_of_constr args.(1)) + | "Rplus" -> + rplus (rational_of_constr args.(0)) + (rational_of_constr args.(1)) + | "Rminus" -> + rminus (rational_of_constr args.(0)) + (rational_of_constr args.(1)) + | _ -> failwith "not a rational") + | Const kn -> + (match (string_of_R_constant kn) with + "R1" -> r1 + |"R0" -> r0 + | _ -> failwith "not a rational") + | _ -> failwith "not a rational" +;; + +let rec flin_of_constr c = + try( + match kind_of_term c with + | Cast (c,_,_) -> (flin_of_constr c) + | App (c,args) -> + (match (string_of_R_constr c) with + "Ropp" -> + flin_emult (rop r1) (flin_of_constr args.(0)) + | "Rplus"-> + flin_plus (flin_of_constr args.(0)) + (flin_of_constr args.(1)) + | "Rminus"-> + flin_minus (flin_of_constr args.(0)) + (flin_of_constr args.(1)) + | "Rmult"-> + (try (let a=(rational_of_constr args.(0)) in + try (let b = (rational_of_constr args.(1)) in + (flin_add_cste (flin_zero()) (rmult a b))) + with _-> (flin_add (flin_zero()) + args.(1) + a)) + with _-> (flin_add (flin_zero()) + args.(0) + (rational_of_constr args.(1)))) + | "Rinv"-> + let a=(rational_of_constr args.(0)) in + flin_add_cste (flin_zero()) (rinv a) + | "Rdiv"-> + (let b=(rational_of_constr args.(1)) in + try (let a = (rational_of_constr args.(0)) in + (flin_add_cste (flin_zero()) (rdiv a b))) + with _-> (flin_add (flin_zero()) + args.(0) + (rinv b))) + |_->assert false) + | Const c -> + (match (string_of_R_constant c) with + "R1" -> flin_one () + |"R0" -> flin_zero () + |_-> assert false) + |_-> assert false) + with _ -> flin_add (flin_zero()) + c + r1 +;; + +let flin_to_alist f = + let res=ref [] in + Hashtbl.iter (fun x c -> res:=(c,x)::(!res)) f; + !res +;; + +(* Représentation des hypothèses qui sont des inéquations ou des équations. +*) +type hineq={hname:constr; (* le nom de l'hypothèse *) + htype:string; (* Rlt, Rgt, Rle, Rge, eqTLR ou eqTRL *) + hleft:constr; + hright:constr; + hflin:flin; + hstrict:bool} +;; + +(* Transforme une hypothese h:t en inéquation flin<0 ou flin<=0 +*) +let ineq1_of_constr (h,t) = + match (kind_of_term t) with + App (f,args) -> + (match kind_of_term f with + Const c when Array.length args = 2 -> + let t1= args.(0) in + let t2= args.(1) in + (match (string_of_R_constant c) with + "Rlt" -> [{hname=h; + htype="Rlt"; + hleft=t1; + hright=t2; + hflin= flin_minus (flin_of_constr t1) + (flin_of_constr t2); + hstrict=true}] + |"Rgt" -> [{hname=h; + htype="Rgt"; + hleft=t2; + hright=t1; + hflin= flin_minus (flin_of_constr t2) + (flin_of_constr t1); + hstrict=true}] + |"Rle" -> [{hname=h; + htype="Rle"; + hleft=t1; + hright=t2; + hflin= flin_minus (flin_of_constr t1) + (flin_of_constr t2); + hstrict=false}] + |"Rge" -> [{hname=h; + htype="Rge"; + hleft=t2; + hright=t1; + hflin= flin_minus (flin_of_constr t2) + (flin_of_constr t1); + hstrict=false}] + |_->assert false) + | Ind (kn,i) -> + if IndRef(kn,i) = Coqlib.glob_eq then + let t0= args.(0) in + let t1= args.(1) in + let t2= args.(2) in + (match (kind_of_term t0) with + Const c -> + (match (string_of_R_constant c) with + "R"-> + [{hname=h; + htype="eqTLR"; + hleft=t1; + hright=t2; + hflin= flin_minus (flin_of_constr t1) + (flin_of_constr t2); + hstrict=false}; + {hname=h; + htype="eqTRL"; + hleft=t2; + hright=t1; + hflin= flin_minus (flin_of_constr t2) + (flin_of_constr t1); + hstrict=false}] + |_-> assert false) + |_-> assert false) + else + assert false + |_-> assert false) + |_-> assert false +;; + +(* Applique la méthode de Fourier à une liste d'hypothèses (type hineq) +*) + +let fourier_lineq lineq1 = + let nvar=ref (-1) in + let hvar=Hashtbl.create 50 in (* la table des variables des inéquations *) + List.iter (fun f -> + Hashtbl.iter (fun x _ -> if not (Hashtbl.mem hvar x) then begin + nvar:=(!nvar)+1; + Hashtbl.add hvar x (!nvar) + end) + f.hflin.fhom) + lineq1; + let sys= List.map (fun h-> + let v=Array.create ((!nvar)+1) r0 in + Hashtbl.iter (fun x c -> v.(Hashtbl.find hvar x)<-c) + h.hflin.fhom; + ((Array.to_list v)@[rop h.hflin.fcste],h.hstrict)) + lineq1 in + unsolvable sys +;; + +(*********************************************************************) +(* Defined constants *) + +let get = Lazy.force +let constant = Coqlib.gen_constant "Fourier" + +(* Standard library *) +open Coqlib +let coq_sym_eqT = lazy (build_coq_eq_sym ()) +let coq_False = lazy (build_coq_False ()) +let coq_not = lazy (build_coq_not ()) +let coq_eq = lazy (build_coq_eq ()) + +(* Rdefinitions *) +let constant_real = constant ["Reals";"Rdefinitions"] + +let coq_Rlt = lazy (constant_real "Rlt") +let coq_Rgt = lazy (constant_real "Rgt") +let coq_Rle = lazy (constant_real "Rle") +let coq_Rge = lazy (constant_real "Rge") +let coq_R = lazy (constant_real "R") +let coq_Rminus = lazy (constant_real "Rminus") +let coq_Rmult = lazy (constant_real "Rmult") +let coq_Rplus = lazy (constant_real "Rplus") +let coq_Ropp = lazy (constant_real "Ropp") +let coq_Rinv = lazy (constant_real "Rinv") +let coq_R0 = lazy (constant_real "R0") +let coq_R1 = lazy (constant_real "R1") + +(* RIneq *) +let coq_Rinv_1 = lazy (constant ["Reals";"RIneq"] "Rinv_1") + +(* Fourier_util *) +let constant_fourier = constant ["fourier";"Fourier_util"] + +let coq_Rlt_zero_1 = lazy (constant_fourier "Rlt_zero_1") +let coq_Rlt_zero_pos_plus1 = lazy (constant_fourier "Rlt_zero_pos_plus1") +let coq_Rle_zero_pos_plus1 = lazy (constant_fourier "Rle_zero_pos_plus1") +let coq_Rlt_mult_inv_pos = lazy (constant_fourier "Rlt_mult_inv_pos") +let coq_Rle_zero_zero = lazy (constant_fourier "Rle_zero_zero") +let coq_Rle_zero_1 = lazy (constant_fourier "Rle_zero_1") +let coq_Rle_mult_inv_pos = lazy (constant_fourier "Rle_mult_inv_pos") +let coq_Rnot_lt0 = lazy (constant_fourier "Rnot_lt0") +let coq_Rle_not_lt = lazy (constant_fourier "Rle_not_lt") +let coq_Rfourier_gt_to_lt = lazy (constant_fourier "Rfourier_gt_to_lt") +let coq_Rfourier_ge_to_le = lazy (constant_fourier "Rfourier_ge_to_le") +let coq_Rfourier_eqLR_to_le = lazy (constant_fourier "Rfourier_eqLR_to_le") +let coq_Rfourier_eqRL_to_le = lazy (constant_fourier "Rfourier_eqRL_to_le") + +let coq_Rfourier_not_ge_lt = lazy (constant_fourier "Rfourier_not_ge_lt") +let coq_Rfourier_not_gt_le = lazy (constant_fourier "Rfourier_not_gt_le") +let coq_Rfourier_not_le_gt = lazy (constant_fourier "Rfourier_not_le_gt") +let coq_Rfourier_not_lt_ge = lazy (constant_fourier "Rfourier_not_lt_ge") +let coq_Rfourier_lt = lazy (constant_fourier "Rfourier_lt") +let coq_Rfourier_le = lazy (constant_fourier "Rfourier_le") +let coq_Rfourier_lt_lt = lazy (constant_fourier "Rfourier_lt_lt") +let coq_Rfourier_lt_le = lazy (constant_fourier "Rfourier_lt_le") +let coq_Rfourier_le_lt = lazy (constant_fourier "Rfourier_le_lt") +let coq_Rfourier_le_le = lazy (constant_fourier "Rfourier_le_le") +let coq_Rnot_lt_lt = lazy (constant_fourier "Rnot_lt_lt") +let coq_Rnot_le_le = lazy (constant_fourier "Rnot_le_le") +let coq_Rlt_not_le_frac_opp = lazy (constant_fourier "Rlt_not_le_frac_opp") + +(****************************************************************************** +Construction de la preuve en cas de succès de la méthode de Fourier, +i.e. on obtient une contradiction. +*) +let is_int x = (x.den)=1 +;; + +(* fraction = couple (num,den) *) +let rec rational_to_fraction x= (x.num,x.den) +;; + +(* traduction -3 -> (Ropp (Rplus R1 (Rplus R1 R1))) +*) +let int_to_real n = + let nn=abs n in + if nn=0 + then get coq_R0 + else + (let s=ref (get coq_R1) in + for i=1 to (nn-1) do s:=mkApp (get coq_Rplus,[|get coq_R1;!s|]) done; + if n<0 then mkApp (get coq_Ropp, [|!s|]) else !s) +;; +(* -1/2 -> (Rmult (Ropp R1) (Rinv (Rplus R1 R1))) +*) +let rational_to_real x = + let (n,d)=rational_to_fraction x in + mkApp (get coq_Rmult, + [|int_to_real n;mkApp(get coq_Rinv,[|int_to_real d|])|]) +;; + +(* preuve que 0<n*1/d +*) +let tac_zero_inf_pos gl (n,d) = + let tacn=ref (apply (get coq_Rlt_zero_1)) in + let tacd=ref (apply (get coq_Rlt_zero_1)) in + for i=1 to n-1 do + tacn:=(tclTHEN (apply (get coq_Rlt_zero_pos_plus1)) !tacn); done; + for i=1 to d-1 do + tacd:=(tclTHEN (apply (get coq_Rlt_zero_pos_plus1)) !tacd); done; + (tclTHENS (apply (get coq_Rlt_mult_inv_pos)) [!tacn;!tacd]) +;; + +(* preuve que 0<=n*1/d +*) +let tac_zero_infeq_pos gl (n,d)= + let tacn=ref (if n=0 + then (apply (get coq_Rle_zero_zero)) + else (apply (get coq_Rle_zero_1))) in + let tacd=ref (apply (get coq_Rlt_zero_1)) in + for i=1 to n-1 do + tacn:=(tclTHEN (apply (get coq_Rle_zero_pos_plus1)) !tacn); done; + for i=1 to d-1 do + tacd:=(tclTHEN (apply (get coq_Rlt_zero_pos_plus1)) !tacd); done; + (tclTHENS (apply (get coq_Rle_mult_inv_pos)) [!tacn;!tacd]) +;; + +(* preuve que 0<(-n)*(1/d) => False +*) +let tac_zero_inf_false gl (n,d) = + if n=0 then (apply (get coq_Rnot_lt0)) + else + (tclTHEN (apply (get coq_Rle_not_lt)) + (tac_zero_infeq_pos gl (-n,d))) +;; + +(* preuve que 0<=(-n)*(1/d) => False +*) +let tac_zero_infeq_false gl (n,d) = + (tclTHEN (apply (get coq_Rlt_not_le_frac_opp)) + (tac_zero_inf_pos gl (-n,d))) +;; + +let create_meta () = mkMeta(Evarutil.new_meta());; + +let my_cut c gl= + let concl = pf_concl gl in + apply_type (mkProd(Anonymous,c,concl)) [create_meta()] gl +;; + +let exact = exact_check;; + +let tac_use h = match h.htype with + "Rlt" -> exact h.hname + |"Rle" -> exact h.hname + |"Rgt" -> (tclTHEN (apply (get coq_Rfourier_gt_to_lt)) + (exact h.hname)) + |"Rge" -> (tclTHEN (apply (get coq_Rfourier_ge_to_le)) + (exact h.hname)) + |"eqTLR" -> (tclTHEN (apply (get coq_Rfourier_eqLR_to_le)) + (exact h.hname)) + |"eqTRL" -> (tclTHEN (apply (get coq_Rfourier_eqRL_to_le)) + (exact h.hname)) + |_->assert false +;; + +(* +let is_ineq (h,t) = + match (kind_of_term t) with + App (f,args) -> + (match (string_of_R_constr f) with + "Rlt" -> true + | "Rgt" -> true + | "Rle" -> true + | "Rge" -> true +(* Wrong:not in Rdefinitions: *) | "eqT" -> + (match (string_of_R_constr args.(0)) with + "R" -> true + | _ -> false) + | _ ->false) + |_->false +;; +*) + +let list_of_sign s = List.map (fun (x,_,z)->(x,z)) s;; + +let mkAppL a = + let l = Array.to_list a in + mkApp(List.hd l, Array.of_list (List.tl l)) +;; + +(* Résolution d'inéquations linéaires dans R *) +let rec fourier gl= + Coqlib.check_required_library ["Coq";"fourier";"Fourier"]; + let goal = strip_outer_cast (pf_concl gl) in + let fhyp=id_of_string "new_hyp_for_fourier" in + (* si le but est une inéquation, on introduit son contraire, + et le but à prouver devient False *) + try (let tac = + match (kind_of_term goal) with + App (f,args) -> + (match (string_of_R_constr f) with + "Rlt" -> + (tclTHEN + (tclTHEN (apply (get coq_Rfourier_not_ge_lt)) + (intro_using fhyp)) + fourier) + |"Rle" -> + (tclTHEN + (tclTHEN (apply (get coq_Rfourier_not_gt_le)) + (intro_using fhyp)) + fourier) + |"Rgt" -> + (tclTHEN + (tclTHEN (apply (get coq_Rfourier_not_le_gt)) + (intro_using fhyp)) + fourier) + |"Rge" -> + (tclTHEN + (tclTHEN (apply (get coq_Rfourier_not_lt_ge)) + (intro_using fhyp)) + fourier) + |_->assert false) + |_->assert false + in tac gl) + with _ -> + (* les hypothèses *) + let hyps = List.map (fun (h,t)-> (mkVar h,t)) + (list_of_sign (pf_hyps gl)) in + let lineq =ref [] in + List.iter (fun h -> try (lineq:=(ineq1_of_constr h)@(!lineq)) + with _ -> ()) + hyps; + (* lineq = les inéquations découlant des hypothèses *) + if !lineq=[] then Util.error "No inequalities"; + let res=fourier_lineq (!lineq) in + let tac=ref tclIDTAC in + if res=[] + then Util.error "fourier failed" + (* l'algorithme de Fourier a réussi: on va en tirer une preuve Coq *) + else (match res with + [(cres,sres,lc)]-> + (* lc=coefficients multiplicateurs des inéquations + qui donnent 0<cres ou 0<=cres selon sres *) + (*print_string "Fourier's method can prove the goal...";flush stdout;*) + let lutil=ref [] in + List.iter + (fun (h,c) -> + if c<>r0 + then (lutil:=(h,c)::(!lutil)(*; + print_rational(c);print_string " "*))) + (List.combine (!lineq) lc); + (* on construit la combinaison linéaire des inéquation *) + (match (!lutil) with + (h1,c1)::lutil -> + let s=ref (h1.hstrict) in + let t1=ref (mkAppL [|get coq_Rmult; + rational_to_real c1; + h1.hleft|]) in + let t2=ref (mkAppL [|get coq_Rmult; + rational_to_real c1; + h1.hright|]) in + List.iter (fun (h,c) -> + s:=(!s)||(h.hstrict); + t1:=(mkAppL [|get coq_Rplus; + !t1; + mkAppL [|get coq_Rmult; + rational_to_real c; + h.hleft|] |]); + t2:=(mkAppL [|get coq_Rplus; + !t2; + mkAppL [|get coq_Rmult; + rational_to_real c; + h.hright|] |])) + lutil; + let ineq=mkAppL [|if (!s) then get coq_Rlt else get coq_Rle; + !t1; + !t2 |] in + let tc=rational_to_real cres in + (* puis sa preuve *) + let tac1=ref (if h1.hstrict + then (tclTHENS (apply (get coq_Rfourier_lt)) + [tac_use h1; + tac_zero_inf_pos gl + (rational_to_fraction c1)]) + else (tclTHENS (apply (get coq_Rfourier_le)) + [tac_use h1; + tac_zero_inf_pos gl + (rational_to_fraction c1)])) in + s:=h1.hstrict; + List.iter (fun (h,c)-> + (if (!s) + then (if h.hstrict + then tac1:=(tclTHENS (apply (get coq_Rfourier_lt_lt)) + [!tac1;tac_use h; + tac_zero_inf_pos gl + (rational_to_fraction c)]) + else tac1:=(tclTHENS (apply (get coq_Rfourier_lt_le)) + [!tac1;tac_use h; + tac_zero_inf_pos gl + (rational_to_fraction c)])) + else (if h.hstrict + then tac1:=(tclTHENS (apply (get coq_Rfourier_le_lt)) + [!tac1;tac_use h; + tac_zero_inf_pos gl + (rational_to_fraction c)]) + else tac1:=(tclTHENS (apply (get coq_Rfourier_le_le)) + [!tac1;tac_use h; + tac_zero_inf_pos gl + (rational_to_fraction c)]))); + s:=(!s)||(h.hstrict)) + lutil; + let tac2= if sres + then tac_zero_inf_false gl (rational_to_fraction cres) + else tac_zero_infeq_false gl (rational_to_fraction cres) + in + tac:=(tclTHENS (my_cut ineq) + [tclTHEN (change_in_concl None + (mkAppL [| get coq_not; ineq|] + )) + (tclTHEN (apply (if sres then get coq_Rnot_lt_lt + else get coq_Rnot_le_le)) + (tclTHENS (Equality.replace + (mkAppL [|get coq_Rminus;!t2;!t1|] + ) + tc) + [tac2; + (tclTHENS + (Equality.replace + (mkApp (get coq_Rinv, + [|get coq_R1|])) + (get coq_R1)) +(* en attendant Field, ça peut aider Ring de remplacer 1/1 par 1 ... *) + + [tclORELSE + (Ring.polynom []) + tclIDTAC; + (tclTHEN (apply (get coq_sym_eqT)) + (apply (get coq_Rinv_1)))] + + ) + ])); + !tac1]); + tac:=(tclTHENS (cut (get coq_False)) + [tclTHEN intro (contradiction None); + !tac]) + |_-> assert false) |_-> assert false + ); +(* ((tclTHEN !tac (tclFAIL 1 (* 1 au hasard... *))) gl) *) + (!tac gl) +(* ((tclABSTRACT None !tac) gl) *) + +;; + +(* +let fourier_tac x gl = + fourier gl +;; + +let v_fourier = add_tactic "Fourier" fourier_tac +*) + diff --git a/plugins/fourier/fourier_plugin.mllib b/plugins/fourier/fourier_plugin.mllib new file mode 100644 index 00000000..0383b1a8 --- /dev/null +++ b/plugins/fourier/fourier_plugin.mllib @@ -0,0 +1,4 @@ +Fourier +FourierR +G_fourier +Fourier_plugin_mod diff --git a/plugins/fourier/g_fourier.ml4 b/plugins/fourier/g_fourier.ml4 new file mode 100644 index 00000000..b952851f --- /dev/null +++ b/plugins/fourier/g_fourier.ml4 @@ -0,0 +1,17 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) +(* \VV/ **************************************************************) +(* // * This file is distributed under the terms of the *) +(* * GNU Lesser General Public License Version 2.1 *) +(************************************************************************) + +(*i camlp4deps: "parsing/grammar.cma" i*) + +(* $Id$ *) + +open FourierR + +TACTIC EXTEND fourier + [ "fourierz" ] -> [ fourier ] +END diff --git a/plugins/fourier/vo.itarget b/plugins/fourier/vo.itarget new file mode 100644 index 00000000..87d82dac --- /dev/null +++ b/plugins/fourier/vo.itarget @@ -0,0 +1,2 @@ +Fourier_util.vo +Fourier.vo |