micromega : bug fixes and optimisations

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@11318 85f007b7-540e-0410-9357-904b9bb8a0f7

7 files changed, 1450 insertions, 1468 deletions

diff --git a/Makefile.common b/Makefile.common
index af380dd63..6b7da1a57 100644
--- a/Makefile.common
+++ b/Makefile.common
@@ -236,7 +236,7 @@ ROMEGACMO:=\
   contrib/romega/refl_omega.cmo contrib/romega/g_romega.cmo
 
 MICROMEGACMO:=\
-  contrib/micromega/mutils.cmo contrib/micromega/vector.cmo \
+  contrib/micromega/mutils.cmo  \
   contrib/micromega/micromega.cmo contrib/micromega/mfourier.cmo \
   contrib/micromega/certificate.cmo \
   contrib/micromega/coq_micromega.cmo contrib/micromega/g_micromega.cmo
@@ -390,8 +390,7 @@ PARSERCMX:= $(PARSERREQUIRESCMX) $(PARSERCODE:.cmo=.cmx)
 INTERFACERC:= contrib/interface/vernacrc
 
 CSDPCERTCMO:= contrib/micromega/mutils.cmo contrib/micromega/micromega.cmo \
-              contrib/micromega/vector.cmo contrib/micromega/mfourier.cmo \
-              contrib/micromega/certificate.cmo \
+	      contrib/micromega/mfourier.cmo  contrib/micromega/certificate.cmo \
               contrib/micromega/sos.cmo contrib/micromega/csdpcert.cmo
 CSDPCERTCMX:= $(CSDPCERTCMO:.cmo=.cmx)
 
diff --git a/contrib/micromega/certificate.ml b/contrib/micromega/certificate.ml
index f4efcd083..a297c1633 100644
--- a/contrib/micromega/certificate.ml
+++ b/contrib/micromega/certificate.ml
@@ -76,7 +76,7 @@ struct
 end
 
 
-module  Poly :
+module Poly :
  (* A polynomial is a map of monomials *)
  (* 
     This is probably a naive implementation 
@@ -96,6 +96,7 @@ sig
  val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a
  val pp : out_channel -> t -> unit
  val compare : t -> t -> int
+ val is_null : t -> bool
 end =
 struct
  (*normalisation bug : 0*x ... *)
@@ -156,6 +157,8 @@ struct
 
  let fold = P.fold
 
+ let is_null p = fold (fun mn vl b -> b & sign_num vl = 0) p  true
+
  let compare = compare compare_num
 end
 
@@ -332,10 +335,10 @@ let factorise_linear_cone c =
 open Mfourier
    (*module Fourier = Fourier(Vector.VList)(SysSet(Vector.VList))*)
    (*module Fourier = Fourier(Vector.VSparse)(SysSetAlt(Vector.VSparse))*)
-module Fourier = Mfourier.Fourier(Vector.VSparse)(*(SysSetAlt(Vector.VMap))*)
+(*module Fourier = Mfourier.Fourier(Vector.VSparse)(*(SysSetAlt(Vector.VMap))*)*)
 
-module Vect = Fourier.Vect
-open Fourier.Cstr
+(*module Vect = Fourier.Vect*)
+(*open Fourier.Cstr*)
 
 (* fold_left followed by a rev ! *)
 
@@ -388,8 +391,8 @@ let build_linear_system l =
  let strict = {
   coeffs = Vect.from_list ((Big_int unit_big_int)::
                             (List.map (fun (x,y) -> 
-                             match y with Mc.Strict -> 
-                              Big_int unit_big_int 
+			      match y with Mc.Strict -> 
+				  Big_int unit_big_int 
                               | _ -> Big_int zero_big_int) l));
   op = Ge ; cst = Big_int unit_big_int } in
   (* Add the positivity constraint *)
@@ -433,13 +436,41 @@ let make_certificate n_spec cert li =
 
 
 exception Found of Monomial.t
- 
-let raw_certificate l = 
+
+exception Strict
+
+let primal l = 
+  let vr = ref 0 in
+  let module Mmn = Map.Make(Monomial) in
+
+  let vect_of_poly map p =
+    Poly.fold (fun mn vl (map,vect) -> 
+      if mn = Monomial.const 
+      then (map,vect)
+      else 
+	let (mn,m) = try (Mmn.find mn map,map) with Not_found -> let res = (!vr, Mmn.add mn !vr map) in incr vr ; res in
+	  (m,if sign_num vl = 0 then vect else (mn,vl)::vect)) p (map,[]) in
+    
+  let op_op = function Mc.NonStrict -> Ge |Mc.Equal -> Eq | _ -> raise Strict in
+
+  let cmp x y = Pervasives.compare (fst x) (fst y) in
+
+   snd (List.fold_right (fun  (p,op) (map,l) -> 
+      let (mp,vect) = vect_of_poly map p in  
+      let cstr = {coeffs = List.sort cmp vect; op = op_op op ; cst = minus_num (Poly.get Monomial.const p)} in
+
+	(mp,cstr::l)) l (Mmn.empty,[]))
+
+let dual_raw_certificate (l:  (Poly.t * Mc.op1) list) = 
+(*  List.iter (fun (p,op) -> Printf.fprintf stdout "%a %s 0\n" Poly.pp p (string_of_op op) ) l ; *)
+    
+    
  let sys = build_linear_system l in
-  try 
+
+   try 
    match Fourier.find_point sys with
-    | None -> None
-    | Some cert ->  Some (rats_to_ints (Vect.to_list cert)) 
+    | Inr _ -> None
+    | Inl cert ->  Some (rats_to_ints (Vect.to_list cert)) 
        (* should not use rats_to_ints *)
   with x -> 
    if debug 
@@ -448,6 +479,21 @@ let raw_certificate l =
    None
 
 
+let raw_certificate l = 
+  try 
+    let p = primal l in
+      match Fourier.find_point p with
+	| Inr prf -> 
+	    if debug then Printf.printf "AProof : %a\n" pp_proof prf ; 
+	    let cert = List.map (fun (x,n) -> x+1,n) (fst (List.hd (Proof.mk_proof p prf))) in
+	      if debug then Printf.printf "CProof : %a" Vect.pp_vect cert ; 
+	      Some (rats_to_ints (Vect.to_list cert))
+	| Inl _   -> None
+  with Strict -> 
+    (* Fourier elimination should handle > *)
+    dual_raw_certificate l 
+
+
 let simple_linear_prover to_constant l =
  let (lc,li) = List.split l in
   match raw_certificate lc with
@@ -495,7 +541,6 @@ let make_linear_system l =
     ,monomials)
 
 
-open Interval 
 let pplus x y = Mc.PEadd(x,y)
 let pmult x y = Mc.PEmul(x,y)
 let pconst x = Mc.PEc x
@@ -560,11 +605,10 @@ let rec xzlinear_prover planes sys =
 
      let smallest_interval = 
       match List.fold_left (fun (x1,i1) (x2,i2) -> 
-       if smaller_itv i1 i2 
-       then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates 
+       if Itv.smaller_itv i1 i2 
+       then (x1,i1) else (x2,i2)) (Vect.null,(None,None)) candidates 
       with
-       | (x,Itv(Some i, Some j))  -> Some(i,x,j)
-       | (x,Point n) ->  Some(n,x,n)
+       | (x,(Some i, Some j))  -> Some(i,x,j)
        |   x        ->   None (* This might be a cutting plane *)
       in
       match smallest_interval with
@@ -613,7 +657,9 @@ and zlinear_enum  planes expr clb cub l =
 let zlinear_prover sys =  
  let candidates = candidates sys in
  (*  Printf.printf "candidates %d" (List.length candidates) ; *) 
- xzlinear_prover candidates sys
+ (*let t0 = Sys.time () in*)
+ let res = xzlinear_prover candidates sys in
+   (*Printf.printf "Time prover : %f" (Sys.time () -. t0) ;*) res
 
 open Sos
 
diff --git a/contrib/micromega/coq_micromega.ml b/contrib/micromega/coq_micromega.ml
index 02ff61a19..847ae7baa 100644
--- a/contrib/micromega/coq_micromega.ml
+++ b/contrib/micromega/coq_micromega.ml
@@ -23,43 +23,68 @@ let time str f x =
 		     flush stdout); 
   res
 
-type ('a,'b) formula =
+
+type tag = int
+type 'cst atom = 'cst Micromega.formula
+
+type 'cst formula =
   | TT 
   | FF 
-  | X of 'b 
-  | A of 'a * Names.name
-  | C of ('a,'b) formula * ('a,'b) formula * Names.name
-  | D of ('a,'b) formula * ('a,'b) formula * Names.name
-  | N of ('a,'b) formula * Names.name
-  | I of ('a,'b) formula * ('a,'b) formula * Names.name
+  | X of Term.constr
+  | A of 'cst atom * tag * Term.constr
+  | C of 'cst  formula * 'cst formula 
+  | D of 'cst formula * 'cst formula 
+  | N of 'cst formula 
+  | I of 'cst formula * Names.identifier option * 'cst formula 
+
+
+let rec pp_formula o f = 
+  match f with
+    | TT -> output_string  o "tt"
+    | FF -> output_string  o "ff"
+    | X c -> output_string o "X " 
+    | A(_,t,_) -> Printf.fprintf o "A(%i)" t
+    | C(f1,f2) -> Printf.fprintf o "C(%a,%a)" pp_formula f1 pp_formula f2
+    | D(f1,f2) -> Printf.fprintf o "D(%a,%a)" pp_formula f1 pp_formula f2
+    | I(f1,n,f2) -> Printf.fprintf o "I(%a%s,%a)" 
+	pp_formula f1 
+	  (match n with 
+	    | Some id -> Names.string_of_id id 
+	    | None -> "") pp_formula f2
+    | N(f) -> Printf.fprintf o "N(%a)" pp_formula f 
+
+let rec ids_of_formula f =
+  match f with
+    | I(f1,Some id,f2) -> id::(ids_of_formula f2)
+    | _                -> []
+
+(* obsolete *)
+
+let tag_formula t f = f (* to be removed *)
+(* obsolete *)
+
+module ISet = Set.Make(struct type t = int let compare : int -> int -> int = Pervasives.compare end)
 
-let none = Names.Anonymous
+type 'cst clause = ('cst Micromega.nFormula * tag) list 
 
-let tag_formula t f = 
- match f with
-  | A(x,_) -> A(x,t)
-  | C(x,y,_) -> C(x,y,t)
-  | D(x,y,_) -> D(x,y,t)
-  | N(x,_) -> N(x,t)
-  | I(x,y,_) -> I(x,y,t)
-  |   _      -> f
+type 'cst cnf = ('cst clause) list
 
-let tt = []
-let ff = [ [] ]
+let tt : 'cst cnf = []
+let ff : 'cst cnf = [ [] ]
 
 
-type ('constant,'contr) sentence =  
-  ('constant Micromega.formula, 'contr) formula
+type 'cst mc_cnf = ('cst Micromega.nFormula)  Micromega.list Micromega.list
 
-let cnf negate normalise f = 
- let negate a = 
-  CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> x) cl) (negate a) in
+let cnf (negate: 'cst atom -> 'cst mc_cnf) (normalise:'cst atom -> 'cst mc_cnf) (f:'cst formula) = 
+ let negate a t = 
+  CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> (x,t)) cl) (negate a) in
 
- let normalise a = 
-  CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> x) cl) (normalise a) in
+ let normalise a t = 
+     CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> (x,t)) cl) (normalise a) in
 
  let and_cnf x y = x @ y in
- let or_clause_cnf t f =  List.map (fun x -> t@x ) f in
+
+ let or_clause_cnf t f =  List.map (fun x -> t@x) f in
   
  let rec or_cnf f f'  =
   match f with
@@ -71,19 +96,17 @@ let cnf negate normalise f =
    | TT -> if pol then tt else ff (* ?? *)
    | FF  -> if pol then ff else tt (* ?? *)
    | X p -> if pol then ff else ff (* ?? *)
-   | A(x,t) -> if pol then normalise  x else negate x
-   | N(e,t)  -> xcnf  (not pol) e
-   | C(e1,e2,t) -> 
+   | A(x,t,_) -> if pol then normalise  x t else negate x t
+   | N(e)  -> xcnf  (not pol) e
+   | C(e1,e2) -> 
       (if pol then and_cnf else or_cnf) (xcnf  pol e1) (xcnf  pol e2)
-   | D(e1,e2,t)  -> 
+   | D(e1,e2)  -> 
       (if pol then or_cnf else and_cnf) (xcnf  pol e1) (xcnf  pol e2)
-   | I(e1,e2,t) -> 
+   | I(e1,_,e2) -> 
       (if pol then or_cnf else and_cnf) (xcnf  (not pol) e1) (xcnf  pol e2) in
 
   xcnf  true f
 
-
-
 module M =
 struct
  open Coqlib
@@ -415,11 +438,11 @@ let parse_q term =
    Printf.fprintf o "%s%a%s" op _pp l cl 
 
 
-
  let pp_var  = pp_positive
  let dump_var = dump_positive
 
- let rec pp_expr o e = 
+ let pp_expr pp_z o e = 
+ let rec pp_expr  o e = 
   match e with
    | Mc.PEX n -> Printf.fprintf o "V %a" pp_var n
    | Mc.PEc z -> pp_z o z
@@ -427,7 +450,8 @@ let parse_q term =
    | Mc.PEmul(e1,e2) -> Printf.fprintf o "%a*(%a)" pp_expr e1 pp_expr e2
    | Mc.PEopp e -> Printf.fprintf o "-(%a)" pp_expr e
    | Mc.PEsub(e1,e2) -> Printf.fprintf o "(%a)-(%a)" pp_expr e1 pp_expr e2
-   | Mc.PEpow(e,n) -> Printf.fprintf o "(%a)^(%a)" pp_expr e pp_n  n
+   | Mc.PEpow(e,n) -> Printf.fprintf o "(%a)^(%a)" pp_expr e pp_n  n in
+   pp_expr o e
 
 
  let  dump_expr typ dump_z e =
@@ -476,9 +500,9 @@ let parse_q term =
     | Mc.S_In n -> 
        Printf.fprintf o "(S_In %a)%%nat" pp_nat n
     | Mc.S_Ideal(e,c) -> 
-       Printf.fprintf o "(S_Ideal %a %a)" pp_expr e pp_cone c
+       Printf.fprintf o "(S_Ideal %a %a)" (pp_expr pp_z) e pp_cone c
     | Mc.S_Square e -> 
-       Printf.fprintf o "(S_Square %a)" pp_expr e
+       Printf.fprintf o "(S_Square %a)" (pp_expr pp_z) e
     | Mc.S_Monoid l -> 
        Printf.fprintf o "(S_Monoid %a)" (pp_list "[" "]" pp_nat) l
     | Mc.S_Add(e1,e2) -> 
@@ -514,8 +538,8 @@ let parse_q term =
 
 
 
