helm/gTopLevel/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 (* Un peu de calcul sur les rationnels...
  25 Les opérations rendent des rationnels normalisés,
  26 i.e. le numérateur et le dénominateur sont premiers entre eux.
  27 *)
  28 type rational = {num:int;
  29                  den:int}
  30 ;;
  31 let print_rational x =
  32         print_int x.num;
  33         print_string "/";
  34         print_int x.den
  35 ;;
  36
  37 let rec pgcd x y = if y = 0 then x else pgcd y (x mod y);;
  38
  39
  40 let r0 = {num=0;den=1};;
  41 let r1 = {num=1;den=1};;
  42
  43 let rnorm x = let x = (if x.den<0 then {num=(-x.num);den=(-x.den)} else x) in
  44               if x.num=0 then r0
  45               else (let d=pgcd x.num x.den in
  46                     let d= (if d<0 then -d else d) in
  47                     {num=(x.num)/d;den=(x.den)/d});;
  48
  49 let rop x = rnorm {num=(-x.num);den=x.den};;
  50
  51 let rplus x y = rnorm {num=x.num*y.den + y.num*x.den;den=x.den*y.den};;
  52
  53 let rminus x y = rnorm {num=x.num*y.den - y.num*x.den;den=x.den*y.den};;
  54
  55 let rmult x y = rnorm {num=x.num*y.num;den=x.den*y.den};;
  56
  57 let rinv x = rnorm {num=x.den;den=x.num};;
  58
  59 let rdiv x y = rnorm {num=x.num*y.den;den=x.den*y.num};;
  60
  61 let rinf x y = x.num*y.den < y.num*x.den;;
  62 let rinfeq x y = x.num*y.den <= y.num*x.den;;
  63
  64 (* {coef;hist;strict}, où coef=[c1; ...; cn; d], représente l'inéquation
  65 c1x1+...+cnxn < d si strict=true, <= sinon,
  66 hist donnant les coefficients (positifs) d'une combinaison linéaire qui permet d'obtenir l'inéquation à partir de celles du départ.
  67 *)
  68
  69 type ineq = {coef:rational list;
  70              hist:rational list;
  71              strict:bool};;
  72
  73 let pop x l = l:=x::(!l);;
  74
  75 (* sépare la liste d'inéquations s selon que leur premier coefficient est
  76 négatif, nul ou positif. *)
  77 let partitionne s =
  78    let lpos=ref [] in
  79    let lneg=ref [] in
  80    let lnul=ref [] in
  81    List.iter (fun ie -> match ie.coef with
  82                         [] ->  raise (Failure "empty ineq")
  83                        |(c::r) -> if rinf c r0
  84                                   then pop ie lneg
  85                                   else if rinf r0 c then pop ie lpos
  86                                               else pop ie lnul)
  87              s;
  88    [!lneg;!lnul;!lpos]
  89 ;;
  90 (* initialise les histoires d'une liste d'inéquations données par leurs listes de coefficients et leurs strictitudes (!):
  91 (add_hist [(equation 1, s1);...;(équation n, sn)])
  92 =
  93 [{équation 1, [1;0;...;0], s1};
  94  {équation 2, [0;1;...;0], s2};
  95  ...
  96  {équation n, [0;0;...;1], sn}]
  97 *)
  98 let add_hist le =
  99    let n = List.length le in
 100    let i=ref 0 in
 101    List.map (fun (ie,s) ->
 102               let h =ref [] in
 103               for k=1 to (n-(!i)-1) do pop r0 h; done;
 104               pop r1 h;
 105               for k=1 to !i do pop r0 h; done;
 106               i:=!i+1;
 107               {coef=ie;hist=(!h);strict=s})
 108              le
 109 ;;
 110 (* additionne deux inéquations *)
 111 let ie_add ie1 ie2 = {coef=List.map2 rplus ie1.coef ie2.coef;
 112                       hist=List.map2 rplus ie1.hist ie2.hist;
 113                       strict=ie1.strict || ie2.strict}
 114 ;;
 115 (* multiplication d'une inéquation par un rationnel (positif) *)
 116 let ie_emult a ie = {coef=List.map (fun x -> rmult a x) ie.coef;
 117                      hist=List.map (fun x -> rmult a x) ie.hist;
 118                      strict= ie.strict}
 119 ;;
 120 (* on enlève le premier coefficient *)
 121 let ie_tl ie = {coef=List.tl ie.coef;hist=ie.hist;strict=ie.strict}
 122 ;;
 123 (* le premier coefficient: "tête" de l'inéquation *)
 124 let hd_coef ie = List.hd ie.coef
 125 ;;
 126
 127 (* calcule toutes les combinaisons entre inéquations de tête négative et inéquations de tête positive qui annulent le premier coefficient.
 128 *)
 129 let deduce_add lneg lpos =
 130    let res=ref [] in
 131    List.iter (fun i1 ->
 132                  List.iter (fun i2 ->
 133                                 let a = rop (hd_coef i1) in
 134                                 let b = hd_coef i2 in
 135                                 pop (ie_tl (ie_add (ie_emult b i1)
 136                                                    (ie_emult a i2))) res)
 137                            lpos)
 138              lneg;
 139    !res
 140 ;;
 141 (* élimination de la première variable à partir d'une liste d'inéquations:
 142 opération qu'on itère dans l'algorithme de Fourier.
 143 *)
 144 let deduce1 s =
 145     match (partitionne s) with
 146       [lneg;lnul;lpos] ->
 147     let lnew = deduce_add lneg lpos in
 148     (List.map ie_tl lnul)@lnew
 149      |_->assert false
 150 ;;
 151 (* algorithme de Fourier: on élimine successivement toutes les variables.
 152 *)
 153 let deduce lie =
 154    let n = List.length (fst (List.hd lie)) in
 155    let lie=ref (add_hist lie) in
 156    for i=1 to n-1 do
 157       lie:= deduce1 !lie;
 158    done;
 159    !lie
 160 ;;
 161
 162 (* donne [] si le système a des solutions,
 163 sinon donne [c,s,lc]
 164 où lc est la combinaison linéaire des inéquations de départ
 165 qui donne 0 < c si s=true
 166        ou 0 <= c sinon
 167 cette inéquation étant absurde.
 168 *)
 169 let unsolvable lie =
 170    let lr = deduce lie in
 171    let res = ref [] in
 172    (try (List.iter (fun e ->
 173          match e with
 174            {coef=[c];hist=lc;strict=s} ->
 175                   if (rinf c r0 && (not s)) || (rinfeq c r0 && s)
 176                   then (res := [c,s,lc];
 177                         raise (Failure "contradiction found"))
 178           |_->assert false)
 179                              lr)
 180    with _ -> ());
 181    !res
 182 ;;
 183
 184 (* Exemples:
 185
 186 let test1=[[r1;r1;r0],true;[rop r1;r1;r1],false;[r0;rop r1;rop r1],false];;
 187 deduce test1;;
 188 unsolvable test1;;
 189
 190 let test2=[
 191 [r1;r1;r0;r0;r0],false;
 192 [r0;r1;r1;r0;r0],false;
 193 [r0;r0;r1;r1;r0],false;
 194 [r0;r0;r0;r1;r1],false;
 195 [r1;r0;r0;r0;r1],false;
 196 [rop r1;rop r1;r0;r0;r0],false;
 197 [r0;rop r1;rop r1;r0;r0],false;
 198 [r0;r0;rop r1;rop r1;r0],false;
 199 [r0;r0;r0;rop r1;rop r1],false;
 200 [rop r1;r0;r0;r0;rop r1],false
 201 ];;
 202 deduce test2;;
 203 unsolvable test2;;
 204
 205 *)