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

(* We take as input a list of polynomials [p1...pn] and return an unfeasibility
   certificate polynomial. *)

(*open Micromega.Polynomial*)
open Big_int
open Num

module Mc = Micromega
module Ml2C = Mutils.CamlToCoq
module C2Ml = Mutils.CoqToCaml

let (<+>) = add_num
let (<->) = minus_num
let (<*>) = mult_num

type var = Mc.positive

module Monomial :
sig
 type t
 val const : t
 val var : var -> t
 val find : var -> t -> int
 val mult : var -> t -> t
 val prod : t -> t -> t
 val compare : t -> t -> int
 val pp : out_channel -> t -> unit
 val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a
end
 =
struct
 (* A monomial is represented by a multiset of variables  *)
 module Map = Map.Make(struct type t = var let compare = Pervasives.compare end)
 open Map
 
 type t = int Map.t

 (* The monomial that corresponds to a constant *)
 let const = Map.empty
  
 (* The monomial 'x' *)
 let var x = Map.add x 1 Map.empty

 (* Get the degre of a variable in a monomial *)
 let find x m = try find x m with Not_found -> 0
  
 (* Multiply a monomial by a variable *)
 let mult x m = add x (  (find x m) + 1) m
  
 (* Product of monomials *)
 let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2
  
 (* Total ordering of monomials *)
 let compare m1 m2 = Map.compare Pervasives.compare m1 m2

 let pp o m = Map.iter (fun k v -> 
  if v = 1 then Printf.fprintf o "x%i." (C2Ml.index k)
  else     Printf.fprintf o "x%i^%i." (C2Ml.index k) v) m

 let fold = fold

end


module  Poly :
 (* A polynomial is a map of monomials *)
 (* 
    This is probably a naive implementation 
    (expected to be fast enough - Coq is probably the bottleneck)
    *The new ring contribution is using a sparse Horner representation.
 *)
sig
 type t
 val get : Monomial.t -> t -> num
 val variable : var -> t
 val add : Monomial.t -> num -> t -> t
 val constant : num -> t
 val mult : Monomial.t -> num -> t -> t
 val product : t -> t -> t
 val addition : t -> t -> t
 val uminus : t -> t
 val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a
 val pp : out_channel -> t -> unit
 val compare : t -> t -> int
end =
struct
 (*normalisation bug : 0*x ... *)
 module P = Map.Make(Monomial)
 open P

 type t = num P.t

 let pp o p = P.iter (fun k v -> 
  if compare_num v (Int 0) <> 0
  then 
   if Monomial.compare Monomial.const k = 0
   then         Printf.fprintf o "%s " (string_of_num v)
   else Printf.fprintf o "%s*%a " (string_of_num v) Monomial.pp k) p  

 (* Get the coefficient of monomial mn *)
 let get : Monomial.t -> t -> num = 
  fun mn p -> try find mn p with Not_found -> (Int 0)


 (* The polynomial 1.x *)
 let variable : var -> t =
  fun  x ->  add (Monomial.var x) (Int 1) empty
   
 (*The constant polynomial *)
 let constant : num -> t =
  fun c -> add (Monomial.const) c empty

 (* The addition of a monomial *)

 let add : Monomial.t -> num -> t -> t =
  fun mn v p -> 
   let vl = (get mn p) <+> v in
    add mn vl p


 (** Design choice: empty is not a polynomial 
     I do not remember why .... 
 **)

 (* The product by a monomial *)
 let mult : Monomial.t -> num -> t -> t =
  fun mn v p -> 
   fold (fun mn' v' res -> P.add (Monomial.prod mn mn') (v<*>v') res) p empty


 let  addition : t -> t -> t =
  fun p1 p2 -> fold (fun mn v p -> add mn v p) p1 p2
   

 let product : t -> t -> t =
  fun p1 p2 -> 
   fold (fun mn v res -> addition (mult mn v p2) res ) p1 empty


 let uminus : t -> t =
  fun p -> map (fun v -> minus_num v) p

 let fold = P.fold

 let compare = compare compare_num