- let pp_cstr o {Mc.flhs = l ; Mc.fop = op ; Mc.frhs = r } =
-  Printf.fprintf o"(%a %a %a)" pp_expr l pp_op op pp_expr r
+ let pp_cstr pp_z o {Mc.flhs = l ; Mc.fop = op ; Mc.frhs = r } =
+  Printf.fprintf o"(%a %a %a)" (pp_expr pp_z) l pp_op op (pp_expr pp_z) r
 
  let dump_cstr typ dump_constant {Mc.flhs = e1 ; Mc.fop = o ; Mc.frhs = e2} =
    Term.mkApp(Lazy.force coq_Build,
@@ -743,21 +767,21 @@ let parse_rexpr =
   | TT  -> Mc.TT
   | FF  -> Mc.FF
   | X _  -> Mc.X
-  | A (x,_) -> Mc.A x
-  | C (a,b,_)  -> Mc.Cj(f2f a,f2f b)
-  | D (a,b,_)  -> Mc.D(f2f a,f2f b)
-  | N (a,_) -> Mc.N(f2f a)
-  | I(a,b,_)  -> Mc.I(f2f a,f2f b)
+  | A (x,_,_) -> Mc.A x
+  | C (a,b)  -> Mc.Cj(f2f a,f2f b)
+  | D (a,b)  -> Mc.D(f2f a,f2f b)
+  | N (a) -> Mc.N(f2f a)
+  | I(a,_,b)  -> Mc.I(f2f a,f2f b)
 
  let is_prop t = 
   match t with
    | Names.Anonymous -> true (* Not quite right *)
    | Names.Name x    -> false
 
- let mkC f1 f2 = C(f1,f2,none)
- let mkD f1 f2 = D(f1,f2,none)
- let mkIff f1 f2 = C(I(f1,f2,none),I(f2,f2,none),none)
- let mkI f1 f2 = I(f1,f2,none)
+ let mkC f1 f2 = C(f1,f2)
+ let mkD f1 f2 = D(f1,f2)
+ let mkIff f1 f2 = C(I(f1,None,f2),I(f2,None,f1))
+ let mkI f1 f2 = I(f1,None,f2)
 
  let mkformula_binary g term f1 f2 =
    match f1 , f2 with
@@ -765,34 +789,36 @@ let parse_rexpr =
    |   _         -> g f1 f2
 
  let parse_formula parse_atom env term =
-  let parse_atom env t = try let (at,env) = parse_atom env t in (A(at,none), env) with _ -> (X(t),env) in
 
-  let rec xparse_formula env term =
+  let parse_atom env tg t = try let (at,env) = parse_atom env t in 
+				  (A(at,tg,t), env,tg+1) with _ -> (X(t),env,tg) in
+
+  let rec xparse_formula env tg term =
    match kind_of_term term with
     | App(l,rst) ->
         (match rst with
          | [|a;b|] when l = Lazy.force coq_and ->
-	     let f,env = xparse_formula env a in
-	     let g,env = xparse_formula env b in
-             mkformula_binary mkC term f g,env
+	     let f,env,tg = xparse_formula env tg a in
+	     let g,env, tg = xparse_formula env tg b  in
+             mkformula_binary mkC term f g,env,tg
          | [|a;b|] when l = Lazy.force coq_or ->
-	     let f,env = xparse_formula env a in
-	     let g,env = xparse_formula env b in
-             mkformula_binary mkD term f g,env
+	     let f,env,tg = xparse_formula env tg a in
+	     let g,env,tg  = xparse_formula env tg b in
+             mkformula_binary mkD term f g,env,tg
          | [|a|] when l = Lazy.force coq_not ->
-             let (f,env) = xparse_formula env a in (N(f,none), env)
+             let (f,env,tg) = xparse_formula env tg a in (N(f), env,tg)
          | [|a;b|] when l = Lazy.force coq_iff ->
-	     let f,env = xparse_formula env a in
-	     let g,env = xparse_formula env b in
-             mkformula_binary mkIff term f g,env
-         | _ -> parse_atom env term)
+	     let f,env,tg = xparse_formula env tg a in
+	     let g,env,tg = xparse_formula env tg b in
+             mkformula_binary mkIff term f g,env,tg
+         | _ -> parse_atom env tg term)
     | Prod(typ,a,b) when not (Termops.dependent (mkRel 1) b) ->
-	let f,env = xparse_formula env a in
-	let g,env = xparse_formula env b in
-        mkformula_binary mkI term f g,env
-    | _ when term = Lazy.force coq_True -> (TT,env)
-    | _ when term = Lazy.force coq_False -> (FF,env)
-    | _  -> X(term),env in
+	let f,env,tg = xparse_formula env tg a in
+	let g,env,tg = xparse_formula env tg b in
+        mkformula_binary mkI term f g,env,tg
+    | _ when term = Lazy.force coq_True -> (TT,env,tg)
+    | _ when term = Lazy.force coq_False -> (FF,env,tg)
+    | _  -> X(term),env,tg in
   xparse_formula env term
 
  let coq_TT = lazy 
@@ -828,11 +854,11 @@ let parse_rexpr =
    match f with
     | TT  -> mkApp(Lazy.force coq_TT,[| typ|])
     | FF  -> mkApp(Lazy.force coq_FF,[| typ|])
-    | C(x,y,_) -> mkApp(Lazy.force coq_And,[| typ ; xdump x ; xdump y|])
-    | D(x,y,_) -> mkApp(Lazy.force coq_Or,[| typ ; xdump x ; xdump y|])
-    | I(x,y,_) -> mkApp(Lazy.force coq_Impl,[| typ ; xdump x ; xdump y|])
-    | N(x,_) -> mkApp(Lazy.force coq_Neg,[| typ ; xdump x|])
-    | A(x,_) -> mkApp(Lazy.force coq_Atom,[| typ ; dump_atom x|]) 
+    | C(x,y) -> mkApp(Lazy.force coq_And,[| typ ; xdump x ; xdump y|])
+    | D(x,y) -> mkApp(Lazy.force coq_Or,[| typ ; xdump x ; xdump y|])
+    | I(x,_,y) -> mkApp(Lazy.force coq_Impl,[| typ ; xdump x ; xdump y|])
+    | N(x) -> mkApp(Lazy.force coq_Neg,[| typ ; xdump x|])
+    | A(x,_,_) -> mkApp(Lazy.force coq_Atom,[| typ ; dump_atom x|]) 
     | X(t) -> mkApp(Lazy.force coq_X,[| typ ; t|])  in
 
    xdump f
@@ -978,22 +1004,22 @@ let pp_q o q = Printf.fprintf o "%a/%a" pp_z q.Micromega.qnum pp_positive q.Micr
  
  
 let rec pp_proof_term o = function
- | Micromega.RatProof cone -> Printf.fprintf o "R[%a]" (pp_cone pp_z) cone
+ | Micromega.RatProof cone -> Printf.fprintf o "R[%a]" (pp_cone  pp_z) cone
  | Micromega.CutProof(e,q,_,p) -> failwith "not implemented"
  | Micromega.EnumProof(q1,e1,q2,c1,c2,rst) -> 
     Printf.fprintf o "EP[%a,%a,%a,%a,%a,%a]" 
-     pp_q q1 pp_expr e1 pp_q q2 (pp_cone pp_z) c1 (pp_cone pp_z) c2 
+     pp_q q1 (pp_expr pp_z) e1 pp_q q2 (pp_cone pp_z) c1 (pp_cone pp_z) c2 
      (pp_list "[" "]" pp_proof_term) rst
 
-let rec parse_hyps parse_arith env hyps =
+let rec parse_hyps parse_arith env tg hyps =
  match hyps with
-  | [] -> ([],env)
+  | [] -> ([],env,tg)
   | (i,t)::l -> 
-     let (lhyps,env) = parse_hyps parse_arith env l in
+     let (lhyps,env,tg) = parse_hyps parse_arith env tg l in
       try 
-       let (c,env) = parse_formula parse_arith env  t in
-	((i,c)::lhyps, env)
-      with _ -> 		(lhyps,env)
+       let (c,env,tg) = parse_formula parse_arith env  tg t in
+	((i,c)::lhyps, env,tg)
+      with _ -> 		(lhyps,env,tg)
        (*(if debug then Printf.printf "parse_arith : %s\n" x);*)
 
 
@@ -1001,18 +1027,18 @@ exception ParseError
 
 let parse_goal parse_arith env hyps term = 
  (*  try*)
- let (f,env) = parse_formula parse_arith env term in
- let (lhyps,env) = parse_hyps parse_arith env hyps in
+ let (f,env,tg) = parse_formula parse_arith env 0 term in
+ let (lhyps,env,tg) = parse_hyps parse_arith env tg hyps in
   (lhyps,f,env)
    (*  with Failure x -> raise ParseError*)
 
 
