1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
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: fourier.ml 5920 2004-07-16 20:01:26Z herbelin $ *)
(* 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;;
*)
|