end

open Mutils
type 'a number_spec = {
 bigint_to_number : big_int -> 'a;
 number_to_num : 'a  -> num;
 zero : 'a;
 unit : 'a;
 mult : 'a -> 'a -> 'a;
 eqb  : 'a -> 'a -> Mc.bool
}

let z_spec = {
 bigint_to_number = Ml2C.bigint ;
 number_to_num = (fun x -> Big_int (C2Ml.z_big_int x));
 zero = Mc.Z0;
 unit = Mc.Zpos Mc.XH;
 mult = Mc.zmult;
 eqb  = Mc.zeq_bool
}
 

let q_spec = {
 bigint_to_number = (fun x -> {Mc.qnum = Ml2C.bigint x; Mc.qden = Mc.XH});
 number_to_num = C2Ml.q_to_num;
 zero = {Mc.qnum = Mc.Z0;Mc.qden = Mc.XH};
 unit = {Mc.qnum =  (Mc.Zpos Mc.XH) ; Mc.qden = Mc.XH};
 mult = Mc.qmult;
 eqb  = Mc.qeq_bool
}

let r_spec = z_spec




let dev_form n_spec  p =
 let rec dev_form p = 
  match p with
   | Mc.PEc z ->  Poly.constant (n_spec.number_to_num z)
   | Mc.PEX v ->  Poly.variable v
   | Mc.PEmul(p1,p2) -> 
      let p1 = dev_form p1 in
      let p2 = dev_form p2 in
       Poly.product p1 p2 
   | Mc.PEadd(p1,p2) -> Poly.addition (dev_form p1) (dev_form p2)
   | Mc.PEopp p ->  Poly.uminus (dev_form p)
   | Mc.PEsub(p1,p2) ->  Poly.addition (dev_form p1) (Poly.uminus (dev_form p2))
   | Mc.PEpow(p,n)   ->  
      let p = dev_form p in
      let n = C2Ml.n n in
      let rec pow n = 
       if n = 0 
       then Poly.constant (n_spec.number_to_num n_spec.unit)
       else Poly.product p (pow (n-1)) in
       pow n in
  dev_form p


let monomial_to_polynomial mn = 
 Monomial.fold 
  (fun v i acc -> 
   let mn = if i = 1 then Mc.PEX v else Mc.PEpow (Mc.PEX v ,Ml2C.n i) in
    if acc = Mc.PEc (Mc.Zpos Mc.XH)
    then mn 
    else Mc.PEmul(mn,acc))
  mn 
  (Mc.PEc (Mc.Zpos Mc.XH))
  
let list_to_polynomial vars l = 
 assert (List.for_all (fun x -> ceiling_num x =/ x) l);
 let var x = monomial_to_polynomial (List.nth vars x)  in 
 let rec xtopoly p i = function
  | [] -> p
  | c::l -> if c =/  (Int 0) then xtopoly p (i+1) l 
    else let c = Mc.PEc (Ml2C.bigint (numerator c)) in
    let mn = 
     if c =  Mc.PEc (Mc.Zpos Mc.XH)
     then var i
     else Mc.PEmul (c,var i) in
    let p' = if p = Mc.PEc Mc.Z0 then mn else
      Mc.PEadd (mn, p) in
     xtopoly p' (i+1) l in
  
  xtopoly (Mc.PEc Mc.Z0) 0 l

let rec fixpoint f x =
 let y' = f x in
  if y' = x then y'
  else fixpoint f y'








let  rec_simpl_cone n_spec e = 
 let simpl_cone = 
  Mc.simpl_cone n_spec.zero n_spec.unit n_spec.mult n_spec.eqb in

 let rec rec_simpl_cone  = function
 | Mc.S_Mult(t1, t2) -> 
    simpl_cone  (Mc.S_Mult (rec_simpl_cone t1, rec_simpl_cone t2))
 | Mc.S_Add(t1,t2)  -> 
    simpl_cone (Mc.S_Add (rec_simpl_cone t1, rec_simpl_cone t2))
 |  x           -> simpl_cone x in
  rec_simpl_cone e
   
   
