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/
27 (* Un peu de calcul sur les rationnels...
28 Les opérations rendent des rationnels normalisés,
29 i.e. le numérateur et le dénominateur sont premiers entre eux.
34 type rational = {num:int;
37 let print_rational x =
43 let rec pgcd x y = if y = 0 then x else pgcd y (x mod y);;
46 let r0 = {num=0;den=1};;
47 let r1 = {num=1;den=1};;
49 let rnorm x = let x = (if x.den<0 then {num=(-x.num);den=(-x.den)} else x) in
51 else (let d=pgcd x.num x.den in
52 let d= (if d<0 then -d else d) in
53 {num=(x.num)/d;den=(x.den)/d});;
55 let rop x = rnorm {num=(-x.num);den=x.den};;
57 let rplus x y = rnorm {num=x.num*y.den + y.num*x.den;den=x.den*y.den};;
59 let rminus x y = rnorm {num=x.num*y.den - y.num*x.den;den=x.den*y.den};;
61 let rmult x y = rnorm {num=x.num*y.num;den=x.den*y.den};;
63 let rinv x = rnorm {num=x.den;den=x.num};;
65 let rdiv x y = rnorm {num=x.num*y.den;den=x.den*y.num};;
67 let rinf x y = x.num*y.den < y.num*x.den;;
68 let rinfeq x y = x.num*y.den <= y.num*x.den;;
70 (* {coef;hist;strict}, où coef=[c1; ...; cn; d], représente l'inéquation
71 c1x1+...+cnxn < d si strict=true, <= sinon,
72 hist donnant les coefficients (positifs) d'une combinaison linéaire qui permet d'obtenir l'inéquation à partir de celles du départ.
75 type ineq = {coef:rational list;
79 let pop x l = l:=x::(!l);;
81 (* sépare la liste d'inéquations s selon que leur premier coefficient est
82 négatif, nul ou positif. *)
87 List.iter (fun ie -> match ie.coef with
88 [] -> raise (Failure "empty ineq")
89 |(c::r) -> if rinf c r0
91 else if rinf r0 c then pop ie lpos
96 (* initialise les histoires d'une liste d'inéquations données par leurs listes de coefficients et leurs strictitudes (!):
97 (add_hist [(equation 1, s1);...;(équation n, sn)])
99 [{équation 1, [1;0;...;0], s1};
100 {équation 2, [0;1;...;0], s2};
102 {équation n, [0;0;...;1], sn}]
105 let n = List.length le in
107 List.map (fun (ie,s) ->
109 for k=1 to (n-(!i)-1) do pop r0 h; done;
111 for k=1 to !i do pop r0 h; done;
113 {coef=ie;hist=(!h);strict=s})
116 (* additionne deux inéquations *)
117 let ie_add ie1 ie2 = {coef=List.map2 rplus ie1.coef ie2.coef;
118 hist=List.map2 rplus ie1.hist ie2.hist;
119 strict=ie1.strict || ie2.strict}
121 (* multiplication d'une inéquation par un rationnel (positif) *)
122 let ie_emult a ie = {coef=List.map (fun x -> rmult a x) ie.coef;
123 hist=List.map (fun x -> rmult a x) ie.hist;
126 (* on enlève le premier coefficient *)
127 let ie_tl ie = {coef=List.tl ie.coef;hist=ie.hist;strict=ie.strict}
129 (* le premier coefficient: "tête" de l'inéquation *)
130 let hd_coef ie = List.hd ie.coef
133 (* calcule toutes les combinaisons entre inéquations de tête négative et inéquations de tête positive qui annulent le premier coefficient.
135 let deduce_add lneg lpos =
139 let a = rop (hd_coef i1) in
140 let b = hd_coef i2 in
141 pop (ie_tl (ie_add (ie_emult b i1)
142 (ie_emult a i2))) res)
147 (* élimination de la première variable à partir d'une liste d'inéquations:
148 opération qu'on itère dans l'algorithme de Fourier.
151 match (partitionne s) with
153 let lnew = deduce_add lneg lpos in
154 (List.map ie_tl lnul)@lnew
157 (* algorithme de Fourier: on élimine successivement toutes les variables.
160 let n = List.length (fst (List.hd lie)) in
161 let lie=ref (add_hist lie) in
168 (* donne [] si le système a des solutions,
170 où lc est la combinaison linéaire des inéquations de départ
171 qui donne 0 < c si s=true
173 cette inéquation étant absurde.
176 let lr = deduce lie in
178 (try (List.iter (fun e ->
180 {coef=[c];hist=lc;strict=s} ->
181 if (rinf c r0 && (not s)) || (rinfeq c r0 && s)
182 then (res := [c,s,lc];
183 raise (Failure "contradiction found"))
192 let test1=[[r1;r1;r0],true;[rop r1;r1;r1],false;[r0;rop r1;rop r1],false];;
197 [r1;r1;r0;r0;r0],false;
198 [r0;r1;r1;r0;r0],false;
199 [r0;r0;r1;r1;r0],false;
200 [r0;r0;r0;r1;r1],false;
201 [r1;r0;r0;r0;r1],false;
202 [rop r1;rop r1;r0;r0;r0],false;
203 [r0;rop r1;rop r1;r0;r0],false;
204 [r0;r0;rop r1;rop r1;r0],false;
205 [r0;r0;r0;rop r1;rop r1],false;
206 [rop r1;r0;r0;r0;rop r1],false