(***********************************************************************)
(*  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;;

*)