let simplify_cone n_spec c = fixpoint (rec_simpl_cone n_spec) c
 
type cone_prod = 
  Const of cone 
  | Ideal of cone *cone 
  | Mult of cone * cone 
  | Other of cone
and cone =   Mc.zWitness



let factorise_linear_cone c =
 
 let rec cone_list  c l = 
  match c with
   | Mc.S_Add (x,r) -> cone_list  r (x::l)
   |  _        ->  c :: l in
  
 let factorise c1 c2 =
  match c1 , c2 with
   | Mc.S_Ideal(x,y) , Mc.S_Ideal(x',y') -> 
      if x = x' then Some (Mc.S_Ideal(x, Mc.S_Add(y,y'))) else None
   | Mc.S_Mult(x,y) , Mc.S_Mult(x',y') -> 
      if x = x' then Some (Mc.S_Mult(x, Mc.S_Add(y,y'))) else None
   |  _     -> None in
  
 let rec rebuild_cone l pending  =
  match l with
   | [] -> (match pending with
      | None -> Mc.S_Z
      | Some p -> p
     )
   | e::l -> 
      (match pending with
       | None -> rebuild_cone l (Some e) 
       | Some p -> (match factorise p e with
          | None -> Mc.S_Add(p, rebuild_cone l (Some e))
          | Some f -> rebuild_cone l (Some f) )
      ) in

  (rebuild_cone (List.sort Pervasives.compare (cone_list c [])) None)



(* The binding with Fourier might be a bit obsolete 
   -- how does it handle equalities ? *)

(* Certificates are elements of the cone such that P = 0  *)

(* To begin with, we search for certificates of the form:
   a1.p1 + ... an.pn + b1.q1 +... + bn.qn + c = 0   
   where pi >= 0 qi > 0
   ai >= 0 
   bi >= 0
   Sum bi + c >= 1
   This is a linear problem: each monomial is considered as a variable.
   Hence, we can use fourier.

   The variable c is at index 0
*)

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 Vect = Fourier.Vect
open Fourier.Cstr

(* fold_left followed by a rev ! *)

let constrain_monomial mn l  = 
 let coeffs = List.fold_left (fun acc p -> (Poly.get mn p)::acc) [] l in
  if mn = Monomial.const
  then  
   { coeffs = Vect.from_list ((Big_int unit_big_int):: (List.rev coeffs)) ; 
     op = Eq ; 
     cst = Big_int zero_big_int  }
  else
   { coeffs = Vect.from_list ((Big_int zero_big_int):: (List.rev coeffs)) ; 
     op = Eq ; 
     cst = Big_int zero_big_int  }

    
let positivity l = 
 let rec xpositivity i l = 
  match l with
   | [] -> []
   | (_,Mc.Equal)::l -> xpositivity (i+1) l
   | (_,_)::l -> 
      {coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ; 
       op = Ge ; 
       cst = Int 0 }  :: (xpositivity (i+1) l)
 in
  xpositivity 0 l


let string_of_op = function
 | Mc.Strict -> "> 0" 
 | Mc.NonStrict -> ">= 0" 
 | Mc.Equal -> "= 0"
 | Mc.NonEqual -> "<> 0"



(* If the certificate includes at least one strict inequality, 
   the obtained polynomial can also be 0 *)
