1 (***********************************************************************)
2 (* v * The Coq Proof Assistant / The Coq Development Team *)
3 (* <O___,, * INRIA-Rocquencourt & LRI-CNRS-Orsay *)
4 (* \VV/ *************************************************************)
5 (* // * This file is distributed under the terms of the *)
6 (* * GNU Lesser General Public License Version 2.1 *)
7 (***********************************************************************)
11 (* Méthode d'élimination de Fourier *)
13 Auteur(s) : Fourier, Jean-Baptiste-Joseph
15 Titre(s) : Oeuvres de Fourier [Document électronique]. Tome second. Mémoires publiés dans divers recueils / publ. par les soins de M. Gaston Darboux,...
17 Publication : Numérisation BnF de l'édition de Paris : Gauthier-Villars, 1890
21 http://gallica.bnf.fr/
24 (* Un peu de calcul sur les rationnels...
25 Les opérations rendent des rationnels normalisés,
26 i.e. le numérateur et le dénominateur sont premiers entre eux.
28 type rational = {num:int;
31 let print_rational x =
37 let rec pgcd x y = if y = 0 then x else pgcd y (x mod y);;
40 let r0 = {num=0;den=1};;
41 let r1 = {num=1;den=1};;
43 let rnorm x = let x = (if x.den<0 then {num=(-x.num);den=(-x.den)} else x) in
45 else (let d=pgcd x.num x.den in
46 let d= (if d<0 then -d else d) in
47 {num=(x.num)/d;den=(x.den)/d});;
49 let rop x = rnorm {num=(-x.num);den=x.den};;
51 let rplus x y = rnorm {num=x.num*y.den + y.num*x.den;den=x.den*y.den};;
53 let rminus x y = rnorm {num=x.num*y.den - y.num*x.den;den=x.den*y.den};;
55 let rmult x y = rnorm {num=x.num*y.num;den=x.den*y.den};;
57 let rinv x = rnorm {num=x.den;den=x.num};;
59 let rdiv x y = rnorm {num=x.num*y.den;den=x.den*y.num};;
61 let rinf x y = x.num*y.den < y.num*x.den;;
62 let rinfeq x y = x.num*y.den <= y.num*x.den;;
64 (* {coef;hist;strict}, où coef=[c1; ...; cn; d], représente l'inéquation
65 c1x1+...+cnxn < d si strict=true, <= sinon,
66 hist donnant les coefficients (positifs) d'une combinaison linéaire qui permet d'obtenir l'inéquation à partir de celles du départ.
69 type ineq = {coef:rational list;
73 let pop x l = l:=x::(!l);;
75 (* sépare la liste d'inéquations s selon que leur premier coefficient est
76 négatif, nul ou positif. *)
81 List.iter (fun ie -> match ie.coef with
82 [] -> raise (Failure "empty ineq")
83 |(c::r) -> if rinf c r0
85 else if rinf r0 c then pop ie lpos
90 (* initialise les histoires d'une liste d'inéquations données par leurs listes de coefficients et leurs strictitudes (!):
91 (add_hist [(equation 1, s1);...;(équation n, sn)])
93 [{équation 1, [1;0;...;0], s1};
94 {équation 2, [0;1;...;0], s2};
96 {équation n, [0;0;...;1], sn}]
99 let n = List.length le in
101 List.map (fun (ie,s) ->
103 for k=1 to (n-(!i)-1) do pop r0 h; done;
105 for k=1 to !i do pop r0 h; done;
107 {coef=ie;hist=(!h);strict=s})
110 (* additionne deux inéquations *)
111 let ie_add ie1 ie2 = {coef=List.map2 rplus ie1.coef ie2.coef;
112 hist=List.map2 rplus ie1.hist ie2.hist;
113 strict=ie1.strict || ie2.strict}
115 (* multiplication d'une inéquation par un rationnel (positif) *)
116 let ie_emult a ie = {coef=List.map (fun x -> rmult a x) ie.coef;
117 hist=List.map (fun x -> rmult a x) ie.hist;
120 (* on enlève le premier coefficient *)
121 let ie_tl ie = {coef=List.tl ie.coef;hist=ie.hist;strict=ie.strict}
123 (* le premier coefficient: "tête" de l'inéquation *)
124 let hd_coef ie = List.hd ie.coef
127 (* calcule toutes les combinaisons entre inéquations de tête négative et inéquations de tête positive qui annulent le premier coefficient.
129 let deduce_add lneg lpos =
133 let a = rop (hd_coef i1) in
134 let b = hd_coef i2 in
135 pop (ie_tl (ie_add (ie_emult b i1)
136 (ie_emult a i2))) res)
141 (* élimination de la première variable à partir d'une liste d'inéquations:
142 opération qu'on itère dans l'algorithme de Fourier.
145 match (partitionne s) with
147 let lnew = deduce_add lneg lpos in
148 (List.map ie_tl lnul)@lnew
151 (* algorithme de Fourier: on élimine successivement toutes les variables.
154 let n = List.length (fst (List.hd lie)) in
155 let lie=ref (add_hist lie) in
162 (* donne [] si le système a des solutions,
164 où lc est la combinaison linéaire des inéquations de départ
165 qui donne 0 < c si s=true
167 cette inéquation étant absurde.
170 let lr = deduce lie in
172 (try (List.iter (fun e ->
174 {coef=[c];hist=lc;strict=s} ->
175 if (rinf c r0 && (not s)) || (rinfeq c r0 && s)
176 then (res := [c,s,lc];
177 raise (Failure "contradiction found"))
186 let test1=[[r1;r1;r0],true;[rop r1;r1;r1],false;[r0;rop r1;rop r1],false];;
191 [r1;r1;r0;r0;r0],false;
192 [r0;r1;r1;r0;r0],false;
193 [r0;r0;r1;r1;r0],false;
194 [r0;r0;r0;r1;r1],false;
195 [r1;r0;r0;r0;r1],false;
196 [rop r1;rop r1;r0;r0;r0],false;
197 [r0;rop r1;rop r1;r0;r0],false;
198 [r0;r0;rop r1;rop r1;r0],false;
199 [r0;r0;r0;rop r1;rop r1],false;
200 [rop r1;r0;r0;r0;rop r1],false