--- /dev/null
+(***********************************************************************)
+(* v * The Coq Proof Assistant / The Coq Development Team *)
+(* <O___,, * INRIA-Rocquencourt & LRI-CNRS-Orsay *)
+(* \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/
+*)
+
+(** @author The Coq Development Team *)
+
+
+(* 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 for coefficents *)
+type rational = {
+num:int; (** Numerator *)
+den:int; (** Denumerator *)
+};;
+
+(** Debug function.
+ @param x the rational to print*)
+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);;
+
+(** The constant 0*)
+let r0 = {num=0;den=1};;
+(** The constant 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});;
+
+(** Calculates the opposite of a rational.
+ @param x The rational
+ @return -x*)
+let rop x = rnorm {num=(-x.num);den=x.den};;
+
+(** Sums two rationals.
+ @param x A rational
+ @param y Another rational
+ @return x+y*)
+let rplus x y = rnorm {num=x.num*y.den + y.num*x.den;den=x.den*y.den};;
+(** Substracts two rationals.
+ @param x A rational
+ @param y Another rational
+ @return x-y*)
+let rminus x y = rnorm {num=x.num*y.den - y.num*x.den;den=x.den*y.den};;
+(** Multiplyes two rationals.
+ @param x A rational
+ @param y Another rational
+ @return x*y*)
+let rmult x y = rnorm {num=x.num*y.num;den=x.den*y.den};;
+(** Inverts arational.
+ @param x A rational
+ @return x{^ -1} *)
+let rinv x = rnorm {num=x.den;den=x.num};;
+(** Divides two rationals.
+ @param x A rational
+ @param y Another rational
+ @return x/y*)
+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 i=
+ match (partitionne s) with
+ [lneg;lnul;lpos] ->
+ let lnew = deduce_add lneg lpos in
+ (match lneg with [] -> print_string("non posso ridurre "^string_of_int i^"\n")|_->();
+ match lpos with [] -> print_string("non posso ridurre "^string_of_int i^"\n")|_->());
+ (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 i;
+ done;
+ !lie
+;;
+
+(* donne [] si le système a des find 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.
+*)
+(** Tryes to find if the system admits solutions.
+ @param lie the list of inequations
+ @return a list that can be empty if the system has solutions. Otherwise it returns a
+ one elements list [\[(c,s,lc)\]]. {b c} is the rational that can be obtained solving the system,
+ {b s} is true if the inequation that proves that the system is absurd is of type [c < 0], false if
+ [c <= 0], {b lc} is a list of rational that represents the liear combination to obtain the
+ absurd inequation *)
+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;;
+
+*)
projects / helm.git / blobdiff