let build_linear_system l =

 (* Gather the monomials:  HINT add up of the polynomials *)
 let l' = List.map fst l in
 let monomials = 
  List.fold_left (fun acc p -> Poly.addition p acc) (Poly.constant (Int 0)) l'
 in  (* For each monomial, compute a constraint *)
 let s0 = 
  Poly.fold (fun mn _ res -> (constrain_monomial mn l')::res) monomials [] in
  (* I need at least something strictly positive *)
 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 
                              | _ -> Big_int zero_big_int) l));
  op = Ge ; cst = Big_int unit_big_int } in
  (* Add the positivity constraint *)
  {coeffs = Vect.from_list ([Big_int unit_big_int]) ; 
   op = Ge ; 
   cst = Big_int zero_big_int}::(strict::(positivity l)@s0)


let big_int_to_z = Ml2C.bigint
 
(* For Q, this is a pity that the certificate has been scaled 
   -- at a lower layer, certificates are using nums... *)
let make_certificate n_spec cert li = 
 let bint_to_cst = n_spec.bigint_to_number in
  match cert with
   | [] -> None
   | e::cert' -> 
      let cst = match compare_big_int e zero_big_int with
       | 0 -> Mc.S_Z
       | 1 ->  Mc.S_Pos (bint_to_cst e) 
       | _ -> failwith "positivity error" 
      in
      let rec scalar_product cert l =
       match cert with
        | [] -> Mc.S_Z
        | c::cert -> match l with
           | [] -> failwith "make_certificate(1)"
           | i::l ->  
              let r = scalar_product cert l in
               match compare_big_int c  zero_big_int with
                | -1 -> Mc.S_Add (
                   Mc.S_Ideal (Mc.PEc ( bint_to_cst c), Mc.S_In (Ml2C.nat i)), 
                   r)
                | 0  -> r
                | _ ->  Mc.S_Add (
                   Mc.S_Mult (Mc.S_Pos (bint_to_cst c), Mc.S_In (Ml2C.nat i)),
                   r) in
       
      Some ((factorise_linear_cone 
              (simplify_cone n_spec (Mc.S_Add (cst, scalar_product cert' li)))))


exception Found of Monomial.t
 
let raw_certificate l = 
 let sys = build_linear_system l in
  try 
   match Fourier.find_point sys with
    | None -> None
    | Some cert ->  Some (rats_to_ints (Vect.to_list cert)) 
       (* should not use rats_to_ints *)
  with x -> 
   if debug 
   then (Printf.printf "raw certificate %s" (Printexc.to_string x);   
         flush stdout) ;
   None


let simple_linear_prover to_constant l =
 let (lc,li) = List.split l in
  match raw_certificate lc with
   |  None -> None (* No certificate *)
   | Some cert -> make_certificate to_constant cert li
      
      

let linear_prover n_spec l  =
 let li = List.combine l (interval 0 (List.length l -1)) in
 let (l1,l') = List.partition 
  (fun (x,_) -> if snd' x = Mc.NonEqual then true else false) li in
 let l' = List.map 
  (fun (c,i) -> let (Mc.Pair(x,y)) = c in 
                 match y with 
                   Mc.NonEqual -> failwith "cannot happen" 
                  |  y -> ((dev_form n_spec x, y),i)) l' in
  
  simple_linear_prover n_spec l' 


let linear_prover n_spec l  =
 try linear_prover n_spec l with
   x -> (print_string (Printexc.to_string x); None)

(* zprover.... *)

(* I need to gather the set of variables --->
   Then go for fold 
   Once I have an interval, I need a certificate : 2 other fourier elims.
   (I could probably get the certificate directly 
   as it is done in the fourier contrib.)
*)

let make_linear_system l =
 let l' = List.map fst l in
 let monomials = List.fold_left (fun acc p -> Poly.addition p acc) 
  (Poly.constant (Int 0)) l' in
 let monomials = Poly.fold 
  (fun mn _ l -> if mn = Monomial.const then l else mn::l) monomials [] in
  (List.map (fun (c,op) -> 
   {coeffs = Vect.from_list (List.map (fun mn ->  (Poly.get mn c)) monomials) ; 
    op = op ; 
    cst = minus_num ( (Poly.get Monomial.const c))}) 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
let popp x = Mc.PEopp x
 
let debug = false
 
(* keep track of enumerated vectors *)
let rec mem p x  l = 
 match l with  [] -> false | e::l -> if p x e then true else mem p x l

let rec remove_assoc p x l = 
 match l with [] -> [] | e::l -> if p x (fst e) then
  remove_assoc p x l else e::(remove_assoc p x l) 

let eq x y = Vect.compare x y = 0

let  remove e  l  = List.fold_left (fun l x -> if eq x e then l else x::l) [] l


(* The prover is (probably) incomplete --
   only searching for naive cutting planes *)

let candidates sys = 
 let ll = List.fold_right (
  fun (Mc.Pair(e,k)) r -> 
   match k with 
    | Mc.NonStrict -> (dev_form z_spec e , Ge)::r
    | Mc.Equal     ->  (dev_form z_spec e , Eq)::r 
       (* we already know the bound -- don't compute it again *)
    | _     -> failwith "Cannot happen candidates") sys [] in

 let (sys,var_mn) = make_linear_system ll in
 let vars = mapi (fun _ i -> Vect.set i (Int 1) Vect.null) var_mn in
  (List.fold_left (fun l cstr ->   
   let gcd = Big_int (Vect.gcd cstr.coeffs) in
    if gcd =/ (Int 1) && cstr.op = Eq 
    then l 
    else  (Vect.mul (Int 1 // gcd) cstr.coeffs)::l) [] sys) @ vars


let rec xzlinear_prover planes sys = 
 match linear_prover z_spec sys with
  | Some prf -> Some (Mc.RatProof prf)
  | None     -> (* find the candidate with the smallest range *)
     (* Grrr - linear_prover is also calling 'make_linear_system' *)
     let ll = List.fold_right (fun (Mc.Pair(e,k)) r -> match k with 
       Mc.NonEqual -> r  
      | k -> (dev_form z_spec e , 
             match k with
               Mc.NonStrict -> Ge 
              | Mc.Equal              -> Eq
              | Mc.Strict | Mc.NonEqual -> failwith "Cannot happen") :: r) sys [] in
     let (ll,var) = make_linear_system ll in
     let candidates = List.fold_left (fun  acc vect -> 
      match Fourier.optimise vect ll with
       | None -> acc
       | Some i -> 
(*	  Printf.printf "%s in %s\n" (Vect.string vect) (string_of_intrvl i) ; *)
	  flush stdout ; 
	  (vect,i) ::acc)  [] planes in

     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 
      with
       | (x,Itv(Some i, Some j))  -> Some(i,x,j)
       | (x,Point n) ->  Some(n,x,n)
       |   x        ->   None (* This might be a cutting plane *)
      in
      match smallest_interval with
       | Some (lb,e,ub) -> 
          let (lbn,lbd) = 
           (Ml2C.bigint (sub_big_int (numerator lb)  unit_big_int),
           Ml2C.bigint (denominator lb)) in
          let (ubn,ubd) = 
           (Ml2C.bigint (add_big_int unit_big_int (numerator ub)) , 
           Ml2C.bigint (denominator ub)) in
          let expr = list_to_polynomial var (Vect.to_list e) in
           (match 
             (*x <= ub ->  x  > ub *)
             linear_prover  z_spec 
              (Mc.Pair(pplus (pmult (pconst ubd) expr) (popp (pconst  ubn)),
                      Mc.NonStrict) :: sys),
            (* lb <= x  -> lb > x *)
            linear_prover z_spec 
             (Mc.Pair( pplus  (popp (pmult  (pconst lbd) expr)) (pconst lbn) ,  
                     Mc.NonStrict)::sys) 
            with
             | Some cub , Some clb  ->   
                (match zlinear_enum  (remove e planes)   expr 
                  (ceiling_num lb)  (floor_num ub) sys 
                 with
                  | None -> None
                  | Some prf -> 
                     Some (Mc.EnumProof(Ml2C.q lb,expr,Ml2C.q ub,clb,cub,prf)))
             | _ -> None
           )
       |  _ -> None
and zlinear_enum  planes expr clb cub l = 
 if clb >/  cub
 then Some Mc.Nil
 else 
  let pexpr = pplus (popp (pconst (Ml2C.bigint (numerator clb)))) expr in
  let sys' = (Mc.Pair(pexpr, Mc.Equal))::l in
   (*let enum =  *)
   match xzlinear_prover planes sys' with
    | None -> if debug then print_string "zlp?"; None
    | Some prf -> if debug then print_string "zlp!";
    match zlinear_enum planes expr (clb +/ (Int 1)) cub l with
     | None -> None
     | Some prfl -> Some (Mc.Cons(prf,prfl))

let zlinear_prover sys =  
 let candidates = candidates sys in
 (*  Printf.printf "candidates %d" (List.length candidates) ; *) 
 xzlinear_prover candidates sys

open Sos

let rec scale_term t = 
 match t with
  | Zero    -> unit_big_int , Zero
  | Const n ->  (denominator n) , Const (Big_int (numerator n))
  | Var n   -> unit_big_int , Var n
  | Inv _   -> failwith "scale_term : not implemented"
  | Opp t   -> let s, t = scale_term t in s, Opp t
  | Add(t1,t2) -> let s1,y1 = scale_term t1 and s2,y2 = scale_term t2 in
    let g = gcd_big_int s1 s2 in
    let s1' = div_big_int s1 g in
    let s2' = div_big_int s2 g in
    let e = mult_big_int g (mult_big_int s1' s2') in
     if (compare_big_int e unit_big_int) = 0
     then (unit_big_int, Add (y1,y2))
     else 	e, Add (Mul(Const (Big_int s2'), y1),
		       Mul (Const (Big_int s1'), y2))
  | Sub _ -> failwith "scale term: not implemented"
  | Mul(y,z) ->       let s1,y1 = scale_term y and s2,y2 = scale_term z in
		       mult_big_int s1 s2 , Mul (y1, y2)
  |  Pow(t,n) -> let s,t = scale_term t in
		  power_big_int_positive_int s  n , Pow(t,n)
  |   _ -> failwith "scale_term : not implemented"

let scale_term t =
 let (s,t') = scale_term t in
  s,t'


let get_index_of_ith_match f i l  =
 let rec get j res l =
  match l with
   | [] -> failwith "bad index"
   | e::l -> if f e
     then 
       (if j = i then res else get (j+1) (res+1) l )
     else get j (res+1) l in
  get 0 0 l


let rec scale_certificate pos = match pos with
 | Axiom_eq i ->  unit_big_int , Axiom_eq i
 | Axiom_le i ->  unit_big_int , Axiom_le i
 | Axiom_lt i ->   unit_big_int , Axiom_lt i
 | Monoid l   -> unit_big_int , Monoid l
 | Rational_eq n ->  (denominator n) , Rational_eq (Big_int (numerator n))
 | Rational_le n ->  (denominator n) , Rational_le (Big_int (numerator n))
 | Rational_lt n ->  (denominator n) , Rational_lt (Big_int (numerator n))
 | Square t -> let s,t' =  scale_term t in 
		mult_big_int s s , Square t'
 | Eqmul (t, y) -> let s1,y1 = scale_term t and s2,y2 = scale_certificate y in
		    mult_big_int s1 s2 , Eqmul (y1,y2)
 | Sum (y, z) -> let s1,y1 = scale_certificate y 
   and s2,y2 = scale_certificate z in
   let g = gcd_big_int s1 s2 in
   let s1' = div_big_int s1 g in
   let s2' = div_big_int s2 g in
    mult_big_int g (mult_big_int s1' s2'), 
    Sum (Product(Rational_le (Big_int s2'), y1),
	Product (Rational_le (Big_int s1'), y2))
 | Product (y, z) -> 
    let s1,y1 = scale_certificate y and s2,y2 = scale_certificate z in
     mult_big_int s1 s2 , Product (y1,y2)


open Micromega
 let rec term_to_q_expr = function
  | Const n ->  PEc (Ml2C.q n)
  | Zero   ->  PEc ( Ml2C.q (Int 0))
  | Var s   ->  PEX (Ml2C.index 
		      (int_of_string (String.sub s 1 (String.length s - 1))))
  | Mul(p1,p2) ->  PEmul(term_to_q_expr p1, term_to_q_expr p2)
  | Add(p1,p2) ->   PEadd(term_to_q_expr p1, term_to_q_expr p2)
  | Opp p ->   PEopp (term_to_q_expr p)
  | Pow(t,n) ->  PEpow (term_to_q_expr t,Ml2C.n n)
  | Sub(t1,t2) ->  PEsub (term_to_q_expr t1,  term_to_q_expr t2)
  | _ -> failwith "term_to_q_expr: not implemented"

let  q_cert_of_pos  pos = 
 let rec _cert_of_pos = function
   Axiom_eq i ->  Mc.S_In (Ml2C.nat i)
  | Axiom_le i ->  Mc.S_In (Ml2C.nat i)
  | Axiom_lt i ->  Mc.S_In (Ml2C.nat i)
  | Monoid l  -> Mc.S_Monoid (Ml2C.list Ml2C.nat l)
  | Rational_eq n | Rational_le n | Rational_lt n -> 
     if compare_num n (Int 0) = 0 then Mc.S_Z else
      Mc.S_Pos (Ml2C.q   n)
  | Square t -> Mc.S_Square (term_to_q_expr  t)
  | Eqmul (t, y) -> Mc.S_Ideal(term_to_q_expr t, _cert_of_pos y)
  | Sum (y, z) -> Mc.S_Add  (_cert_of_pos y, _cert_of_pos z)
  | Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in
  simplify_cone q_spec (_cert_of_pos pos)


 let rec term_to_z_expr = function
  | Const n ->  PEc (Ml2C.bigint (big_int_of_num n))
  | Zero   ->  PEc ( Z0)
  | Var s   ->  PEX (Ml2C.index 
		      (int_of_string (String.sub s 1 (String.length s - 1))))
  | Mul(p1,p2) ->  PEmul(term_to_z_expr p1, term_to_z_expr p2)
  | Add(p1,p2) ->   PEadd(term_to_z_expr p1, term_to_z_expr p2)
  | Opp p ->   PEopp (term_to_z_expr p)
  | Pow(t,n) ->  PEpow (term_to_z_expr t,Ml2C.n n)
  | Sub(t1,t2) ->  PEsub (term_to_z_expr t1,  term_to_z_expr t2)
  | _ -> failwith "term_to_z_expr: not implemented"

let  z_cert_of_pos  pos = 
 let s,pos = (scale_certificate pos) in
 let rec _cert_of_pos = function
   Axiom_eq i ->  Mc.S_In (Ml2C.nat i)
  | Axiom_le i ->  Mc.S_In (Ml2C.nat i)
  | Axiom_lt i ->  Mc.S_In (Ml2C.nat i)
  | Monoid l  -> Mc.S_Monoid (Ml2C.list Ml2C.nat l)
  | Rational_eq n | Rational_le n | Rational_lt n -> 
     if compare_num n (Int 0) = 0 then Mc.S_Z else
      Mc.S_Pos (Ml2C.bigint (big_int_of_num  n))
  | Square t -> Mc.S_Square (term_to_z_expr  t)
  | Eqmul (t, y) -> Mc.S_Ideal(term_to_z_expr t, _cert_of_pos y)
  | Sum (y, z) -> Mc.S_Add  (_cert_of_pos y, _cert_of_pos z)
  | Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in
  simplify_cone z_spec (_cert_of_pos pos)