-type ('a, 'b) domain_spec = { 
+type ('d, 'prf) domain_spec = { 
  typ : Term.constr; (* Z, Q , R *)
  coeff : Term.constr ; (* Z, Q *)
- dump_coeff : 'a -> Term.constr ; 
+ dump_coeff : 'd -> Term.constr ; 
  proof_typ : Term.constr ; 
- dump_proof : 'b -> Term.constr
+ dump_proof : 'prf -> Term.constr
 }
 
 let zz_domain_spec  = lazy {
@@ -1040,6 +1066,34 @@ let rz_domain_spec = lazy {
 }
 
 
+let abstract_formula hyps f = 
+  
+  let rec xabs f = 
+    match f with
+      | X c -> X c
+      | A(a,t,term) -> if ISet.mem t hyps then A(a,t,term) else X(term)
+      | C(f1,f2) -> 
+	  (match xabs f1 , xabs f2 with
+	    |   X a1    ,  X a2   -> X (Term.mkApp(Lazy.force coq_and, [|a1;a2|]))
+	    |    f1     , f2      -> C(f1,f2) )
+      | D(f1,f2) -> 
+	  (match xabs f1 , xabs f2 with
+	    |   X a1    ,  X a2   -> X (Term.mkApp(Lazy.force coq_or, [|a1;a2|]))
+	    |    f1     , f2      -> D(f1,f2) )
+      | N(f) -> 
+	  (match xabs f with
+	    |   X a    -> X (Term.mkApp(Lazy.force coq_not, [|a|]))
+	    |     f     -> N f)
+      | I(f1,hyp,f2) -> 
+	  (match xabs f1 , hyp, xabs f2 with
+	    | X a1      , Some _ , af2    ->  af2
+	    | X a1      , None   , X a2   -> X (Term.mkArrow a1 a2)
+	    |   af1     ,  _     , af2    -> I(af1,hyp,af2)
+	  )
+      | x -> x in
+    xabs f
+
+
 
 
 let micromega_order_change spec cert cert_typ env  ff gl = 
@@ -1062,42 +1116,27 @@ let micromega_order_change spec cert cert_typ env  ff gl =
    gl 
    
 
-let detect_duplicates cnf wit = 
- let cnf = CoqToCaml.list (fun x -> x) cnf in
- let wit = CoqToCaml.list (fun x -> x) wit in
-
- let rec xdup cnf wit = 
-  match wit with
-   | [] -> []
-   | w :: wit -> 
-      let sg = sig_of_cone w in
-       match cnf with
-	| [] -> []
-	| e::cnf -> 
-	   let (dups,cnf) = (List.partition (fun x -> same_proof sg e x) cnf) in
-	    dups@(xdup cnf wit) in
-  xdup cnf wit
+type ('a,'prf) prover = {
+  name : string ; (* name of the prover *)
+  prover : 'a list -> 'prf option; (* the prover itself *)
+  hyps : 'prf -> ISet.t; (* extract the indexes of the hypotheses really used in the proof *)
+  compact : 'prf -> (int -> int) -> 'prf; (* remap the hyp indexes according to function *)
+  pp_prf : out_channel -> 'prf -> unit ;(* pretting printing of proof *)
+  pp_f   : out_channel -> 'a   -> unit
+}
 
-let find_witness prover    polys1 =
- try_any prover polys1
+let find_witness provers    polys1 =
+  
+  let provers = List.map (fun p -> 
+    (fun l -> 
+      match p.prover l with
+        | None -> None
+        | Some prf -> Some(prf,p)) , p.name) provers in
+                           
+    try_any provers (List.map fst polys1)
 
 
-let witness_list_with_tags prover l = 
- 
- let rec xwitness_list l = 
-  match l with
-   | [] -> Some([])
-   | e::l -> 
-      match find_witness prover  (List.map fst e)  with
-       | None -> None
-       | Some w -> 
-	  (match xwitness_list l with
-	   | None -> None
-	   | Some l -> Some (w::l)
-	  ) in
-  xwitness_list l
-
-let witness_list_without_tags prover l = 
+let witness_list_tags prover l = 
  
  let rec xwitness_list l = 
   match l with
@@ -1127,10 +1166,43 @@ let witness_list prover l =
   xwitness_list l
 
 
+let is_singleton = function [] -> true | [e] -> true | _ -> false
+
+let pp_ml_list pp_elt o l = 
+  output_string o "[" ;
+  List.iter (fun x -> Printf.fprintf o "%a ;" pp_elt x) l ; 
+  output_string o "]" 
+
+let compact_proofs (cnf_ff: 'cst cnf) res (cnf_ff': 'cst cnf) = 
+
+  let compact_proof (old_cl:'cst clause) (prf,prover) (new_cl:'cst clause) = 
+    let new_cl = Mutils.mapi (fun (f,_) i -> (f,i)) new_cl in
+    let remap i = 
+      (*Printf.printf "nth (len:%i) %i\n" (List.length old_cl) i;*)
+      let formula = try fst (List.nth old_cl i) with Failure _ -> failwith "bad old index" in
+	List.assoc formula new_cl  in
+    let res = prover.compact prf remap  in
+
+      if debug then 
+	begin
+	  Printf.printf "\ncompact_proof : %a %a %a -> %a\n" 
+	    (pp_ml_list prover.pp_f) (List.map fst old_cl) 
+	    prover.pp_prf prf 
+	    (pp_ml_list prover.pp_f) (List.map fst new_cl)  
+	    prover.pp_prf res ;
+	  flush stdout
+	end
+	  ; 
+      res in
+
+    (* This seems overkill when List.length cnf_ff = List.length cnf_ff'. Isn't there a 1-1-mapping *)
+  let cnf_res = List.combine cnf_ff res in
+  let l = 
+    List.map (fun x -> let (o,p) = List.find (fun lp -> is_sublist x (fst lp)) cnf_res in (o,p,x)) cnf_ff' in
+
+   List.map (fun (o,p,n) -> compact_proof o p n) l
 
 
-let is_singleton = function [] -> true | [e] -> true | _ -> false
- 
 
 let micromega_tauto negate normalise spec prover env polys1 polys2  gl = 
  let spec = Lazy.force spec in 
@@ -1139,34 +1211,63 @@ let micromega_tauto negate normalise spec prover env polys1 polys2  gl =
    (fun (id,f)  (cc,ids) -> 
     match f with 
       X _ -> (cc,ids) 
-     | _ -> (I(tag_formula (Names.Name id) f,cc,none), id::ids)) 
+     | _ -> (I(f,Some id,cc), id::ids)) 
    polys1 (polys2,[]) in
 
  let cnf_ff = cnf negate normalise ff in
 
   if debug then 
-   (Pp.pp (Pp.str "Formula....\n") ;
-     let formula_typ = (Term.mkApp( Lazy.force coq_Cstr,[| spec.coeff|])) in
-     let ff = dump_formula formula_typ 
-      (dump_cstr spec.typ spec.dump_coeff) ff in
-      Pp.pp (Printer.prterm  ff) ;  Pp.pp_flush ()) ;
-
-   match witness_list_without_tags prover  cnf_ff with
+    begin
+      Pp.pp (Pp.str "Formula....\n") ;
+      let formula_typ = (Term.mkApp( Lazy.force coq_Cstr,[| spec.coeff|])) in
+      let ff = dump_formula formula_typ 
+        (dump_cstr spec.typ spec.dump_coeff) ff in
+        Pp.pp (Printer.prterm  ff) ;  Pp.pp_flush ()
+    end;
+
+   match witness_list_tags prover  cnf_ff with
     | None -> Tacticals.tclFAIL 0 (Pp.str "Cannot find witness") gl
     | Some res -> (*Printf.printf "\nList %i" (List.length res); *)
-       let (ff,res,ids) = (ff,res,List.map Term.mkVar ids) in
-       let res' = dump_ml_list (spec.proof_typ) spec.dump_proof  res in
-	(Tacticals.tclTHENSEQ
-	  [
-	   Tactics.generalize ids;
-	   micromega_order_change spec res' 
-	    (Term.mkApp(Lazy.force coq_list,[|  spec.proof_typ|]))  env  ff ;
-	   ]) gl
-
-
-let micromega_gen parse_arith negate normalise spec  prover   gl =
- let concl = Tacmach.pf_concl gl in
- let hyps  = Tacmach.pf_hyps_types gl in
+        let hyps = List.fold_left (fun s (cl,(prf,p)) -> 
+	  let tags = ISet.fold (fun i s -> try ISet.add (snd (List.nth cl i)) s  with Invalid_argument _ -> s) (p.hyps prf) ISet.empty in
+	  ISet.union s tags) ISet.empty (List.combine cnf_ff res) in
+
+	  if debug then (Printf.printf "TForm : %a\n" pp_formula ff ; flush stdout;
+			 Printf.printf "Hyps : %a\n" (fun o s -> ISet.fold (fun i _ -> Printf.fprintf o "%i " i) s ()) hyps) ;  
+	  
+        let ff'   = abstract_formula hyps ff in
+	  
+	  if debug then
+	    begin
+	      Pp.pp (Pp.str "\nAFormula\n") ; 
+	      let formula_typ = (Term.mkApp( Lazy.force coq_Cstr,[| spec.coeff|])) in
+	      let ff' = dump_formula formula_typ 
+		(dump_cstr spec.typ spec.dump_coeff) ff' in
+		Pp.pp (Printer.prterm  ff') ;  Pp.pp_flush ()
+	    end;
+
+
+        let cnf_ff' = cnf negate normalise ff' in
+        let res' = compact_proofs cnf_ff res cnf_ff' in
+
+        let (ff',res',ids) = (ff',res',List.map Term.mkVar (ids_of_formula ff')) in
+
+       let res' = dump_ml_list (spec.proof_typ) spec.dump_proof   res' in
+        (Tacticals.tclTHENSEQ
+          [
+           Tactics.generalize ids;
+           micromega_order_change spec res' 
+            (Term.mkApp(Lazy.force coq_list,[|  spec.proof_typ|]))  env  ff' ;
+           ]) gl
+
+
+let micromega_gen 
+    parse_arith 
+    (negate:'cst atom -> 'cst mc_cnf) 
+    (normalise:'cst atom -> 'cst mc_cnf) 
+    spec  prover   gl =
+  let concl = Tacmach.pf_concl gl in
+  let hyps  = Tacmach.pf_hyps_types gl in
   try
    let (hyps,concl,env) = parse_goal parse_arith Env.empty hyps concl in
    let env = Env.elements env in
@@ -1235,52 +1336,174 @@ let call_csdpcert_z provername poly =
 	| Mc.False ->  (print_string "buggy certificate" ; flush stdout) ;None
 
 
+let hyps_of_cone  prf = 
+  let rec xtract e acc = 
+    match e with
+    | Mc.S_Pos _ | Mc.S_Z | Mc.S_Square _ -> acc
+    | Mc.S_In n -> ISet.add (CoqToCaml.nat n) acc
+    | Mc.S_Ideal(_,c) -> xtract  c acc
+    | Mc.S_Monoid Mc.Nil   -> acc
+    | Mc.S_Monoid (Mc.Cons(i,l)) -> xtract (Mc.S_Monoid l) (ISet.add (CoqToCaml.nat i) acc)
+    | Mc.S_Add(e1,e2) |  Mc.S_Mult(e1,e2) -> xtract e1 (xtract e2 acc) in
+
+    xtract prf ISet.empty
+
+let compact_cone prf f  =
+  let np n = CamlToCoq.nat (f (CoqToCaml.nat n)) in
+
+  let rec xinterp prf = 
+    match prf with
+    | Mc.S_Pos _ | Mc.S_Z | Mc.S_Square _ -> prf
+    | Mc.S_In n -> Mc.S_In (np n)
+    | Mc.S_Ideal(e,c) -> Mc.S_Ideal(e,xinterp c)
+    | Mc.S_Monoid l  -> Mc.S_Monoid (Mc.map np l)
+    | Mc.S_Add(e1,e2) -> Mc.S_Add(xinterp e1,xinterp e2)
+    | Mc.S_Mult(e1,e2) -> Mc.S_Mult(xinterp e1,xinterp e2) in
+
+    xinterp prf
+
+
+let hyps_of_pt pt = 
+  let translate x s = ISet.fold (fun i s -> ISet.add (i - x) s) s ISet.empty in
+
+  let rec xhyps ofset pt acc = 
+    match pt with
+      | Mc.RatProof p -> ISet.union acc (translate ofset (hyps_of_cone p))
+      | Mc.CutProof(_,_,c,pt) -> xhyps (ofset+1) pt (ISet.union (translate (ofset+1) (hyps_of_cone c)) acc)
+      | Mc.EnumProof(_,_,_,c1,c2,l) -> 
+	  let l = CoqToCaml.list (fun x -> x) l in
+	  let s = ISet.union acc (translate (ofset +1) (ISet.union (hyps_of_cone c1) (hyps_of_cone c2))) in
+	    List.fold_left (fun s x -> xhyps (ofset +  1) x s) s l in
+    
+    xhyps 0 pt ISet.empty
+	  
+  
+let compact_pt pt f = 
+  let translate ofset x =
+    if x < ofset then x
+    else (f (x-ofset)  + ofset) in
+
+  let rec compact_pt ofset pt = 
+    match pt with
+      | Mc.RatProof p -> Mc.RatProof (compact_cone p (translate ofset))
+      | Mc.CutProof(x,y,c,pt) -> Mc.CutProof(x,y,compact_cone c (translate (ofset+1)), compact_pt (ofset+1) pt )
+      | Mc.EnumProof(a,b,c,c1,c2,l) -> Mc.EnumProof(a,b,c,compact_cone c1 (translate (ofset+1)), compact_cone c2 (translate (ofset+1)),
+						   Mc.map (fun x -> compact_pt (ofset+1) x) l) in
+    compact_pt 0 pt
+
+(** Definition of provers *)
+      
+let linear_prover_Z = {
+  name    = "linear prover" ; 
+  prover  = lift_ratproof (Certificate.linear_prover Certificate.z_spec) ; 
+  hyps    = hyps_of_pt ;
+  compact = compact_pt ;
+  pp_prf  = pp_proof_term;
+  pp_f    = fun o x -> pp_expr pp_z o (fst' x)
+}
+
+let linear_prover_Q = {
+  name    = "linear prover";
+  prover  = (Certificate.linear_prover Certificate.q_spec) ; 
+  hyps    = hyps_of_cone ; 
+  compact = compact_cone ;
+  pp_prf  = pp_cone pp_q ;
+  pp_f   = fun o x -> pp_expr pp_q o  (fst' x)
+}
+
+let linear_prover_R = {
+  name    = "linear prover";
+  prover  = (Certificate.linear_prover Certificate.z_spec) ; 
+  hyps    = hyps_of_cone ; 
+  compact = compact_cone ;
+  pp_prf  = pp_cone pp_z ;
+  pp_f    =  fun o x -> pp_expr pp_z o (fst' x)
+}
+
+let non_linear_prover_Q str o = {
+  name    = "real nonlinear prover";
+  prover  = call_csdpcert_q (str, o);
+  hyps    = hyps_of_cone;
+  compact = compact_cone ;
+  pp_prf  = pp_cone pp_q ;
+  pp_f    = fun o x -> pp_expr pp_q o  (fst' x)
+}
+
+let non_linear_prover_R str o = {
+  name    = "real nonlinear prover";
+  prover  = call_csdpcert_z (str, o);
+  hyps    = hyps_of_cone;
+  compact = compact_cone;
+  pp_prf  = pp_cone pp_z;
+  pp_f    = fun o x -> pp_expr pp_z o  (fst' x)
+}
+
+let non_linear_prover_Z str o  = {
+  name    = "real nonlinear prover";
+  prover  = lift_ratproof (call_csdpcert_z (str, o));
+  hyps    = hyps_of_pt;
+  compact = compact_pt;
+  pp_prf  = pp_proof_term;
+  pp_f    =  fun o x -> pp_expr pp_z o (fst' x)
+}
+
+let linear_Z =   {
+  name = "lia";
+  prover = Certificate.zlinear_prover ; 
+  hyps   = hyps_of_pt;
+  compact = compact_pt;
+  pp_prf      = pp_proof_term;
+  pp_f    = fun o x -> pp_expr pp_z o (fst' x)
+}
+
+
 
+(** Instantiation of the tactics *)
 
 let psatzl_Z gl = 
  micromega_gen parse_zarith  Mc.negate Mc.normalise zz_domain_spec
-  [lift_ratproof 
-    (Certificate.linear_prover Certificate.z_spec), "fourier refutation" ] gl
-
+  [linear_prover_Z ] gl
 
 let psatzl_Q gl = 
  micromega_gen  parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec 
-  [ Certificate.linear_prover Certificate.q_spec, "fourier refutation" ]   gl
+  [ linear_prover_Q ]   gl
+
 
 let psatz_Q i gl = 
  micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec
-  [ call_csdpcert_q ("real_nonlinear_prover", Some i), "fourier refutation" ]  gl
+  [ non_linear_prover_Q "real_nonlinear_prover" (Some i) ]  gl
+
 
 let psatzl_R gl = 
  micromega_gen  parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec 
-  [ Certificate.linear_prover Certificate.z_spec, "fourier refutation" ]   gl
-
+  [ linear_prover_R ]   gl
 
 let psatz_R i gl = 
  micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec
-  [ call_csdpcert_z  ("real_nonlinear_prover", Some i), "fourier refutation" ]  gl
+  [ non_linear_prover_R "real_nonlinear_prover" (Some i)]  gl
+
 
-  
 let psatz_Z i gl = 
  micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec 
- [lift_ratproof (call_csdpcert_z ("real_nonlinear_prover",Some i)), 
-  "fourier refutation" ] gl
+ [non_linear_prover_Z "real_nonlinear_prover" (Some i) ] gl
 
 
 let sos_Z gl = 
  micromega_gen parse_zarith Mc.negate Mc.normalise  zz_domain_spec 
-  [lift_ratproof (call_csdpcert_z ("pure_sos", None)), "pure sos refutation"]  gl
+  [non_linear_prover_Z "pure_sos" None]  gl
 
 let sos_Q gl = 
  micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise  qq_domain_spec 
-  [call_csdpcert_q ("pure_sos", None), "pure sos refutation"]  gl
+  [non_linear_prover_Q "pure_sos" None]  gl
+
 
 let sos_R gl = 
  micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise  rz_domain_spec 
-  [call_csdpcert_z ("pure_sos", None), "pure sos refutation"]  gl
+  [non_linear_prover_R "pure_sos" None]  gl
 
 
 
 let xlia gl = 
  micromega_gen parse_zarith Mc.negate Mc.normalise   zz_domain_spec 
-  [Certificate.zlinear_prover, "zprover"] gl
+  [linear_Z] gl
+
diff --git a/contrib/micromega/mfourier.ml b/contrib/micromega/mfourier.ml
index 415d3a3e2..c2eb572c9 100644
--- a/contrib/micromega/mfourier.ml
+++ b/contrib/micromega/mfourier.ml
@@ -1,667 +1,1014 @@
-(************************************************************************)
-(*  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        *)
-(************************************************************************)
-(*                                                                      *)
-(* Micromega: A reflexive tactic using the Positivstellensatz           *)
-(*                                                                      *)
-(*  Frédéric Besson (Irisa/Inria) 2006-2008                             *)
-(*                                                                      *)
-(************************************************************************)
-
-(* Yet another implementation of Fourier *)
 open Num
+module Utils = Mutils
+
+let map_option = Utils.map_option
+let from_option = Utils.from_option
+
+let debug = false
+type ('a,'b) lr = Inl of 'a | Inr of 'b
+
+
+module Vect = 
+  struct 
+    (** [t] is the type of vectors.
+	A vector [(x1,v1) ; ... ; (xn,vn)] is such that:
+	- variables indexes are ordered (x1 < ... < xn
+	- values are all non-zero
+    *)
+    type var = int
+    type t = (var * num) list 
+
+(** [equal v1 v2 = true] if the vectors are syntactically equal. 
+    ([num] is not handled by  [Pervasives.equal] *)
+
+    let rec equal v1 v2 = 
+      match v1 , v2 with
+	| [] , []   -> true
+	| [] , _    -> false
+	| _::_ , [] -> false
+	| (i1,n1)::v1 , (i2,n2)::v2 -> 
+	    (i1 = i2) && n1 =/ n2 && equal v1 v2
+
+    let hash v = 
+      let rec hash i = function 
+	| [] -> i
+	| (vr,vl)::l ->  hash (i + (Hashtbl.hash (vr, float_of_num vl))) l in
+	Hashtbl.hash (hash 0 v )
 
-module Cmp =
- (* How to compare pairs, lists ... *)
-struct
- let rec compare_lexical l = 
-  match l with
-   | [] -> 0 (* Equal *)
-   | f::l -> 
-      let cmp = f () in
-       if cmp = 0 	then  compare_lexical l else cmp
 	
- let rec compare_list cmp l1 l2 = 
-  match l1 , l2 with
-   | []  , [] -> 0
-   | []  , _  -> -1
-   | _   , [] -> 1
-   | e1::l1 , e2::l2 -> 
-      let c = cmp e1 e2 in
-       if c = 0 then compare_list cmp l1 l2 else c
+    let null = []
+
+    let pp_vect o vect = 
+      List.iter (fun (v,n) -> Printf.printf "%sx%i + " (string_of_num n) v) vect
+	
+    let from_list (l: num list) = 
+      let rec xfrom_list i l = 
+	match l with
+	  | [] -> []
+	  | e::l ->  
+	      if e <>/ Int 0 
+	      then (i,e)::(xfrom_list (i+1) l)
+	      else xfrom_list (i+1) l in
 	
- let hash_list hash l = 
-  let rec xhash res l = 
-   match l with
-    | []  -> res 
-    | e::l -> xhash ((hash e)  lxor res) l in
-   xhash  (Hashtbl.hash []) l
+	xfrom_list 0 l
 
-end
- 
-module Interval = 
-struct
- (** The type of intervals.  **)
- type intrvl = Empty | Point of num | Itv  of  num option * num option
+    let zero_num = Int 0
+    let unit_num = Int 1
+      
+      
+    let to_list m = 
+      let rec xto_list i l =
+	match l with
+	  | [] -> []
+	  |	(x,v)::l' -> 
+		  if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in
+	xto_list 0 m
+	  
+
+    let cons i v rst = if v =/ Int 0 then rst else (i,v)::rst
+	
+    let rec update i f t = 
+      match t with
+	| [] -> cons i (f zero_num) []
+	| (k,v)::l -> 
+	    match Pervasives.compare i k with
+	  | 0 -> cons k (f v) l
+	  | -1 -> cons i (f zero_num) t
+	  |  1 -> (k,v) ::(update i f l)
+	  |  _  -> failwith "compare_num"
+	       
+    let rec set i n t =
+      match t with
+	| [] -> cons i n []
+	| (k,v)::l -> 
+	    match Pervasives.compare i k with
+	      | 0 -> cons k n l
+	      | -1 -> cons i n t
+	      |  1 -> (k,v) :: (set i n l)
+	      |  _  -> failwith "compare_num"
+		   
+    let gcd m = 
+      let res = List.fold_left (fun x (i,e) -> Big_int.gcd_big_int  x (Utils.numerator e))  Big_int.zero_big_int m in
+	if Big_int.compare_big_int res Big_int.zero_big_int = 0 
+	then Big_int.unit_big_int else res
+	  
+    let rec mul z t = 
+      match z with
+	| Int 0 -> []
+	| Int 1 -> t
+	|  _    -> List.map (fun (i,n) -> (i, mult_num z n)) t
+
+    let compare : t -> t -> int = Utils.Cmp.compare_list (fun x y -> Utils.Cmp.compare_lexical 
+      [
+	(fun () -> Pervasives.compare (fst x) (fst y));
+	(fun () -> compare_num (snd x) (snd y))]) 
+
+    (** [tail v vect] returns
+	- [None] if [v] is not a variable of the vector [vect]
+	- [Some(vl,rst)]  where [vl] is the value of [v] in vector [vect]
+        and [rst] is the remaining of the vector
+	We exploit that vectors are ordered lists
+    *)
+    let rec tail (v:var) (vect:t) = 
+      match vect with
+	| [] -> None
+	| (v',vl)::vect' -> 
+	    match Pervasives.compare v' v with
+	      | 0 -> Some (vl,vect) (* Ok, found *)
+	      | -1 -> tail v vect' (* Might be in the tail *)
+	      | _  -> None (* Hopeless *)
+		  
+    let get v vect = 
+      match tail v vect with
+	| None -> None
+	| Some(vl,_) -> Some vl
+
+
+    let rec fresh v =
+	match v with
+	  | [] -> 1
+	  | [v,_] -> v + 1
+	  | _::v  -> fresh v
+
+  end
+open Vect
+
+(** Implementation of intervals *)
+module Itv = 
+struct 
   
- (**
-    Different intervals can denote the same set of variables e.g., 
-    Point n && Itv (Some n, Some n)
-    Itv (Some x) (Some y) && Empty if x > y
-    see the 'belongs_to' function.
- **)
-
- (* The set of numerics that belong to an interval *)
- let belongs_to n = function
-  | Empty -> false
-  | Point x -> n =/ x
-  | Itv(Some x, Some y) -> x <=/ n && n <=/ y
-  | Itv(None,Some y) -> n <=/ y
-  | Itv(Some x,None) -> x <=/ n
-  | Itv(None,None) -> true
-
- let string_of_bound = function
-  | None -> "oo"
-  | Some n -> Printf.sprintf "Bd(%s)"  (string_of_num n)
-
- let string_of_intrvl = function
-  | Empty -> "[]"
-  | Point n -> Printf.sprintf "[%s]" (string_of_num n)
-  | Itv(bd1,bd2) -> 
-     Printf.sprintf "[%s,%s]" (string_of_bound bd1) (string_of_bound bd2)
-
- let pick_closed_to_zero = function
-  | Empty -> None
-  | Point n -> Some n
-  | Itv(None,None) -> Some (Int 0)
-  | Itv(None,Some i) -> 
-     Some (if  (Int 0) <=/ (floor_num  i) then Int 0 else floor_num i)
-  | Itv(Some i,None) -> 
-     Some (if i <=/ (Int 0) then Int 0 else ceiling_num i)
-  | Itv(Some i,Some j) -> 
-     Some (
-      if i <=/ Int 0 && Int 0 <=/ j 
-      then Int 0
-      else if ceiling_num i <=/ floor_num j 
-      then ceiling_num i (* why not *) else i)
-
- type status = 
-   | O | Qonly | Z | Q 
-       
- let interval_kind = function
-  | Empty -> O
-  | Point n -> if ceiling_num n =/ n then Z else Qonly
-  | Itv(None,None) -> Z
-  | Itv(None,Some i) -> if ceiling_num i <>/ i then Q else Z
-  | Itv(Some i,None) -> if ceiling_num i <>/ i then Q else Z
-  | Itv(Some i,Some j) -> 
-     if ceiling_num i <>/ i or floor_num j <>/ j then Q else Z
-
- let empty_z = function
-  | Empty -> true
-  | Point n ->  ceiling_num n <>/ n 
-  | Itv(None,None) | Itv(None,Some _) | Itv(Some _,None) -> false
-  | Itv(Some i,Some j) -> ceiling_num i >/ floor_num j
-
-
- let normalise b1 b2 =
-  match b1 , b2 with
-   | Some i , Some j -> 
-      (match compare_num i j with
-       | 1 -> Empty
-       | 0 -> Point i
-       | _ -> Itv(b1,b2)
-      )
-   |  _ -> Itv(b1,b2)
-
-       
-
- let min x y = 
-  match x , y with
-   | None , x | x , None -> x
-   | Some i , Some j -> Some (min_num i j)
+  (** The type of intervals is *)
+  type interval = num option * num option
+      (** None models the absence of bound i.e. infinity *)
+      (** As a result, 
+	  - None , None   -> ]-oo,+oo[
+	  - None , Some v -> ]-oo,v]
+	  - Some v, None  -> [v,+oo[
+	  - Some v, Some v' -> [v,v']
+      Intervals needs to be explicitely normalised.
+      *)
+
+  type who = Left | Right 
+
+
+  (** if then interval [itv] is empty, [norm_itv itv] returns [None] 
+      otherwise, it returns [Some itv] *)
+      
+  let norm_itv itv = 
+    match itv with
+      | Some a , Some b -> if a <=/ b then Some itv else None
+      |   _             -> Some itv
+	    
+  (** [opp_itv itv] computes the opposite interval *)
+  let opp_itv itv = 
+  let (l,r) = itv in
+    (map_option minus_num r, map_option minus_num l)
+      
+
+
+
+(** [inter i1 i2 = None] if the intersection of intervals is empty
+    [inter i1 i2 = Some i] if [i] is the intersection of the intervals [i1] and [i2] *)
+  let inter i1 i2 = 
+    let (l1,r1) = i1
+    and (l2,r2) = i2 in
       
- let max x y = 
-  match x , y with
-   | None , x | x , None -> x
-   | Some i , Some j -> Some (max_num i j)
-
- let inter i1 i2 = 
-  match i1,i2 with
-   | Empty , _ -> Empty
-   |  _   , Empty -> Empty
-   |  Point n , Point m -> if n =/ m then i1 else Empty
-   | Point n , Itv (mn,mx) | Itv (mn,mx) , Point n-> 
-      if (match mn with 
-       | None -> true
-       | Some mn ->  mn <=/ n) && 
-       (match mx with 
-	| None -> true
-	| Some mx ->  n <=/ mx) then Point n else Empty
-   | Itv (min1,max1) , Itv (min2,max2) -> 
-      let bmin = max min1 min2
-      and bmax = min max1 max2 in
-       normalise bmin bmax
-
- (* a.x >= b*)
- let bound_of_constraint (a,b) =
-  match compare_num a (Int 0) with
-   | 0 -> 
-      if compare_num b (Int 0) = 1
-      then Empty  
-       (*actually this is a contradiction  failwith "bound_of_constraint" *)
-      else Itv (None,None)
-   | 1 -> Itv (Some (div_num  b a),None)
-   | -1 -> Itv (None, Some (div_num  b a))
-   | x -> failwith "bound_of_constraint(2)"
-
-
- let bounded x = 
-  match x with
-   | Itv(None,_) | Itv(_,None) -> false
-   |   _   -> true
-
-
- let range = function
-  | Empty -> Some (Int 0)
-  | Point n -> Some (Int (if ceiling_num n =/ n then 1 else 0))
-  | Itv(None,_) | Itv(_,None)-> None
-  | Itv(Some i,Some j) -> Some (floor_num j -/ceiling_num i +/ (Int 1))
-
- (* Returns the interval of smallest range *)
- let smaller_itv i1 i2 = 
-  match range i1 , range i2  with
-   | None , _ -> false
-   |  _   , None -> true
-   | Some i , Some j -> i <=/ j
+    let inter f o1 o2 = 
+      match o1 , o2 with
+	| None , None -> None
+	| Some _ , None -> o1
+	| None  , Some _ -> o2 
+	| Some n1 , Some n2 -> Some (f n1 n2) in
+
+    norm_itv (inter max_num l1 l2 , inter min_num r1 r2)
+
+  let range = function
+    | None,_ | _,None -> None
+    | Some i,Some j -> Some (floor_num j -/ceiling_num i +/ (Int 1))
+	
+
+  let smaller_itv i1 i2 = 
+    match  range i1 ,  range i2  with
+      | None , _ -> false
+      |  _   , None -> true
+      | Some i , Some j -> i <=/ j
+
+
+(** [in_bound bnd v] checks whether [v] is within the bounds [bnd] *)
+let in_bound bnd v =
+  let (l,r) = bnd in
+  match l , r with
+    | None , None -> true
+    | None , Some a -> v <=/ a
+    | Some a , None -> a <=/ v
+    | Some a , Some b -> a <=/ v && v <=/ b
 
 end
-open Interval
+open Itv 
+type vector = Vect.t
+
+type cstr = { coeffs : vector ; bound : interval }
+(** 'cstr' is the type of constraints.
+    {coeffs = v ; bound = (l,r) } models the constraints l <= v <= r
+**)
+
+module ISet = Set.Make(struct type t = int let compare = Pervasives.compare end)
+
+
+module PSet = ISet
+
+
+module System = Hashtbl.Make(Vect)
+
+  type proof = 
+      | Hyp of int 
+      | Elim of  var * proof * proof
+      | And of proof * proof
+
+
+
+type system = { 
+  sys : cstr_info ref System.t ; 
+  vars : ISet.t
+} 
+and cstr_info = { 
+  bound : interval ;
+  prf : proof ;
+  pos : int ;
+  neg : int ; 
+}
+
+
+(** A system of constraints has the form [{sys = s ; vars = v}].
+    [s] is a hashtable mapping a normalised vector to a [cstr_info] record where
+    - [bound] is an interval
+    - [prf_idx] is the set of hypothese indexes (i.e. constraints in the initial system) used to obtain the current constraint.
+       In the initial system, each constraint is given an unique singleton proof_idx.
+       When a new constraint c is computed by a function f(c1,...,cn), its proof_idx is ISet.fold union (List.map (fun x -> x.proof_idx) [c1;...;cn]
+    - [pos] is the number of positive values of the vector
+    - [neg] is the number of negative values of the vector
+    ( [neg] + [pos] is therefore the length of the vector) 
+    [v] is an upper-bound of the set of variables which appear in [s].
+*)
+
+(** To be thrown when a system has no solution *)
+exception SystemContradiction of proof
+let hyps prf = 
+  let rec hyps prf acc = 
+    match prf with
+      | Hyp i -> ISet.add i acc
+      | Elim(_,prf1,prf2) 
+      | And(prf1,prf2) -> hyps prf1 (hyps prf2 acc) in
+    hyps prf ISet.empty
+
+
+(** Pretty printing *)
+  let rec pp_proof o prf = 
+    match prf with
+      | Hyp i -> Printf.fprintf o "H%i" i
+      | Elim(v, prf1,prf2) -> Printf.fprintf o "E(%i,%a,%a)" v pp_proof prf1 pp_proof prf2
+      | And(prf1,prf2)   -> Printf.fprintf o "A(%a,%a)"  pp_proof prf1 pp_proof prf2
+	  
+let pp_bound o  = function
+  | None -> output_string o "oo"
+  | Some a -> output_string o (string_of_num a)
+
+let pp_itv o (l,r) = Printf.fprintf o "(%a,%a)" pp_bound l pp_bound r
+
+let rec pp_list f o l = 
+  match l with
+    | [] -> ()
+    | e::l -> f o e ; output_string o ";" ; pp_list f o l
 
-(* A set of constraints *)
-module Sys(V:Vector.S) (* : Vector.SystemS with module Vect = V*) = 
-struct
- 
- module Vect = V
-  
- module Cstr = Vector.Cstr(V)
- open Cstr
- 
- 
- module CMap = Map.Make(
-  struct
-   type t = Vect.t
-   let compare = Vect.compare
-  end)
-  
- module CstrBag = 
- struct 
+let pp_iset o s = 
+  output_string o "{" ;
+  ISet.fold (fun i _ -> Printf.fprintf o "%i " i) s ();
+  output_string o "}" 
 
-  type mut_itv = { mutable itv : intrvl} 
+let pp_pset o s = 
+  output_string o "{" ;
+  PSet.fold (fun i _ -> Printf.fprintf o "%i " i) s ();
+  output_string o "}" 
 
-  type t = mut_itv CMap.t
-    
-  exception Contradiction
-
-  let cstr_to_itv cstr = 
-   let (n,l) = V.normalise cstr.coeffs in 
-    if n =/ (Int 0)
-    then (Vect.null, bound_of_constraint (Int 0,cstr.cst)) (* Might be empty *)
-    else
-     match cstr.op with
-      | Eq -> let n = cstr.cst // n in (l, Point n)
-      | Ge -> 
-	 match  compare_num n (Int 0) with
-	  | 0 -> failwith "intrvl_of_constraint"
-	  | 1 -> (l,Itv (Some (cstr.cst // n), None))
-	  | -1 -> (l, Itv(None,Some (cstr.cst // n)))
-	  |  _ -> failwith "cstr_to_itv"
 
-	      
-  let empty = CMap.empty
+let pp_info o i = pp_itv o i.bound
 
+let pp_cstr o (vect,bnd) = 
+    let (l,r) = bnd in
+      (match l with
+	| None -> ()
+	| Some n -> Printf.fprintf o "%s <= " (string_of_num n))
+      ;
+      pp_vect o vect ; 
+      (match r with
+	      | None -> output_string o"\n"
+	      | Some n -> Printf.fprintf o "<=%s\n" (string_of_num n))
 
 
+let pp_system o sys= 
+  System.iter (fun vect ibnd -> 
+    pp_cstr o (vect,(!ibnd).bound)) sys
 
-  let is_empty = CMap.is_empty
 
-  let find_vect v bag = 
-   try
-    (bag,CMap.find v bag) 
-   with Not_found -> let x = { itv = Itv(None,None)} in (CMap.add v x bag ,x)
 
+let pp_split_cstr o (vl,v,c,_) = 
+  Printf.fprintf o "(val x = %s ,%a,%s)" (string_of_num vl) pp_vect  v  (string_of_num c)
 
-  let add (v,b) bag = 
-   match b with
-    | Empty -> raise Contradiction
-    | Itv(None,None) -> bag
-    | _    -> 
-       let (bag,intrl) = find_vect v bag in
-	match inter b intrl.itv with
-	 | Empty -> raise Contradiction
-	 |  itv  -> intrl.itv <- itv ; bag
+(** [merge_cstr_info] takes:
+    - the intersection of bounds and 
+    - the union of proofs
+    - [pos] and [neg] fields should be identical *)
 
-  exception Found of cstr
+let merge_cstr_info i1 i2 = 
+  let { pos = p1 ; neg = n1 ; bound = i1 ; prf = prf1 } = i1
+  and { pos = p2 ; neg = n2 ; bound = i2 ; prf = prf2 } = i2 in
+    assert (p1 = p2 && n1 = n2) ; 
+    match inter i1 i2 with
+      | None -> None (* Could directly raise a system contradiction exception *)
+      | Some bnd -> 
+	  Some { pos = p1 ; neg = n1 ; bound = bnd ; prf = And(prf1,prf2) }
 
-  let find_equation bag = 
-   try 
-    CMap.fold (fun v i () -> 
-     match i.itv with
-      | Point n -> let e = {coeffs = v ; op = Eq ; cst = n} 
-	in 	raise	(Found e)
-      |    _    -> () ) bag () ; None
-   with Found c -> Some c
+(** [xadd_cstr vect cstr_info] loads an constraint into the system.
+    The constraint is neither redundant nor contradictory.
+    @raise SystemContradiction if [cstr_info] returns [None]
+*)
 
+let xadd_cstr vect cstr_info sys =   
+  if debug && System.length sys mod 1000 = 0 then (print_string "*" ; flush stdout) ; 
+  try 
+    let info = System.find sys vect in
+      match merge_cstr_info cstr_info !info with
+	  | None       -> raise (SystemContradiction  (And(cstr_info.prf, (!info).prf)))
+	  | Some info' -> info := info'
+      with 
+	| Not_found -> System.replace  sys vect (ref cstr_info)
+
+
+type cstr_ext = 
+    | Contradiction (** The constraint is contradictory.
+			Typically, a [SystemContradiction] exception will be raised. *)
+    | Redundant  (** The constrain is redundant.
+		     Typically, the constraint will be dropped *)
+    | Cstr of vector * cstr_info (** Taken alone, the constraint is neither contradictory nor redundant.
+				     Typically, it will be added to the constraint system. *)
+
+(** [normalise_cstr] : vector -> cstr_info -> cstr_ext *)
+let normalise_cstr vect cinfo = 
+  match norm_itv cinfo.bound with
+    | None -> Contradiction
+    | Some (l,r) -> 
+	match vect with
+	  | [] -> if Itv.in_bound (l,r) (Int 0) then Redundant else Contradiction
+	  | (_,n)::_ ->  Cstr(
+	      (if n <>/ Int 1 then List.map (fun (x,nx) -> (x,nx // n)) vect else vect), 
+			     let divn x = x // n in
+			       if sign_num n = 1			       
+			       then{cinfo with bound = (map_option divn l , map_option  divn r) }
+			       else {cinfo with pos = cinfo.neg ; neg = cinfo.pos ; bound = (map_option divn r , map_option divn l)})
+
+(** For compatibility, there an external representation of constraints *)
+
+type cstr_compat = {coeffs : vector ; op : op ; cst : num}
+and op = |Eq | Ge
+
+let string_of_op = function Eq -> "=" | Ge -> ">="
+
+
+let eval_op = function
+  | Eq -> (=/)
+  | Ge ->  (>=/)
+
+let count v = 
+  let rec count n p v =
+    match v with
+      | [] -> (n,p)
+      | (_,vl)::v -> let sg = sign_num vl in 
+		       assert (sg <> 0) ; 
+	  if sg = 1 then count n (p+1) v else count (n+1) p v in
+    count 0 0 v
+
+
+let norm_cstr {coeffs = v ; op = o ; cst = c} idx =
+  let (n,p) = count v in 
+
+  normalise_cstr v {pos = p ; neg = n ; bound = 
+  (match o with 
+	| Eq -> Some c , Some c
+	| Ge -> Some c , None) ;
+	    prf = Hyp idx }
+
+
+(** [load_system l] takes a list of constraints of type [cstr_compat]
+    @return a system of constraints
+    @raise SystemContradiction if a contradiction is found
+*)
+let load_system l = 
     
-  let fold f bag acc = 
-   CMap.fold (fun v itv acc -> 
-    match itv.itv with
-     | Empty | Itv(None,None)  -> failwith "fold Empty"
-     | Itv(None ,Some i) -> 
-	f {coeffs = V.mul (Int (-1)) v ; op = Ge ; cst = minus_num i} acc
-     | Point n           ->   f {coeffs =  v ; op = Eq ; cst = n} acc
-     | Itv(x,y) -> 	
-	(match x with 
-	 | None -> (fun x -> x)
-	 | Some i -> f {coeffs = v ; op = Ge ; cst = i})
-	 (match y with
-	  | None -> acc
-	  | Some i -> 
-	     f {coeffs = V.mul (Int (-1)) v ; op = Ge ; cst = minus_num i} acc
-	 ) ) bag acc
-
-
-  let remove l _ = failwith "remove:Not implemented"
-   
-  module Map = 
-   Map.Make(
-    struct 
-     type t = int 
-     let compare : int -> int -> int = Pervasives.compare 
-    end)
-
-  let split f (t:t) =
-   let res = 
-    fold (fun e m -> let i = f e in
-                      Map.add i (add (cstr_to_itv e) 
-				  (try Map.find i m with 
-				    Not_found -> empty)) m) t Map.empty in
-    (fun i -> try Map.find i res with Not_found -> empty)
-     
-  type map = (int list * int list) Map.t
+  let sys = System.create 1000 in
     
+  let li = Mutils.mapi (fun e i -> (e,i)) l in
+
+  let vars = List.fold_left (fun vrs (cstr,i) -> 
+    match norm_cstr cstr i with
+      | Contradiction -> raise (SystemContradiction (Hyp i))
+      | Redundant      -> vrs
+      | Cstr(vect,info) -> 
+	  xadd_cstr  vect info sys ;
+	  List.fold_left (fun s (v,_) -> ISet.add v s) vrs cstr.coeffs) ISet.empty li   in
+
+    {sys = sys ;vars = vars}
+
+let system_list sys = 
+  let { sys = s ; vars  = v } = sys in  
+    System.fold (fun k bi l -> (k, !bi)::l) s [] 
+
+
+(** [add (v1,c1)  (v2,c2)  ] 
+    precondition:  (c1 <>/ Int 0 && c2 <>/ Int 0)
+    @return a pair [(v,ln)] such that 
+    [v] is the sum of vector [v1] divided by [c1] and vector [v2] divided by [c2]
+    Note that the resulting vector is not normalised.
+*)
+
+let add (v1,c1)  (v2,c2)  = 
+    assert (c1 <>/ Int 0 && c2 <>/ Int 0) ;
+
+  let rec xadd v1 v2  = 
+    match v1 , v2 with
+      | (x1,n1)::v1' , (x2,n2)::v2' -> 
+	  if x1 = x2 
+	  then 
+	    let n' = (n1 // c1) +/ (n2 // c2) in
+	      if n' =/ Int 0 then xadd v1' v2' 
+	      else 
+		let res = xadd v1' v2'  in
+		  (x1,n') ::res
+	  else if x1 < x2
+	  then let res = xadd v1' v2   in
+		 (x1, n1 // c1)::res 
+	  else let res = xadd v1 v2'  in
+		 (x2, n2 // c2)::res
+      |  [] , [] -> []
+      |  [] ,  _ -> List.map (fun (x,vl) -> (x,vl // c2)) v2
+      |  _  , [] -> List.map (fun (x,vl) -> (x,vl // c1)) v1 in
     
-  let status (b:t) = 
-   let _ , map = fold (fun c ( (idx:int),(res: map)) ->
-    ( idx + 1,
-    List.fold_left (fun (res:map) (pos,s) -> 
-     let (lp,ln) = try  Map.find pos res with Not_found -> ([],[]) in
-      match s with
-       | Vect.Pos -> Map.add pos (idx::lp,ln) res
-       | Vect.Neg -> 
-	  Map.add pos (lp, idx::ln) res) res 
-     (Vect.status c.coeffs))) b (0,Map.empty) in
-    Map.fold (fun k e res -> (k,e)::res)  map []
-     
-
-  type it = num CMap.t
+  let res = xadd v1 v2 in
+    (res, count res)
+
+let add (v1,c1)   (v2,c2)  = 
+  let res = add (v1,c1)   (v2,c2)  in
+    (*    Printf.printf "add(%a,%s,%a,%s) -> %a\n" pp_vect v1 (string_of_num c1) pp_vect v2 (string_of_num c2) pp_vect (fst res) ;*)
+    res
+
+type tlr = (num * vector * cstr_info) list
+type tm = (vector * cstr_info ) list
+
+(** To perform Fourier elimination, constraints are categorised depending on the sign of the variable to eliminate. *)
     
-  let iterator x = x
-   
-  let element it = failwith "element:Not implemented"
+(** [split x vect info (l,m,r)]
+    @param v is the variable to eliminate
+    @param l contains constraints such that (e + a*x) // a >= c / a 
+    @param r contains constraints such that (e + a*x) // - a >= c / -a
+    @param m contains constraints which do not mention [x]
+*)
+
+let split x (vect: vector) info (l,m,r) =
+    match get x vect with 
+      | None -> (* The constraint does not mention [x], store it in m *)
+	  (l,(vect,info)::m,r) 
+      | Some vl -> (* otherwise *)
+
+          let cons_bound lst bd = 
+            match  bd with
+              | None -> lst
+              | Some bnd -> (vl,vect,{info with bound = Some bnd,None})::lst in
+        
+          let lb,rb = info.bound in
+            if sign_num vl = 1 
+            then  (cons_bound l lb,m,cons_bound r rb)
+            else (* sign_num vl = -1 *)
+              (cons_bound l rb,m,cons_bound r lb)
+
+
+(** [project vr sys] projects system [sys] over the set of variables [ISet.remove vr sys.vars ].
+    This is a one step Fourier elimination.
+*)
+let project vr sys = 
+  
+  let (l,m,r)  = System.fold (fun vect rf l_m_r -> split vr vect !rf l_m_r) sys.sys  ([],[],[])  in
 
- end
-end
+  let new_sys = System.create (System.length sys.sys) in
+    
+    (* Constraints in [m] belong to the projection - for those [vr] is already projected out *)
+    List.iter (fun  (vect,info) -> System.replace new_sys vect (ref info) )  m ;
+
+    let elim (v1,vect1,info1) (v2,vect2,info2) = 
+      let {neg = n1 ; pos = p1 ; bound = bound1 ; prf = prf1} = info1
+      and {neg = n2 ; pos = p2 ; bound = bound2 ; prf = prf2} = info2 in
+
+      let bnd1 = from_option (fst bound1) 
+      and bnd2 = from_option (fst bound2) in
+      let bound = (bnd1 // v1) +/ (bnd2 // minus_num v2) in
+      let vres,(n,p) = add (vect1,v1) (vect2,minus_num v2)  in
+        (vres,{neg = n ; pos = p ; bound = (Some bound, None); prf = Elim(vr,info1.prf,info2.prf)}) in
+
+      List.iter(fun  l_elem -> List.iter (fun r_elem -> 
+        let (vect,info) = elim l_elem r_elem in
+	  match normalise_cstr vect info with
+	    | Redundant -> ()
+	    | Contradiction -> raise (SystemContradiction  info.prf)
+	    | Cstr(vect,info)   -> xadd_cstr vect info new_sys) r ) l;
+      {sys = new_sys ; vars = ISet.remove  vr sys.vars}
+      
 
-module Fourier(Vect : Vector.S)  =
-struct
- module Vect = Vect
- module Sys = Sys( Vect)
- module Cstr = Sys.Cstr
- module Bag = Sys.CstrBag
+(** [project_using_eq] performs elimination by pivoting using an equation.
+    This is the counter_part of the [elim] sub-function of [!project]. 
+    @param vr is the variable to be used as pivot
+    @param c is the coefficient of variable [vr] in vector [vect]
+    @param len is the length of the equation
+    @param bound is the bound of the equation
+    @param prf is the proof of the equation
+*)
+
+let project_using_eq vr c vect bound  prf (vect',info') = 
+    match get vr vect' with
+    | Some c2 -> 
+	let c1 = if c2 >=/ Int 0 then minus_num c else c in
+	  
+	let c2 = abs_num c2 in
+	  
+	let (vres,(n,p)) = add (vect,c1) (vect', c2)  in
+	  
+	let cst = bound // c1 in
+	  
+	let bndres = 
+	  let f x = cst +/ x // c2 in
+	  let (l,r) = info'.bound in
+	    (map_option f l , map_option f r) in
+	  
+	  (vres,{neg = n ; pos = p ; bound = bndres ; prf = Elim(vr,prf,info'.prf)})
+    | None -> (vect',info')
+
+let elim_var_using_eq vr vect cst  prf sys =
+  let c = from_option (get vr vect) in
+  
+  let elim_var = project_using_eq vr c vect cst prf  in
 
- open Cstr
- open Sys
+  let  new_sys  =  System.create (System.length sys.sys) in
 
- let debug = false
+    System.iter(fun vect iref -> 
+      let (vect',info') = elim_var (vect,!iref) in
+	match normalise_cstr vect' info' with
+	    | Redundant -> ()
+	    | Contradiction -> raise (SystemContradiction info'.prf)
+	    | Cstr(vect,info')   -> xadd_cstr vect info' new_sys) sys.sys ; 
+    
+    {sys = new_sys ; vars = ISet.remove  vr sys.vars}
 
- let print_bag msg b =
-  print_endline msg;
-  CstrBag.fold (fun e () -> print_endline (Cstr.string_of_cstr e)) b () 
+	
+(** [size sys] computes the number of entries in the system of constraints *)
+let size sys = System.fold (fun v iref s -> s + (!iref).neg + (!iref).pos) sys 0
+
+module IMap = Map.Make(struct type t = int let compare : int -> int -> int = Pervasives.compare end)
 
- let print_bag_file file msg b =
-  let f = open_out file in
-   output_string f msg;
-   CstrBag.fold (fun e () -> 
-    Printf.fprintf f "%s\n" (Cstr.string_of_cstr e)) b () 
+let pp_map o map = IMap.fold (fun k elt () -> Printf.fprintf o "%i -> %s\n" k (string_of_num elt)) map ()
 
+(** [eval_vect map vect] evaluates vector [vect] using the values of [map].
+    If [map] binds all the variables of [vect], we get
+    [eval_vect map [(x1,v1);...;(xn,vn)] = (IMap.find x1 map * v1) + ... + (IMap.find xn map) * vn , []]
+    The function returns as second argument, a sub-vector consisting in the variables that are not in [map]. *)
     
- (* A system with only inequations --
- *)
- let partition i m =
-  let splitter cstr =  compare_num (Vect.get i cstr.coeffs ) (Int 0) in
-  let split = CstrBag.split splitter m in
-   (split (-1) , split 0, split 1)
+let  eval_vect map vect = 
+  let rec xeval_vect vect sum rst = 
+    match vect with
+      | [] -> (sum,rst)
+      | (v,vl)::vect -> 
+	  try 
+	    let val_v = IMap.find v map in
+	      xeval_vect vect (sum +/ (val_v */ vl)) rst
+	  with
+	      Not_found -> xeval_vect vect sum ((v,vl)::rst) in
+    xeval_vect vect (Int 0) []
+  
 
+(** [restrict_bound n sum itv] returns the interval of [x]
+    given that (fst itv) <=  x * n + sum <= (snd itv) *)
+let restrict_bound n sum (itv:interval) = 
+  let f x  = (x -/ sum) // n in
+  let l,r = itv in
+    match sign_num n with
+      | 0 -> if in_bound itv sum
+	then (None,None) (* redundant *)
+	else failwith "SystemContradiction"
+      | 1 ->  map_option f l , map_option f r
+      | _ -> map_option f r , map_option f l
+
+
+(** [bound_of_variable map v sys] computes the interval of [v] in
+    [sys] given a mapping [map] binding all the other variables *)
+let bound_of_variable map v sys = 
+  System.fold (fun vect iref bnd  -> 
+    let sum,rst = eval_vect  map vect in
+    let vl = match get v rst with
+      | None -> Int 0
+      | Some v -> v in
+      match inter bnd (restrict_bound vl sum (!iref).bound) with
+	| None -> failwith "bound_of_variable: impossible"
+	| Some itv -> itv) sys  (None,None)
+
+
+(** [pick_small_value bnd] picks a value being closed to zero within the interval *)
+let pick_small_value bnd = 
+  match bnd with
+    | None , None   ->  Int 0
+    | None , Some i ->  if  (Int 0) <=/ (floor_num  i) then Int 0 else floor_num i
+    | Some i,None   ->  if i <=/ (Int 0) then Int 0 else ceiling_num i
+    | Some i,Some j -> 
+	if i <=/ Int 0 && Int 0 <=/ j 
+	then Int 0
+	else if ceiling_num i <=/ floor_num j 
+	then ceiling_num i (* why not *) else i
+
+
+(** [solution s1 sys_l  = Some(sn,[(vn-1,sn-1);...; (v1,s1)]@sys_l)] 
+    then [sn] is a  system which contains only [black_v] -- if it existed in [s1]
+    and [sn+1] is obtained by projecting [vn] out of [sn] 
+    @raise SystemContradiction if  system [s] has no solution 
+*)
+
+let solve_sys black_v choose_eq choose_variable sys sys_l = 
+
+  let rec solve_sys sys sys_l = 
+    if debug then Printf.printf "S #%i size %i\n" (System.length sys.sys) (size sys.sys);
     
- (* op of the result is arbitrary Ge *) 
- let lin_comb n1 c1 n2 c2 =
-  { coeffs = Vect.lin_comb n1 c1.coeffs n2 c2.coeffs ; 
-    op = Ge ; 
-    cst = (n1 */ c1.cst) +/ (n2 */ c2.cst)}
-
- (* BUG? : operator of the result ? *)
-
- let combine_project i c1 c2 =
-  let p = Vect.get  i c1.coeffs  
-  and n = Vect.get i c2.coeffs  in 
-   assert (n </ Int 0 && p >/ Int 0) ;
-   let nopp = minus_num n in 
-   let c =lin_comb nopp c1 p c2  in
-   let op = if c1.op = Ge ||  c2.op = Ge then Ge else Eq in
-    CstrBag.cstr_to_itv {coeffs =  c.coeffs ; op = op ; cst= c.cst }
-
-
- let project i m =
-  let (neg,zero,pos) = partition i m in
-  let project1 cpos acc = 
-   CstrBag.fold (fun cneg res ->
-    CstrBag.add (combine_project i cpos cneg) res) neg acc in
-   (CstrBag.fold project1 pos zero)
-
- (* Given a vector [x1 -> v1; ... ; xn -> vn]
-    and a constraint {x1 ; .... xn >= c }
- *)
- let evaluate_constraint i map cstr =
-  let {coeffs = _coeffs ; op = _op ; cst = _cst} = cstr in
-  let vi = Vect.get  i _coeffs  in
-  let v = Vect.set i (Int 0) _coeffs in
-   (vi, _cst -/ Vect.dotp map v)
-
-
- let rec bounds m itv =
-  match m with
-   | [] -> itv
-   | e::m -> bounds m (inter itv (bound_of_constraint e))
-
-
-
- let compare_status (i,(lp,ln)) (i',(lp',ln')) = 
-  let cmp = Pervasives.compare 
-   ((List.length lp) * (List.length ln))  
-   ((List.length lp') * (List.length ln')) in
-   if cmp = 0
-   then Pervasives.compare i i'
-   else cmp
-
- let cardinal m = CstrBag.fold (fun _ x -> x + 1) m 0
-
- let lightest_projection  l c m  =
-  let bound = c in 
-   if debug then (Printf.printf "l%i" bound; flush stdout) ; 
-   let rec xlight best l =
-    match l with
-     | [] -> best
-     | i::l -> 
-	let proj = (project i m) in
-	let cproj = cardinal proj in
-	 (*Printf.printf " p %i " cproj; flush stdout;*)
-	 match best with
-	  | None -> 
-	     if cproj < bound 
-	     then Some(cproj,proj,i) 
-	     else xlight (Some(cproj,proj,i)) l
-	  | Some (cbest,_,_) -> 		
-	     if cproj < cbest 
-	     then 
-	      if cproj < bound then Some(cproj,proj,i)
-	      else xlight (Some(cproj,proj,i)) l 
-	     else xlight best l in
-    match xlight None l with
-     | None -> None
-     | Some(_,p,i) -> Some (p,i)
-	
+    let eqs = choose_eq sys in
+      try 
+	let (v,vect,cst,ln) =  fst (List.find (fun ((v,_,_,_),_) -> v <> black_v) eqs) in
+	  if debug then 
+	    (Printf.printf "\nE %a = %s variable %i\n" pp_vect vect (string_of_num cst) v ;
+	     flush stdout);
+	  let sys' = elim_var_using_eq v vect cst ln sys in
+	    solve_sys sys' ((v,sys)::sys_l)  
+    with Not_found -> 
+      let vars = choose_variable  sys in
+	try 
+	  let (v,est) =  (List.find (fun (v,_) -> v <> black_v) vars) in
+	    if debug then (Printf.printf "\nV : %i esimate %f\n" v est ; flush stdout) ; 
+	    let sys' =  project v sys in
+	      solve_sys sys' ((v,sys)::sys_l) 
+	with Not_found ->  (* we are done *) Inl (sys,sys_l)  in
+    solve_sys sys sys_l
 
 
- exception Equality of cstr
 
- let find_equality m = Bag.find_equation m
 
-  
+let  solve black_v choose_eq choose_variable cstrs = 
 
- let pivot (n,v) eq ge =
-  assert (eq.op = Eq) ;
-  let res = 
-   match 
-    compare_num v (Int 0), 
-    compare_num (Vect.get n ge.coeffs) (Int 0) 
-   with
-    | 0 , _ -> failwith "Buggy"
-    | _ ,0  -> (CstrBag.cstr_to_itv ge)
-    | 1 , -1 -> combine_project n eq ge
-    | -1 , 1 -> combine_project n ge eq
-    | 1  , 1 -> 
-       combine_project n ge 
-	{coeffs = Vect.mul (Int (-1)) eq.coeffs; 
-	 op = eq.op ; 
-	 cst = minus_num eq.cst}
-    | -1 , -1 -> 
-       combine_project n 
-	{coeffs = Vect.mul (Int (-1)) eq.coeffs; 
-	 op = eq.op ; cst = minus_num eq.cst} ge
-    | _ -> failwith "pivot" in
-   res
-
- let check_cstr v c = 
-  let {coeffs = _coeffs ; op = _op ; cst = _cst} = c in
-  let vl = Vect.dotp v _coeffs in
-   match _op with
-    | Eq -> vl =/ _cst
-    | Ge -> vl >= _cst
-
-
- let forall p  sys = 
   try 
-   CstrBag.fold (fun c () -> if p c then () else raise Not_found) sys (); true
-  with Not_found -> false
-
-
- let check_sys v sys = forall (check_cstr v) sys
-
- let check_null_cstr c = 
-  let {coeffs = _coeffs ; op = _op ; cst = _cst} = c in
-   match _op with
-    | Eq -> (Int 0) =/ _cst
-    | Ge -> (Int 0) >= _cst
-       
- let check_null sys = forall check_null_cstr sys
-
-
- let optimise_ge 
-   quick_check choose choose_idx return_empty return_ge return_eq  m =
-  let c = cardinal m in
-  let bound =  2 * c  in
-   if debug then (Printf.printf "optimise_ge: %i\n" c; flush stdout);
-
-   let rec xoptimise  m =
-    if debug then (Printf.printf "x%i" (cardinal m) ; flush stdout);
-    if debug then (print_bag "xoptimise" m ; flush stdout);
-    if quick_check m
-    then return_empty m
-    else 
-     match find_equality m with
-      | None -> xoptimise_ge m
-      | Some eq -> xoptimise_eq eq m
-	 
-   and xoptimise_ge m =
-    begin
-     let c = cardinal m in
-     let l = List.map fst (List.sort compare_status (CstrBag.status m)) in
-     let  idx = choose bound l c m in
-      match idx with
-       | None -> return_empty m
-       | Some (proj,i) -> 
-	  match xoptimise proj with
-	   | None -> None
-	   | Some mapping -> return_ge m i mapping
-    end
-   and xoptimise_eq eq m = 
-    let l = List.map fst (Vect.status eq.coeffs) in
-     match choose_idx l with
-      | None -> (*if l = [] then None else*) return_empty m
-      | Some i -> 
-	 let p  = (i,Vect.get i eq.coeffs) in
-	 let m' = CstrBag.fold 
-	  (fun ge res -> CstrBag.add (pivot p eq ge) res) m CstrBag.empty in
-	  match xoptimise ( m') with
-	   | None -> None
-	   | Some mapp -> return_eq m eq i mapp in
-    try 
-     let res = xoptimise   m  in res
-    with CstrBag.Contradiction -> (*print_string "contradiction" ;*) None
-
-
-
- let minimise m = 
-  let opt_zero_choose bound l c m = 
-   if c > bound 
-   then lightest_projection l c m
-   else match l with
-    | [] -> None
-    | i::_ -> Some (project i m, i) in
-   
-  let choose_idx = function [] -> None | x::l -> Some x in
-
-  let opt_zero_return_empty m = Some Vect.null in
-
-   
-  let opt_zero_return_ge m i mapping = 
-   let (it:intrvl) = CstrBag.fold (fun cstr itv -> Interval.inter 
-    (bound_of_constraint (evaluate_constraint i mapping cstr)) itv) m 
-    (Itv (None, None)) in
-    match pick_closed_to_zero it with
-     | None -> print_endline "Cannot pick" ; None
-     | Some v -> 
-	let res =  (Vect.set i v mapping) in
-	 if debug 
-	 then Printf.printf "xoptimise res %i [%s]" i (Vect.string res) ;
-	 Some res in
-
-  let opt_zero_return_eq m eq i mapp = 
-   let (a,b) = evaluate_constraint i mapp eq in
-    Some (Vect.set i (div_num b a) mapp) in
-   
-   optimise_ge check_null opt_zero_choose 
-    choose_idx opt_zero_return_empty opt_zero_return_ge opt_zero_return_eq   m
-
- let normalise cstr = [CstrBag.cstr_to_itv cstr]
-
- let find_point l = 
-  (*  List.iter (fun e -> print_endline (Cstr.string_of_cstr e)) l;*)
-  try 
-   let m = List.fold_left (fun sys e -> CstrBag.add (CstrBag.cstr_to_itv e) sys) 
-    CstrBag.empty l in
-    match minimise m with
-     | None -> None
-     | Some res -> 
-	if debug then Printf.printf "[%s]"  (Vect.string res); 
-	Some res
-  with CstrBag.Contradiction -> None
+    let sys = load_system cstrs in
+(*      Printf.printf "solve :\n %a" pp_system sys.sys ; *)
+      solve_sys black_v choose_eq choose_variable sys []
+  with SystemContradiction prf -> Inr prf
+
+
+(** The  purpose of module [EstimateElimVar] is to try to estimate the cost of eliminating a variable.
+    The output is an ordered list of (variable,cost).
+*)
+
+module EstimateElimVar = 
+struct
+  type sys_list = (vector * cstr_info) list
+
+  let abstract_partition (v:int) (l: sys_list) =
+
+    let rec xpart (l:sys_list) (ltl:sys_list) (n:int list) (z:int) (p:int list) = 
+      match l with
+        | [] -> (ltl, n,z,p)
+        | (l1,info) ::rl -> 
+            match  l1 with
+            | [] -> xpart rl (([],info)::ltl) n (info.neg+info.pos+z) p
+            | (vr,vl)::rl1 -> 
+		if v = vr
+		then
+		  let cons_bound lst bd = 
+		    match  bd with
+		      | None -> lst
+		      | Some bnd -> info.neg+info.pos::lst in
+
+		  let lb,rb = info.bound in
+		    if sign_num vl = 1
+		    then  xpart rl ((rl1,info)::ltl) (cons_bound n lb) z (cons_bound p rb)
+		    else  xpart rl ((rl1,info)::ltl) (cons_bound n rb) z (cons_bound p lb)
+		else
+		  (* the variable is greater *)
+		  xpart rl ((l1,info)::ltl) n (info.neg+info.pos+z) p 
+
+    in
+    let (sys',n,z,p) =  xpart l [] [] 0 []  in
+
+    let ln = float_of_int (List.length n) in
+    let sn = float_of_int (List.fold_left (+) 0 n) in
+    let lp = float_of_int (List.length p) in
+    let sp = float_of_int (List.fold_left (+) 0 p) in
+      (sys',  float_of_int z +.  lp *.  sn +.  ln *.  sp -. lp*.ln)
+        
+        
+  let choose_variable   sys =
+    let {sys = s ; vars = v} = sys in
+      
+    let sl = system_list sys in
+
+    let evals  = fst
+      (ISet.fold (fun v (eval,s) -> let ts,vl = abstract_partition v s  in
+                                      ((v,vl)::eval, ts)) v ([],sl)) in
+      
+      List.sort (fun x y -> Pervasives.compare (snd x) (snd y) ) evals
 
 
- let find_q_interval_for x m = 
-  if debug then Printf.printf "find_q_interval_for %i\n" x ;
+end 
+open EstimateElimVar
+
+(** The  module [EstimateElimEq] is similar to [EstimateElimVar] but it orders equations.
+*)
+module EstimateElimEq =
+struct 
+  
+  let itv_point bnd = 
+    match bnd with
+      |(Some a, Some b) -> a =/ b
+      | _   -> false
+
+  let eq_bound bnd c = 
+    match bnd with
+      |(Some a, Some b) -> a =/ b && c =/ b
+      | _   -> false
+    
 
-  let choose bound l c m =
-   let rec xchoose l = 
+  let rec unroll_until v l = 
     match l with
-     | [] -> None
-     | i::l -> if i = x then xchoose l else Some (project i m,i) in
-    xchoose l in
-
-  let rec choose_idx = function 
-    [] -> None 
-   | e::l -> if e = x then choose_idx l else Some e in
-
-  let return_empty m = (* Beurk *)
-   (* returns the interval of x *)
-   Some (CstrBag.fold (fun cstr itv -> 
-    let i = if cstr.op = Eq 
-    then Point (cstr.cst // Vect.get x cstr.coeffs)  
-     else if Vect.is_null (Vect.set x (Int 0) cstr.coeffs) 
-     then  bound_of_constraint (Vect.get x cstr.coeffs  , cstr.cst) 
-     else itv
-   in
-		       Interval.inter i itv) m (Itv (None, None)))  in
-   
-  let return_ge m i res = Some res in
-
-  let return_eq m eq i res = Some res in
-
-   try 
-    optimise_ge 
-     (fun x -> false) choose choose_idx return_empty return_ge return_eq m
-   with CstrBag.Contradiction -> None
-
-
- let find_q_intervals sys =
-  let variables =  
-   List.map fst (List.sort compare_status (CstrBag.status sys)) in
-   List.map (fun x -> (x,find_q_interval_for x sys)) variables
-
- let pp_option f o = function 
-   None -> Printf.fprintf o "None" 
-  | Some x -> Printf.fprintf o "Some %a" f x
-
- let optimise vect sys = 
-  (* we have to modify the system with a dummy variable *)
-  let fresh = 
-   List.fold_left (fun fr c -> Pervasives.max fr (Vect.fresh c.coeffs)) 0 sys in
-   assert (List.for_all (fun x -> Vect.get fresh x.coeffs =/ Int 0) sys);
-   let cstr = {
-    coeffs = Vect.set fresh (Int (-1)) vect ; 
-    op = Eq ; 
-    cst = (Int 0)} in
-    try 
-     find_q_interval_for fresh 
-      (List.fold_left 
-	(fun bg c -> CstrBag.add (CstrBag.cstr_to_itv c) bg)  
-	CstrBag.empty (cstr::sys))
-    with CstrBag.Contradiction -> None
-
-
- let optimise vect sys = 
-  let res = optimise vect sys in 
-   if debug 
-   then Printf.printf "optimise %s -> %a\n" 
-    (Vect.string vect) (pp_option (fun o x -> Printf.printf "%s" (string_of_intrvl x))) res  
-   ; res
-
- let find_Q_interval sys = 
-  try 
-   let sys = 
-    (List.fold_left 
-      (fun bg c -> CstrBag.add (CstrBag.cstr_to_itv c) bg)  CstrBag.empty sys) in
-   let candidates = 
-    List.fold_left 
-     (fun l (x,i) -> match i with 
-       None -> (x,Empty)::l 
-      | Some i ->  (x,i)::l) [] (find_q_intervals sys) in
-    match List.fold_left 
-     (fun (x1,i1) (x2,i2) -> 
-      if smaller_itv i1 i2 
-      then (x1,i1) else (x2,i2)) (-1,Itv(None,None)) candidates 
-    with
-     | (i,Empty) -> None
-     | (x,Itv(Some i, Some j))  -> Some(i,x,j)
-     | (x,Point n) -> Some(n,x,n)
-     |   _        -> None
-  with CstrBag.Contradiction -> None
+      | [] -> (false,[])
+      | (i,_)::rl -> if i = v 
+	then (true,rl) 
+	else if i < v then unroll_until v rl else (false,l)
+
+
+  let choose_primal_equation eqs sys_l = 
+
+    let is_primal_equation_var v = 
+      List.fold_left (fun (nb_eq,nb_cst) (vect,info) -> 
+	if fst (unroll_until v vect) 
+	then if itv_point  info.bound then (nb_eq +  1,nb_cst) else (nb_eq,nb_cst)
+	else (nb_eq,nb_cst)) (0,0) sys_l in
+
+    let rec find_var vect = 
+      match vect with
+	| [] -> None
+	| (i,_)::vect -> 
+	    let (nb_eq,nb_cst) = is_primal_equation_var i in
+	      if nb_eq = 2 && nb_cst = 0
+	      then Some i else find_var vect in
+
+    let rec find_eq_var eqs = 
+      match eqs with
+	| [] -> None
+	| (vect,a,prf,ln)::l -> 
+	    match find_var vect with 
+		| None -> find_eq_var l
+		| Some r -> Some (r,vect,a,prf,ln) 
+    in
+      
+
+      find_eq_var eqs
+
+
+
+
+  let choose_equality_var  sys =
+
+    let sys_l = system_list sys in
+
+    let equalities = List.fold_left 
+      (fun  l (vect,info) -> 
+	match  info.bound with
+	  | Some a , Some b -> 
+	      if a =/ b then (* This an equation *)
+		(vect,a,info.prf,info.neg+info.pos)::l else l
+	  |   _ -> l
+      ) [] sys_l  in
+
+    let rec estimate_cost v ct sysl acc tlsys = 
+      match sysl with
+	| [] -> (acc,tlsys)
+	| (l,info)::rsys ->
+	    let ln = info.pos + info.neg in
+	    let (b,l) = unroll_until v l in
+	    match b with
+	      | true -> 
+		  if itv_point info.bound 
+		  then estimate_cost  v ct rsys (acc+ln) ((l,info)::tlsys) (* this is free *)
+		  else estimate_cost v ct rsys (acc+ln+ct) ((l,info)::tlsys)  (* should be more ? *)
+	      | false -> estimate_cost v ct rsys (acc+ln) ((l,info)::tlsys) in
+
+      match choose_primal_equation equalities sys_l with
+	| None -> 
+	    let cost_eq eq const prf ln acc_costs = 
+	      
+	      let rec cost_eq eqr sysl costs = 
+		match eqr with
+		  | [] -> costs
+		  | (v,_) ::eqr -> let (cst,tlsys) = estimate_cost v (ln-1) sysl 0 [] in
+				     cost_eq eqr tlsys (((v,eq,const,prf),cst)::costs) in
+		cost_eq eq sys_l acc_costs     in
+
+	    let all_costs = List.fold_left (fun all_costs (vect,const,prf,ln) -> cost_eq vect const prf ln all_costs) [] equalities in
+
+	      (*      pp_list (fun o ((v,eq,_,_),cst) -> Printf.fprintf o "((%i,%a),%i)\n" v pp_vect eq cst) stdout all_costs ; *)
+	      
+	      List.sort (fun x y -> Pervasives.compare (snd x) (snd y) ) all_costs
+	| Some (v,vect, const,prf,_) -> [(v,vect,const,prf),0]
 
 
 end
+open EstimateElimEq
+
+
+
+module Fourier =
+struct
+
+  let optimise vect l = 
+    (* We add a dummy (fresh) variable for vector *)
+    let fresh = 
+      List.fold_left (fun fr c -> Pervasives.max fr (Vect.fresh c.coeffs)) 0 l in
+    let cstr = {
+      coeffs = Vect.set fresh (Int (-1)) vect ; 
+      op = Eq ; 
+      cst = (Int 0)} in
+      match solve fresh choose_equality_var choose_variable (cstr::l) with
+	| Inr prf -> None (* This is an unsatisfiability proof *)
+	| Inl (s,_) -> 
+	    try 
+	      Some (bound_of_variable IMap.empty fresh s.sys)
+	    with
+		x -> Printf.printf "optimise Exception : %s" (Printexc.to_string x) ; None
+
+
+  let find_point cstrs = 
+  
+    match solve max_int choose_equality_var choose_variable cstrs with
+      | Inr prf -> Inr prf
+      | Inl (_,l) -> 
+      
+	let rec rebuild_solution l map = 
+	  match l with
+	    | [] -> map
+	    | (v,e)::l -> 
+		let itv = bound_of_variable map v e.sys in
+		let map = IMap.add v (pick_small_value itv) map in
+		  rebuild_solution l map
+	in
+
+	let map = rebuild_solution l IMap.empty in
+	  let vect = List.rev (IMap.fold (fun v i vect -> (v,i)::vect) map []) in
+(*	    Printf.printf "SOLUTION %a" pp_vect vect ; *)
+	  let res = Inl vect in
+	    res
+
+
+end
+
+
+module Proof =
+struct 
+  
+	
+
+
+(** A proof term in the sense of a ZMicromega.RatProof is a positive combination of the hypotheses which leads to a contradiction.
+    The proofs constructed by Fourier elimination are more like execution traces:
+         - certain facts are recorded but are useless
+         - certain inferences are implicit.
+    The following code implements proof reconstruction.
+*)
+  let add x y = fst (add x y)
+
+
+  let forall_pairs f l1 l2 =
+    List.fold_left (fun acc e1 ->
+      List.fold_left (fun acc e2 -> 
+	match f e1 e2 with
+	  | None -> acc
+	  | Some v -> v::acc) acc l2) [] l1
+
+
+  let add_op x y = 
+    match x , y with
+      | Eq , Eq -> Eq
+      |     _   -> Ge
+
+
+  let pivot v (p1,c1) (p2,c2) = 
+    let {coeffs = v1 ; op = op1 ; cst = n1} = c1
+    and {coeffs = v2 ; op = op2 ; cst = n2} = c2 in
+      
+      match Vect.get v v1 , Vect.get v v2 with
+        | None , _ | _ , None -> None
+        | Some a   , Some b   -> 
+            if (sign_num a) * (sign_num b) = -1
+            then Some (add (p1,abs_num a) (p2,abs_num b) , 
+                      {coeffs = add (v1,abs_num a) (v2,abs_num b) ; 
+                       op = add_op op1 op2 ;
+                       cst = n1 // (abs_num a) +/ n2 // (abs_num b) })
+            else if op1 = Eq
+            then Some (add (p1,minus_num (a // b)) (p2,Int 1), 
+                      {coeffs = add (v1,minus_num (a// b)) (v2 ,Int 1) ; 
+                       op     = add_op op1 op2;
+                       cst    = n1 // (minus_num (a// b)) +/ n2 // (Int 1)})
+            else if op2 = Eq
+	    then
+	      Some (add (p2,minus_num (b // a)) (p1,Int 1), 
+                   {coeffs = add (v2,minus_num (b// a)) (v1 ,Int 1) ; 
+                    op     = add_op op1 op2;
+                    cst    = n2 // (minus_num (b// a)) +/ n1 // (Int 1)})
+	    else  None (* op2 could be Eq ... this might happen *) 
+    
+
+  let normalise_proofs l = 
+    List.fold_left (fun acc (prf,cstr) -> 
+      match acc with
+        | Inr _ -> acc (* I already found a contradiction *)
+        | Inl acc -> 
+            match norm_cstr cstr 0 with
+              | Redundant -> Inl acc
+              | Contradiction -> Inr (prf,cstr)
+              | Cstr(v,info)  -> Inl ((prf,cstr,v,info)::acc)) (Inl []) l
+
+
+  type oproof = (vector * cstr_compat * num) option
+
+  let merge_proof (oleft:oproof) (prf,cstr,v,info) (oright:oproof) = 
+    let (l,r) = info.bound in
+
+    let keep p ob bd = 
+      match ob , bd with 
+        | None , None -> None
+        | None , Some b -> Some(prf,cstr,b)
+        | Some _ , None -> ob
+        | Some(prfl,cstrl,bl) , Some b -> if p bl b then Some(prf,cstr, b) else ob in
+
+    let oleft  = keep (<=/) oleft l in
+    let oright = keep (>=/) oright r in
+      (* Now, there might be a contradiction *)
+      match oleft , oright with
+        | None , _ | _ , None -> Inl (oleft,oright)
+        | Some(prfl,cstrl,l) , Some(prfr,cstrr,r) -> 
+            if l <=/ r 
+            then Inl (oleft,oright)
+            else (* There is a contradiction - it should show up by scaling up the vectors - any pivot should do*)
+              match cstrr.coeffs with
+                | [] -> Inr (add (prfl,Int 1) (prfr,Int 1), cstrr) (* this is wrong *)
+                | (v,_)::_ -> 
+                    match pivot v (prfl,cstrl) (prfr,cstrr) with
+                      | None -> failwith "merge_proof : pivot is not possible"
+                      | Some x -> Inr x
+
+let  mk_proof hyps prf = 
+  (* I am keeping list - I might have a proof for the left bound and a proof for the right bound.
+     If I perform aggressive elimination of redundancies, I expect the list to be of length at most 2.
+     For each proof list, all the vectors should be of the form a.v for different constants a.
+  *)
+
+  let rec mk_proof prf = 
+    match prf with
+      | Hyp i -> [ ([i, Int 1] , List.nth hyps i) ]
+
+      | Elim(v,prf1,prf2) ->
+          let prfsl = mk_proof prf1
+          and prfsr = mk_proof prf2 in
+            (* I take only the pairs for which the elimination is meaningfull *)
+            forall_pairs (pivot v) prfsl prfsr
+      | And(prf1,prf2) -> 
+          let prfsl1 = mk_proof prf1 
+          and prfsl2 = mk_proof prf2 in
+          (* detect trivial redundancies and contradictions *)
+            match normalise_proofs (prfsl1@prfsl2) with
+              | Inr x -> [x] (* This is a contradiction - this should be the end of the proof *)
+              | Inl l -> (* All the vectors are the same *)
+                  let prfs = 
+                    List.fold_left (fun acc e -> 
+                      match acc with
+                        | Inr _ -> acc (* I have a contradiction *)
+                        | Inl (oleft,oright) -> merge_proof oleft e oright) (Inl(None,None)) l in
+                    match prfs with
+                      | Inr x -> [x]
+                      | Inl (oleft,oright) ->
+			  match oleft , oright with
+			    | None , None -> []
+			    | None , Some(prf,cstr,_) | Some(prf,cstr,_) , None -> [prf,cstr]
+			    | Some(prf1,cstr1,_) , Some(prf2,cstr2,_) -> [prf1,cstr1;prf2,cstr2] in
+
+    mk_proof prf
+
+
+end 
 
diff --git a/contrib/micromega/mutils.ml b/contrib/micromega/mutils.ml
index 2473608f2..03e3999fd 100644
--- a/contrib/micromega/mutils.ml
+++ b/contrib/micromega/mutils.ml
@@ -17,6 +17,15 @@ let debug = false
 let fst' (Micromega.Pair(x,y)) = x
 let snd' (Micromega.Pair(x,y)) = y
 
+let map_option f x = 
+  match x with
+    | None -> None
+    | Some v -> Some (f v)
+
+let from_option = function
+  | None -> failwith "from_option"
+  | Some v -> v  
+
 let rec try_any l x = 
  match l with
   | [] -> None
@@ -24,6 +33,37 @@ let rec try_any l x =
      | None -> try_any l x
      | x -> x
 
+let iteri f l =
+ let rec xiter i l = 
+  match l with
+   | [] -> ()
+   | e::l -> f i e ; xiter (i+1) l in
+  xiter 0 l
+
+let mapi f l =
+ let rec xmap i l = 
+  match l with
+   | [] -> []
+   | e::l -> (f i e)::xmap (i+1) l in
+  xmap 0 l
+
+let rec map3 f l1 l2 l3 = 
+  match l1 , l2 ,l3 with
+    | [] , [] , [] -> []
+    | e1::l1 , e2::l2 , e3::l3 -> (f e1 e2 e3)::(map3 f l1 l2 l3)
+    |      _   -> raise (Invalid_argument "map3")
+
+
+let rec is_sublist l1 l2 = 
+  match l1 ,l2 with
+    | [] ,_ -> true
+    | e::l1', [] -> false
+    | e::l1' , e'::l2' -> 
+	if e = e' then is_sublist l1' l2'
+	else is_sublist l1 l2'
+		
+
+
 let list_try_find f =
   let rec try_find_f = function
     | [] -> failwith "try_find"
diff --git a/contrib/micromega/vector.ml b/contrib/micromega/vector.ml
deleted file mode 100644
index fee4ebfc1..000000000
--- a/contrib/micromega/vector.ml
+++ /dev/null
@@ -1,674 +0,0 @@
-(************************************************************************)
-(*  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        *)
-(************************************************************************)
-(*                                                                      *)
-(* Micromega: A reflexive tactic using the Positivstellensatz           *)
-(*                                                                      *)
-(*  Frédéric Besson (Irisa/Inria) 2006-2008                             *)
-(*                                                                      *)
-(************************************************************************)
-
-open Num
-
-module type S =
-sig
- type t 
-
- val fresh : t -> int
-
- val null : t
-
- val is_null : t -> bool
-
- val get : int -> t -> num
-
- val update : int -> (num -> num) -> t -> t
-  (* behaviour is undef if index < 0 -- might loop*)
-
- val set : int -> num -> t -> t
-
- (* 
-    For efficiency...
-
-    val get_update : int -> (num -> num) -> t -> num * t  
- *)
-
- val mul : num -> t -> t
-
- val uminus : t -> t
-
- val add : t -> t -> t
-
- val dotp : t -> t -> num
-
- val lin_comb : num -> t -> num -> t -> t
-  (* lin_comb n1 t1 n2 t2 = (n1 * t1) + (n2 * t2) *)
-
- val gcd : t -> Big_int.big_int
-
- val normalise : t -> num * t 
-
- val hash : t -> int
-
- val compare : t -> t -> int 
-
- type it 
-
- val iterator : t -> it
- val element : it -> (num*it) option
-
- val string : t -> string
-
- type status = Pos | Neg
-
- (* the result list is ordered by fst *)
- val status : t -> (int * status) list 
-
- val from_list : num list -> t
- val to_list : t -> num list
-
-end
-
-
-module type SystemS = 
-sig
- 
- module Vect : S
-
- module  Cstr :
- sig
-  type kind = Eq | Ge
-  val string_of_kind : kind -> string
-  type cstr = {coeffs : Vect.t ; op : kind ; cst : num} 
-  val  string_of_cstr : cstr -> string
-  val compare : cstr -> cstr -> int
- end
- open Cstr
-
-
- module CstrBag : 
- sig 
-  type t
-  exception Contradiction
-
-  val empty : t
-
-  val is_empty : t -> bool
-
-  val add : cstr -> t -> t
-   (* c can be deduced from add c t *)
-
-  val find : (cstr -> bool) -> t  ->  cstr option
-
-  val fold : (cstr -> 'a -> 'a) -> t -> 'a -> 'a
-
-  val status : t ->  (int * (int list * int list)) list
-   (* aggregate of vector statuses *)
-
-  val remove : cstr -> t -> t
-
-  (* remove_list the ith element -- it is the ith element visited by 'fold' *)
-
-  val split : (cstr -> int) -> t -> (int ->  t)
-
-  type it 
-  val iterator : t -> it
-  val element : it -> (cstr*it) option
-
- end
-
-end
-
-let zero_num = Int 0
-let unit_num = Int 1
-
-
-
-
-module Cstr(V:S) = 
-struct 
- type kind = Eq | Ge
- let string_of_kind = function Eq -> "Eq" | Ge -> "Ge"
-
- type cstr = {coeffs : V.t ; op : kind ; cst : num} 
-
- let string_of_cstr  {coeffs =a  ; op = b ; cst =c} =
-  Printf.sprintf "{coeffs = %s;op=%s;cst=%s}" (V.string a) (string_of_kind b) (string_of_num c)
-
- type t = cstr
- let compare 
-   {coeffs = v1 ; op = op1 ; cst = c1}
-   {coeffs = v2 ; op = op2 ; cst = c2} =
-  Mutils.Cmp.compare_lexical [
-   (fun () -> V.compare v1 v2);
-   (fun () -> Pervasives.compare op1 op2);
-   (fun () -> compare_num c1 c2)
-  ]
-
-
-end
-
-
-
-module VList : S with type t = num list =
-struct
- type t =  num list
-
- let fresh l = failwith "not implemented"
-
- let null = []
-
- let is_null  = List.for_all ((=/) zero_num) 
-
- let normalise l = failwith "Not implemented"
-  (*match l with (* Buggy : What if the first num is zero! *)
-    | [] -> (Int 0,[])
-    | [n] -> (n,[Int 1])
-    | n::l -> (n, (Int 1) :: List.map (fun x -> x // n) l)
-  *)
-
-
- let get i l = try List.nth l i with _ -> zero_num
-
- (* This is not tail-recursive *)
- let rec update i f t =
-  match t with
-   | []   -> if i = 0 then [f zero_num] else (zero_num)::(update (i-1) f [])
-   | e::t -> if i = 0 then (f e)::t else e::(update (i-1) f t)
-
- let rec set i n t =
-  match t with
-   | []   -> if i = 0 then [n] else (zero_num)::(set (i-1) n [])
-   | e::t -> if i = 0 then (n)::t else e::(set (i-1) n t)
-
-
-
-
- let rec mul z t = 
-  match z with
-   | Int 0 -> null
-   | Int 1 -> t
-   |  _    -> List.map (mult_num z) t
-
- let uminus t = mul (Int (-1)) t
-
- let rec add t1 t2 =
-  match t1,t2 with
-   | [], _  -> t2
-   | _ , [] -> t1
-   | e1::t1,e2::t2 -> (e1 +/ e2 )::(add t1 t2)
-
- let dotp t1 t2 =
-  let rec _dotp t1 t2 acc = 
-   match t1, t2 with
-    | [] , _  -> acc
-    | _  , [] -> acc
-    | e1::t1,e2::t2 -> _dotp t1 t2 (acc +/ (e1 */ e2)) in
-   _dotp t1 t2 zero_num
-
- let add_mul n t1 t2 =
-  match n with
-   | Int 0 -> t2
-   | Int 1  -> add t1 t2
-   |  _     -> 
-       let rec _add_mul t1 t2 = 
-	match t1,t2 with
-	 | [], _  -> t2
-	 | _ , [] -> mul n t1
-	 | e1::t1,e2::t2 -> ( (n */e1) +/ e2 )::(_add_mul t1 t2) in
-	_add_mul t1 t2
-
- let lin_comb n1 t1 n2 t2 = 
-  match n1,n2 with
-   | Int 0 , _  -> mul n2 t2
-   | Int 1 , _ -> add_mul n2 t2 t1
-   |   _   , Int 0 ->  mul n1 t1
-   |   _   , Int 1 -> add_mul n1 t1 t2  
-   |       _       -> 
-	    let rec _lin_comb t1 t2 = 
-	     match t1,t2 with
-	      | [], _  -> mul n2 t2
-	      | _ , [] -> mul n1 t1
-	      | e1::t1,e2::t2 -> ( (n1 */e1) +/ (n2 */ e2 ))::(_lin_comb t1 t2) in
-	     _lin_comb t1 t2
-	      
- (* could be computed on the fly *)
- let gcd  t =Mutils.gcd_list t
-
-
-
-
- let hash = Mutils.Cmp.hash_list int_of_num 
-
- let compare = Mutils.Cmp.compare_list compare_num
-
- type it = t
- let iterator (x:t) : it = x
- let element it = 
-  match it with
-   | [] -> None
-   | e::l -> Some (e,l)
-
- (* TODO: Buffer! *)
- let string l = List.fold_right (fun n s -> (string_of_num n)^";"^s) l ""
-
- type status = Pos | Neg
-
- let status l =
-  let rec xstatus i l = 
-   match l with 
-    | [] -> []
-    | e::l -> 
-       begin
-	match compare_num e (Int 0) with
-	 | 1 -> (i,Pos):: (xstatus (i+1) l)
-	 | 0 -> xstatus (i+1) l
-	 | -1 -> (i,Neg) :: (xstatus (i+1) l)
-	 | _  -> assert false
-       end in
-   xstatus 0 l
-
- let from_list l = l
- let to_list l = l
-  
-end
-
-module VMap : S =
-struct
- module Map = Map.Make(struct type t = int let compare (x:int) (y:int) = Pervasives.compare x y end)
-
- type t =  num Map.t
-
- let null = Map.empty
-
- let fresh m = failwith "not implemented"
-
- let is_null  = Map.is_empty
-
- let normalise m = failwith "Not implemented"
-
-
-
- let get i l = try Map.find i l with _ -> zero_num
-
- let  update i f t =
-  try
-   let res = f (Map.find i t) in
-    if res =/ zero_num
-    then Map.remove i t
-    else Map.add i res t
-  with
-    Not_found -> 
-     let res = f zero_num in
-      if res =/ zero_num then t else Map.add i res t
-
- let set i n t =
-  if n =/ zero_num then Map.remove i t
-  else Map.add i n t
-
-
- let rec mul z t = 
-  match z with
-   | Int 0 -> null
-   | Int 1 -> t
-   |  _    -> Map.map (mult_num z) t
-
- let uminus t = mul (Int (-1)) t
-  
-
- let map2 f m1 m2 = 
-  let res,m2' = 
-   Map.fold (fun k e (res,m2) -> 
-    let v = f e (get k m2) in
-     if v =/ zero_num
-     then (res,Map.remove k m2)
-     else (Map.add k v res,Map.remove k m2)) m1 (Map.empty,m2) in
-   Map.fold (fun k e res -> 
-    let v = f zero_num e in
-     if v =/ zero_num
-     then res else Map.add k v res) m2' res
-    
- let  add t1 t2 = map2 (+/) t1 t2
-  
-
- let dotp t1 t2 =
-  Map.fold (fun k e res -> 
-   res +/ (e */ get k t2)) t1 zero_num
-
-
-
- let add_mul n t1 t2 =
-  match n with
-   | Int 0 -> t2
-   | Int 1  -> add t1 t2
-   |  _     -> map2 (fun x y -> (n */ x) +/ y) t1 t2
-
- let lin_comb n1 t1 n2 t2 = 
-  match n1,n2 with
-   | Int 0 , _  -> mul n2 t2
-   | Int 1 , _ -> add_mul n2 t2 t1
-   |   _   , Int 0 ->  mul n1 t1
-   |   _   , Int 1 -> add_mul n1 t1 t2  
-   |       _       -> map2 (fun x y -> (n1 */ x) +/ (n2 */ y)) t1 t2
-
-
- let hash map = Map.fold (fun k e res -> k lxor (int_of_num e) lxor res) map 0
-
- let compare = Map.compare compare_num
-
- type it = t * int
-
- let iterator (x:t) : it = (x,0)
-  
- let element (mp,id) = 
-  try 
-   Some (Map.find id mp, (mp, id+1))
-  with
-    Not_found -> None
-
- (* TODO: Buffer! *)
- type status = Pos | Neg
-
- let status l = Map.fold (fun k e l -> 
-  match compare_num e (Int 0) with
-   | 1 -> (k,Pos)::l
-   | 0 ->  l
-   | -1 -> (k,Neg) :: l
-   | _  -> assert false) l []
- let from_list l = 
-  let rec from_list i l map = 
-   match l with
-    | [] -> map
-    | e::l -> from_list (i+1) l (if e <>/ Int 0 then Map.add i e map else map) in
-   from_list 0 l Map.empty
-
- let gcd m = 
-  let res = Map.fold (fun _ e x -> Big_int.gcd_big_int  x (Mutils.numerator e)) m Big_int.zero_big_int in
-   if Big_int.compare_big_int res Big_int.zero_big_int = 0 
-   then Big_int.unit_big_int else res
-
-
- let to_list m = 
-  let l = List.rev (Map.fold (fun k e l -> (k,e)::l) m []) in
-  let rec xto_list i l =
-   match l with
-    | [] -> []
-    |	(x,v)::l' -> if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in
-   xto_list 0 l
-
- let string l = VList.string (to_list l)
-
-  
-end
-
-
-module VSparse : S =
-struct
-
- type t = (int*num) list
-
- let null = []
-  
- let fresh l = List.fold_left (fun acc (i,_) -> max (i+1) acc)  0 l
-
- let is_null l = l = []
-  
- let rec is_sorted l =
-  match l with
-   |	[] -> true
-   | [e] -> true
-   | (i,_)::(j,x)::l -> i < j && is_sorted ((j,x)::l)
-      
-      
- let check l = (List.for_all (fun (_,n) -> compare_num n (Int 0) <> 0) l) && (is_sorted l)
-  
- (*	let get i t = 
-	assert (check t);
-	try List.assoc i t with Not_found -> zero_num *)
-
- let rec get (i:int) t = 
-  match t with
-   | [] -> zero_num
-   | (j,n)::t -> 
-      match compare i j with
-       | 0 -> n
-       | 1 -> get i t
-       |  _ -> zero_num
-
- let cons i v rst = if v =/ Int 0 then rst else (i,v)::rst
-   
- let rec update i f t = 
-  match t with
-   | [] -> cons i (f zero_num) []
-   | (k,v)::l -> 
-      match Pervasives.compare i k with
-       | 0 -> cons k (f v) l
-       | -1 -> cons i (f zero_num) t
-       |  1 -> (k,v) ::(update i f l)
-       |  _  -> failwith "compare_num"
-	   
- let update i f t =
-  assert (check t); 
-  let res = update i f t in
-   assert (check t) ; res
-    
-    
- let rec set i n t =
-  match t with
-   | [] -> cons i n []
-   | (k,v)::l -> 
-      match Pervasives.compare i k with
-       | 0 -> cons k n l
-       | -1 -> cons i n t
-       |  1 -> (k,v) :: (set i n l)
-       |  _  -> failwith "compare_num"
-	   
-
- let rec map f l = 
-  match l with
-   | [] -> []
-   | (i,e)::l -> cons i (f e) (map f l)
-      
- let rec mul z t = 
-  match z with
-   | Int 0 -> null
-   | Int 1 -> t
-   |  _    -> List.map (fun (i,n) -> (i, mult_num z n)) t
-       
- let mul z t = 
-  assert (check t) ;
-  let res = mul z t in
-   assert (check res) ; 
-   res
-
- let uminus t = mul (Int (-1)) t
-
-
- let normalise l = 
-  match l with
-   | []   -> (Int 0,[])
-   | (i,n)::_ -> (n, mul ((Int 1) // n) l)
-      
-
- let rec map2 f m1 m2 = 
-  match m1, m2 with
-   | [] , [] -> []
-   | l  , [] -> map (fun x -> f x zero_num) l
-   | [] ,l   -> map (f zero_num) l
-   | (i,e)::l1,(i',e')::l2 -> 
-      match Pervasives.compare i i' with
-       | 0 -> cons i (f e e') (map2 f l1 l2)
-       | -1 -> cons i (f e zero_num) (map2 f l1 m2)
-       |  1 -> cons i' (f zero_num e') (map2 f m1 l2)
-       | _  -> assert false
-	  
- (* let  add t1 t2 = map2 (+/) t1 t2*)
-
- let rec add (m1:t) (m2:t) = 
-  match m1, m2 with
-   | [] , [] -> []
-   | l  , [] -> l
-   | [] ,l   -> l
-   | (i,e)::l1,(i',e')::l2 -> 
-      match Pervasives.compare i i' with
-       | 0 -> cons i ( e +/ e') (add  l1 l2)
-       | -1 -> (i,e) :: (add l1 m2)
-       |  1 -> (i', e') :: (add m1 l2)
-       | _  -> assert false
-
-
-
-
- let add t1 t2 = 
-  assert (check t1 && check t2);
-  let res = add t1 t2 in
-   assert (check res);
-   res
-
-
- let rec dotp (t1:t) (t2:t) =
-  match t1, t2 with
-   | [] , _ -> zero_num
-   | _  , [] -> zero_num
-   | (i,e)::l1 , (i',e')::l2 -> 
-      match Pervasives.compare i i' with
-       | 0 -> (e */ e') +/ (dotp l1 l2)
-       | -1 -> dotp l1 t2
-       |  1 -> dotp t1 l2
-       |  _ -> assert false
-	   
- let dotp t1 t2 = 
-  assert (check t1 && check t2) ; dotp t1 t2
-
- let add_mul n t1 t2 =
-  match n with
-   | Int 0 -> t2
-   | Int 1  -> add t1 t2
-   |  _     -> map2 (fun x y -> (n */ x) +/ y) t1 t2
-
- let add_mul n (t1:t) (t2:t) =
-  match n with
-   | Int 0 -> t2
-   | Int 1  -> add t1 t2
-   |   _    -> 
-	let rec xadd_mul m1 m2 = 
-	 match m1, m2 with
-	  | [] , [] -> []
-	  | _  , [] -> mul n m1
-	  | [] , _   -> m2
-	  | (i,e)::l1,(i',e')::l2 -> 
-	     match Pervasives.compare i i' with
-	      | 0 -> cons i ( n */ e +/ e') (xadd_mul  l1 l2)
-	      | -1 -> (i,n */ e) :: (xadd_mul l1 m2)
-	      |  1 -> (i', e') :: (xadd_mul m1 l2)
-	      | _  -> assert false in
-	 xadd_mul t1 t2
-	  
-
-
-
- let lin_comb n1 t1 n2 t2 = 
-  match n1,n2 with
-   | Int 0 , _  -> mul n2 t2
-   | Int 1 , _ -> add_mul n2 t2 t1
-   |   _   , Int 0 ->  mul n1 t1
-   |   _   , Int 1 -> add_mul n1 t1 t2  
-   |       _       -> (*map2 (fun x y -> (n1 */ x) +/ (n2 */ y)) t1 t2*)
-	    let rec xlin_comb m1 m2 = 
-	     match m1, m2 with
-	      | [] , [] -> []
-	      | _  , [] -> mul n1 m1
-	      | [] , _   -> mul n2 m2
-	      | (i,e)::l1,(i',e')::l2 -> 
-		 match Pervasives.compare i i' with
-		  | 0 -> cons i ( n1 */ e +/ n2 */ e') (xlin_comb  l1 l2)
-		  | -1 -> (i,n1 */ e) :: (xlin_comb l1 m2)
-		  |  1 -> (i', n2 */ e') :: (xlin_comb m1 l2)
-		  | _  -> assert false in
-	     xlin_comb t1 t2
-
-
-
-
-
- let lin_comb n1 t1 n2 t2 = 
-  assert (check t1 && check t2);
-  let res = lin_comb n1 t1 n2 t2 in
-   assert (check res); res
-
- let hash = Mutils.Cmp.hash_list (fun (x,y) -> (Hashtbl.hash x) lxor (int_of_num y))
-
-
- let compare = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical 
-  [
-   (fun () -> Pervasives.compare (fst x) (fst y));
-   (fun () -> compare_num (snd x) (snd y))]) 
-
- (* 
-    let compare (x:t) (y:t) = 
-    let rec xcompare acc1 acc2 x y = 
-    match x , y with
-    | [] , [] -> xcomp acc1 acc2
-    | [] , _  -> -1
-    |  _ , [] -> 1
-    | (i,n1)::l1 , (j,n2)::l2 -> 
-    match Pervasives.compare i j with
-    | 0 -> xcompare (n1::acc1) (n2::acc2) l1 l2
-    | c -> c
-    and xcomp acc1 acc2 = Mutils.Cmp.compare_list compare_num acc1 acc2 in
-    xcompare [] [] x y
- *)
-
- type it = t
-
- let iterator (x:t) : it =  x
-  
- let element l = failwith "Not_implemented"
-
- (* TODO: Buffer! *)
- type status = Pos | Neg
-
- let status l = List.map (fun (i,e) -> 
-  match compare_num e (Int 0) with
-   | 1 -> i,Pos
-   | -1 -> i,Neg
-   | _  -> assert false) l
-
- let from_list (l: num list) = 
-  let rec xfrom_list i l = 
-   match l with
-    | [] -> []
-    | e::l ->  
-       if e <>/ Int 0 
-       then (i,e)::(xfrom_list (i+1) l)
-       else xfrom_list (i+1) l in
-   
-  let res = xfrom_list 0 l in
-   assert (check res) ; res
-
-
- let gcd m = 
-  let res = List.fold_left (fun x (i,e) -> Big_int.gcd_big_int  x (Mutils.numerator e))  Big_int.zero_big_int m in
-   if Big_int.compare_big_int res Big_int.zero_big_int = 0 
-   then Big_int.unit_big_int else res
-
- let to_list m = 
-  let rec xto_list i l =
-   match l with
-    | [] -> []
-    |	(x,v)::l' -> 
-	 if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in
-   xto_list 0 m
-
- let to_list l =
-  assert (check l); 
-  to_list l
-
-
- let string l = VList.string (to_list l)
-  
-end
diff --git a/test-suite/micromega/zomicron.v b/test-suite/micromega/zomicron.v
index 2b40f6c93..3087db1ad 100644
--- a/test-suite/micromega/zomicron.v
+++ b/test-suite/micromega/zomicron.v
@@ -20,8 +20,9 @@ Proof.
   lia.
 Qed.
 
-Lemma omega_nightmare : forall x y, 27 <= 11 * x + 13 * y <= 45 -> 7 * x - 9 * y = 4 -> -10 <= 7 * x - 9 * y <= 4 -> False.
+Lemma omega_nightmare : forall x y, 27 <= 11 * x + 13 * y <= 45 ->  -10 <= 7 * x - 9 * y <= 4 -> False.
 Proof.
   intros ; intuition auto.
   lia.
 Qed.  
+