]> matita.cs.unibo.it Git - helm.git/blob - helm/ocaml/tactics/fourier.ml
ready for 0.1.1 release
[helm.git] / helm / ocaml / tactics / fourier.ml
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 (***********************************************************************)
8
9 (* $Id$ *)
10
11 (* Méthode d'élimination de Fourier *)
12 (* Référence:
13 Auteur(s) : Fourier, Jean-Baptiste-Joseph
14  
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,...
16  
17 Publication : Numérisation BnF de l'édition de Paris : Gauthier-Villars, 1890
18  
19 Pages: 326-327
20
21 http://gallica.bnf.fr/
22 *)
23
24 (** @author The Coq Development Team *)
25
26
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.
30 *)
31
32
33 (** Type for coefficents *)
34 type rational = {
35 num:int; (** Numerator *)
36 den:int;  (** Denumerator *)
37 };;
38
39 (** Debug function.
40     @param x the rational to print*)
41 let print_rational x =
42         print_int x.num;
43         print_string "/";
44         print_int x.den
45 ;;
46
47 let rec pgcd x y = if y = 0 then x else pgcd y (x mod y);;
48
49 (** The constant 0*)
50 let r0 = {num=0;den=1};;
51 (** The constant 1*)
52 let r1 = {num=1;den=1};;
53
54 let rnorm x = let x = (if x.den<0 then {num=(-x.num);den=(-x.den)} else x) in
55               if x.num=0 then r0
56               else (let d=pgcd x.num x.den in
57                     let d= (if d<0 then -d else d) in
58                     {num=(x.num)/d;den=(x.den)/d});;
59
60 (** Calculates the opposite of a rational.
61     @param x The rational
62     @return -x*)
63 let rop x = rnorm {num=(-x.num);den=x.den};;
64
65 (** Sums two rationals.
66     @param x A rational
67     @param y Another rational
68     @return x+y*)
69 let rplus x y = rnorm {num=x.num*y.den + y.num*x.den;den=x.den*y.den};;
70 (** Substracts two rationals.
71     @param x A rational
72     @param y Another rational
73     @return x-y*)
74 let rminus x y = rnorm {num=x.num*y.den - y.num*x.den;den=x.den*y.den};;
75 (** Multiplyes two rationals.
76     @param x A rational
77     @param y Another rational
78     @return x*y*)
79 let rmult x y = rnorm {num=x.num*y.num;den=x.den*y.den};;
80 (** Inverts arational.
81     @param x A rational
82     @return x{^ -1} *)
83 let rinv x = rnorm {num=x.den;den=x.num};;
84 (** Divides two rationals.
85     @param x A rational
86     @param y Another rational
87     @return x/y*)
88 let rdiv x y = rnorm {num=x.num*y.den;den=x.den*y.num};;
89
90 let rinf x y = x.num*y.den < y.num*x.den;;
91 let rinfeq x y = x.num*y.den <= y.num*x.den;;
92
93
94 (* {coef;hist;strict}, où coef=[c1; ...; cn; d], représente l'inéquation
95 c1x1+...+cnxn < d si strict=true, <= sinon,
96 hist donnant les coefficients (positifs) d'une combinaison linéaire qui permet d'obtenir l'inéquation à partir de celles du départ.
97 *)
98
99 type ineq = {coef:rational list;
100              hist:rational list;
101              strict:bool};;
102
103 let pop x l = l:=x::(!l);;
104
105 (* sépare la liste d'inéquations s selon que leur premier coefficient est 
106 négatif, nul ou positif. *)
107 let partitionne s =
108    let lpos=ref [] in
109    let lneg=ref [] in
110    let lnul=ref [] in
111    List.iter (fun ie -> match ie.coef with
112                         [] ->  raise (Failure "empty ineq")
113                        |(c::r) -> if rinf c r0
114                                   then pop ie lneg
115                                   else if rinf r0 c then pop ie lpos
116                                               else pop ie lnul)
117              s;
118    [!lneg;!lnul;!lpos]
119 ;;
120 (* initialise les histoires d'une liste d'inéquations données par leurs listes de coefficients et leurs strictitudes (!):
121 (add_hist [(equation 1, s1);...;(équation n, sn)])
122 =
123 [{équation 1, [1;0;...;0], s1};
124  {équation 2, [0;1;...;0], s2};
125  ...
126  {équation n, [0;0;...;1], sn}]
127 *)
128 let add_hist le =
129    let n = List.length le in
130    let i=ref 0 in
131    List.map (fun (ie,s) -> 
132               let h =ref [] in
133               for k=1 to (n-(!i)-1) do pop r0 h; done;
134               pop r1 h;
135               for k=1 to !i do pop r0 h; done;
136               i:=!i+1;
137               {coef=ie;hist=(!h);strict=s})
138              le
139 ;;
140 (* additionne deux inéquations *)      
141 let ie_add ie1 ie2 = {coef=List.map2 rplus ie1.coef ie2.coef;
142                       hist=List.map2 rplus ie1.hist ie2.hist;
143                       strict=ie1.strict || ie2.strict}
144 ;;
145 (* multiplication d'une inéquation par un rationnel (positif) *)
146 let ie_emult a ie = {coef=List.map (fun x -> rmult a x) ie.coef;
147                      hist=List.map (fun x -> rmult a x) ie.hist;
148                      strict= ie.strict}
149 ;;
150 (* on enlève le premier coefficient *)
151 let ie_tl ie = {coef=List.tl ie.coef;hist=ie.hist;strict=ie.strict}
152 ;;
153 (* le premier coefficient: "tête" de l'inéquation *)
154 let hd_coef ie = List.hd ie.coef
155 ;;
156
157 (* calcule toutes les combinaisons entre inéquations de tête négative et inéquations de tête positive qui annulent le premier coefficient.
158 *)
159 let deduce_add lneg lpos =
160    let res=ref [] in
161    List.iter (fun i1 ->
162                  List.iter (fun i2 ->
163                                 let a = rop (hd_coef i1) in
164                                 let b = hd_coef i2 in
165                                 pop (ie_tl (ie_add (ie_emult b i1)
166                                                    (ie_emult a i2))) res)
167                            lpos)
168              lneg;
169    !res
170 ;;
171 (* élimination de la première variable à partir d'une liste d'inéquations:
172 opération qu'on itère dans l'algorithme de Fourier.
173 *)
174 let deduce1 s i=
175     match (partitionne s) with 
176       [lneg;lnul;lpos] ->
177          let lnew = deduce_add lneg lpos in
178          (match lneg with [] -> print_string("non posso ridurre "^string_of_int i^"\n")|_->();
179           match lpos with [] -> print_string("non posso ridurre "^string_of_int i^"\n")|_->());
180          (List.map ie_tl lnul)@lnew
181    |_->assert false
182 ;;
183 (* algorithme de Fourier: on élimine successivement toutes les variables.
184 *)
185 let deduce lie =
186    let n = List.length (fst (List.hd lie)) in
187    let lie=ref (add_hist lie) in
188    for i=1 to n-1 do
189       lie:= deduce1 !lie i;
190    done;
191    !lie
192 ;;
193
194 (* donne [] si le système a des  find solutions,
195 sinon donne [c,s,lc]
196 où lc est la combinaison linéaire des inéquations de départ
197 qui donne 0 < c si s=true
198        ou 0 <= c sinon
199 cette inéquation étant absurde.
200 *)
201 (** Tryes to find if the system admits solutions.
202     @param lie the list of inequations 
203     @return a list that can be empty if the system has solutions. Otherwise it returns a
204             one elements list [\[(c,s,lc)\]]. {b c} is the rational that can be obtained solving the system,
205             {b s} is true if the inequation that proves that the system is absurd is of type [c < 0], false if 
206             [c <= 0], {b lc} is a list of rational that represents the liear combination to obtain the
207             absurd inequation *)
208 let unsolvable lie =
209    let lr = deduce lie in
210    let res = ref [] in
211    (try (List.iter (fun e ->
212          match e with
213            {coef=[c];hist=lc;strict=s} ->
214                   if (rinf c r0 && (not s)) || (rinfeq c r0 && s) 
215                   then (res := [c,s,lc];
216                         raise (Failure "contradiction found"))
217           |_->assert false)
218                              lr)
219    with _ -> ());
220    !res
221 ;;
222
223 (* Exemples:
224
225 let test1=[[r1;r1;r0],true;[rop r1;r1;r1],false;[r0;rop r1;rop r1],false];;
226 deduce test1;;
227 unsolvable test1;;
228
229 let test2=[
230 [r1;r1;r0;r0;r0],false;
231 [r0;r1;r1;r0;r0],false;
232 [r0;r0;r1;r1;r0],false;
233 [r0;r0;r0;r1;r1],false;
234 [r1;r0;r0;r0;r1],false;
235 [rop r1;rop r1;r0;r0;r0],false;
236 [r0;rop r1;rop r1;r0;r0],false;
237 [r0;r0;rop r1;rop r1;r0],false;
238 [r0;r0;r0;rop r1;rop r1],false;
239 [rop r1;r0;r0;r0;rop r1],false
240 ];;
241 deduce test2;;
242 unsolvable test2;;
243
244 *)