(************************************************************************)
(*         *   The Coq Proof Assistant / The Coq Development Team       *)
(*  v      *         Copyright INRIA, CNRS and contributors             *)
(* <O___,, * (see version control and CREDITS file for authors & dates) *)
(*   \VV/  **************************************************************)
(*    //   *    This file is distributed under the terms of the         *)
(*         *     GNU Lesser General Public License Version 2.1          *)
(*         *     (see LICENSE file for the text of the license)         *)
(************************************************************************)
(*                                                                      *)
(* Micromega: A reflexive tactic using the Positivstellensatz           *)
(*                                                                      *)
(*  Frédéric Besson (Irisa/Inria) 2006-20018                            *)
(*                                                                      *)
(************************************************************************)

open NumCompat
open Q.Notations
open Mutils
module Mc = Micromega

let max_nb_cstr = ref max_int

type var = int

let debug = false
let ( <+> ) = ( +/ )
let ( <*> ) = ( */ )

module Monomial : sig
  type t

  val const : t
  val is_const : t -> bool
  val var : var -> t
  val is_var : t -> bool
  val get_var : t -> var option
  val prod : t -> t -> t
  val exp : t -> int -> t
  val div : t -> t -> t * int
  val compare : t -> t -> int
  val pp : out_channel -> t -> unit
  val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a
  val sqrt : t -> t option
  val variables : t -> ISet.t
  val degree : t -> int
end = struct
  (* A monomial is represented by a multiset of variables  *)
  module Map = Map.Make (Int)
  open Map

  type t = int Map.t

  let degree m = Map.fold (fun _ i d -> i + d) m 0

  let is_singleton m =
    try
      let k, v = choose m in
      let l, e, r = split k m in
      if is_empty l && is_empty r then Some (k, v) else None
    with Not_found -> None

  let pp o m =
    let pp_elt o (k, v) =
      if v = 1 then Printf.fprintf o "x%i" k else Printf.fprintf o "x%i^%i" k v
    in
    let rec pp_list o l =
      match l with
      | [] -> ()
      | [e] -> pp_elt o e
      | e :: l -> Printf.fprintf o "%a*%a" pp_elt e pp_list l
    in
    pp_list o (Map.bindings m)

  (* The monomial that corresponds to a constant *)
  let const = Map.empty
  let sum_degree m = Map.fold (fun _ n s -> s + n) m 0

  (* Total ordering of monomials *)
  let compare : t -> t -> int =
   fun m1 m2 ->
    let s1 = sum_degree m1 and s2 = sum_degree m2 in
    if Int.equal s1 s2 then Map.compare Int.compare m1 m2 else Int.compare s1 s2

  let is_const m = m = Map.empty

  (* The monomial 'x' *)
  let var x = Map.add x 1 Map.empty

  let is_var m =
    match is_singleton m with None -> false | Some (_, i) -> i = 1

  let get_var m =
    match is_singleton m with
    | None -> None
    | Some (k, i) -> if i = 1 then Some k else None

  let sqrt m =
    if is_const m then None
    else
      try
        Some
          (Map.fold
             (fun v i acc ->
               let i' = i / 2 in
               if i mod 2 = 0 then add v i' acc else raise Not_found)
             m const)
      with Not_found -> None

  (* Get the degre of a variable in a monomial *)
  let find x m = try find x m with Not_found -> 0

  (* Product of monomials *)
  let prod m1 m2 = Map.fold (fun k d m -> add k (find k m + d) m) m1 m2

  let exp m n =
    let rec exp acc n = if n = 0 then acc else exp (prod acc m) (n - 1) in
    exp const n

  (*  [div m1 m2 = mr,n] such that mr * (m2)^n = m1 *)
  let div m1 m2 =
    let n =
      fold
        (fun x i n ->
          let i' = find x m1 in
          let nx = i' / i in
          min n nx)
        m2 max_int
    in
    let mr =
      fold
        (fun x i' m ->
          let i = find x m2 in
          let ir = i' - (i * n) in
          if ir = 0 then m else add x ir m)
        m1 empty
    in
    (mr, n)

  let variables m = fold (fun v i acc -> ISet.add v acc) m ISet.empty
  let fold = fold
end

module MonMap = struct
  include Map.Make (Monomial)

  let union f =
    merge (fun x v1 v2 ->
        match (v1, v2) with
        | None, None -> None
        | Some v, None | None, Some v -> Some v
        | Some v1, Some v2 -> f x v1 v2)
end

let pp_mon o (m, i) =
  if Monomial.is_const m then
    if Q.zero =/ i then () else Printf.fprintf o "%s" (Q.to_string i)
  else if Q.one =/ i then Monomial.pp o m
  else if Q.minus_one =/ i then Printf.fprintf o "-%a" Monomial.pp m
  else if Q.zero =/ i then ()
  else Printf.fprintf o "%s*%a" (Q.to_string i) Monomial.pp m

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 pp : out_channel -> t -> unit
  val get : Monomial.t -> t -> Q.t
  val variable : var -> t
  val add : Monomial.t -> Q.t -> t -> t
  val constant : Q.t -> t
  val product : t -> t -> t
  val addition : t -> t -> t
  val uminus : t -> t
  val fold : (Monomial.t -> Q.t -> 'a -> 'a) -> t -> 'a -> 'a
  val factorise : var -> t -> t * t
end = struct
  (*normalisation bug : 0*x ... *)
  module P = Map.Make (Monomial)
  open P

  type t = Q.t P.t

  let pp o p = P.iter (fun mn i -> Printf.fprintf o "%a + " pp_mon (mn, i)) p

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

  (* The polynomial 1.x *)
  let variable : var -> t = fun x -> add (Monomial.var x) Q.one empty

  (*The constant polynomial *)
  let constant : Q.t -> t = fun c -> add Monomial.const c empty

  (* The addition of a monomial *)

  let add : Monomial.t -> Q.t -> t -> t =
   fun mn v p ->
    if Q.sign v = 0 then p
    else
      let vl = get mn p <+> v in
      if Q.sign vl = 0 then remove mn p else 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 -> Q.t -> t -> t =
   fun mn v p ->
    if Q.sign v = 0 then constant Q.zero
    else
      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 -> Q.neg v) p
  let fold = P.fold

  let factorise x p =
    let x = Monomial.var x in
    P.fold
      (fun m v (px, cx) ->
        let m1, i = Monomial.div m x in
        if i = 0 then (px, add m v cx)
        else
          let mx = Monomial.prod m1 (Monomial.exp x (i - 1)) in
          (add mx v px, cx))
      p
      (constant Q.zero, constant Q.zero)
end

type vector = Vect.t

type cstr = {coeffs : vector; op : op; cst : Q.t}

and op = Eq | Ge | Gt

exception Strict

let is_strict c = c.op = Gt
let eval_op = function Eq -> ( =/ ) | Ge -> ( >=/ ) | Gt -> ( >/ )
let string_of_op = function Eq -> "=" | Ge -> ">=" | Gt -> ">"

let compare_op o1 o2 =
  match (o1, o2) with
  | Eq, Eq -> 0
  | Eq, _ -> -1
  | _, Eq -> 1
  | Ge, Ge -> 0
  | Ge, _ -> -1
  | _, Ge -> 1
  | Gt, Gt -> 0

let output_cstr o {coeffs; op; cst} =
  Printf.fprintf o "%a %s %s" Vect.pp coeffs (string_of_op op) (Q.to_string cst)

let opMult o1 o2 =
  match (o1, o2) with Eq, _ | _, Eq -> Eq | Ge, _ | _, Ge -> Ge | Gt, Gt -> Gt

let opAdd o1 o2 =
  match (o1, o2) with Eq, x | x, Eq -> x | Gt, x | x, Gt -> Gt | Ge, Ge -> Ge

module LinPoly = struct
  (** A linear polynomial a0 + a1.x1 + ... + an.xn
      By convention, the constant a0 is the coefficient of the variable 0.
   *)

  type t = Vect.t

  module MonT = struct
    module MonoMap = Map.Make (Monomial)
    module IntMap = Map.Make (Int)

    (** A hash table might be preferable but requires a hash function. *)
    let (index_of_monomial : int MonoMap.t ref) = ref MonoMap.empty

    let (monomial_of_index : Monomial.t IntMap.t ref) = ref IntMap.empty
    let fresh = ref 0

    let reserve vr =
      if !fresh > vr then failwith (Printf.sprintf "Cannot reserve %i" vr)
      else fresh := vr + 1

    let safe_reserve vr = if !fresh > vr then () else fresh := vr + 1

    let get_fresh () =
      let vr = !fresh in
      incr fresh; vr

    let register m =
      try MonoMap.find m !index_of_monomial
      with Not_found ->
        let res = !fresh in
        index_of_monomial := MonoMap.add m res !index_of_monomial;
        monomial_of_index := IntMap.add res m !monomial_of_index;
        incr fresh;
        res

    let retrieve i = IntMap.find i !monomial_of_index

    let clear () =
      index_of_monomial := MonoMap.empty;
      monomial_of_index := IntMap.empty;
      fresh := 0;
      ignore (register Monomial.const)

    let _ = register Monomial.const
  end

  let var v = Vect.set (MonT.register (Monomial.var v)) Q.one Vect.null

  let of_monomial m =
    let v = MonT.register m in
    Vect.set v Q.one Vect.null

  let linpol_of_pol p =
    Poly.fold
      (fun mon num vct ->
        let vr = MonT.register mon in
        Vect.set vr num vct)
      p Vect.null

  let pol_of_linpol v =
    Vect.fold
      (fun p vr n -> Poly.add (MonT.retrieve vr) n p)
      (Poly.constant Q.zero) v

  let coq_poly_of_linpol cst p =
    let pol_of_mon m =
      Monomial.fold
        (fun x v p ->
          Mc.PEmul (Mc.PEpow (Mc.PEX (CamlToCoq.positive x), CamlToCoq.n v), p))
        m
        (Mc.PEc (cst Q.one))
    in
    Vect.fold
      (fun acc x v ->
        let mn = MonT.retrieve x in
        Mc.PEadd (Mc.PEmul (Mc.PEc (cst v), pol_of_mon mn), acc))
      (Mc.PEc (cst Q.zero))
      p

  let pp_var o vr =
    try Monomial.pp o (MonT.retrieve vr) (* this is a non-linear variable *)
    with Not_found -> Printf.fprintf o "v%i" vr

  let pp o p = Vect.pp_gen pp_var o p
  let constant c = if Q.sign c = 0 then Vect.null else Vect.set 0 c Vect.null

  let is_linear p =
    Vect.for_all
      (fun v _ ->
        let mn = MonT.retrieve v in
        Monomial.is_var mn || Monomial.is_const mn)
      p

  let is_variable p =
    let (x, v), r = Vect.decomp_fst p in
    if Vect.is_null r && v >/ Q.zero then Monomial.get_var (MonT.retrieve x)
    else None

  let factorise x p =
    let px, cx = Poly.factorise x (pol_of_linpol p) in
    (linpol_of_pol px, linpol_of_pol cx)

  let is_linear_for x p =
    let a, b = factorise x p in
    Vect.is_constant a

  let search_all_linear p l =
    Vect.fold
      (fun acc x v ->
        if p v then
          let x' = MonT.retrieve x in
          match Monomial.get_var x' with
          | None -> acc
          | Some x -> if is_linear_for x l then x :: acc else acc
        else acc)
      [] l

  let get_bound p = Vect.Bound.of_vect p

  let min_list (l : int list) =
    match l with [] -> None | e :: l -> Some (List.fold_left min e l)

  let search_linear p l = min_list (search_all_linear p l)

  let product p1 p2 =
    linpol_of_pol (Poly.product (pol_of_linpol p1) (pol_of_linpol p2))

  let addition p1 p2 = Vect.add p1 p2

  let of_vect v =
    Vect.fold
      (fun acc v vl -> addition (product (var v) (constant vl)) acc)
      Vect.null v

  let variables p =
    Vect.fold
      (fun acc v _ -> ISet.union (Monomial.variables (MonT.retrieve v)) acc)
      ISet.empty p

  let monomials p = Vect.fold (fun acc v _ -> ISet.add v acc) ISet.empty p

  let degree v =
    Vect.fold (fun acc v vl -> max acc (Monomial.degree (MonT.retrieve v))) 0 v

  let pp_goal typ o l =
    let vars =
      List.fold_left
        (fun acc p -> ISet.union acc (variables (fst p)))
        ISet.empty l
    in
    let pp_vars o i =
      ISet.iter (fun v -> Printf.fprintf o "(x%i : %s) " v typ) vars
    in
    Printf.fprintf o "forall %a\n" pp_vars vars;
    List.iteri
      (fun i (p, op) ->
        Printf.fprintf o "(H%i : %a %s 0)\n" i pp p (string_of_op op))
      l;
    Printf.fprintf o ", False\n"

  let collect_square p =
    Vect.fold
      (fun acc v _ ->
        let m = MonT.retrieve v in
        match Monomial.sqrt m with None -> acc | Some s -> MonMap.add s m acc)
      MonMap.empty p
end

module ProofFormat = struct
  type prf_rule =
    | Annot of string * prf_rule
    | Hyp of int
    | Def of int
    | Cst of Q.t
    | Zero
    | Square of Vect.t
    | MulC of Vect.t * prf_rule
    | Gcd of Z.t * prf_rule
    | MulPrf of prf_rule * prf_rule
    | AddPrf of prf_rule * prf_rule
    | CutPrf of prf_rule

  type proof =
    | Done
    | Step of int * prf_rule * proof
    | Split of int * Vect.t * proof * proof
    | Enum of int * prf_rule * Vect.t * prf_rule * proof list
    | ExProof of int * int * int * var * var * var * proof

  (* x = z - t, z >= 0, t >= 0 *)

  let rec output_prf_rule o = function
    | Annot (s, p) -> Printf.fprintf o "(%a)@%s" output_prf_rule p s
    | Hyp i -> Printf.fprintf o "Hyp %i" i
    | Def i -> Printf.fprintf o "Def %i" i
    | Cst c -> Printf.fprintf o "Cst %s" (Q.to_string c)
    | Zero -> Printf.fprintf o "Zero"
    | Square s -> Printf.fprintf o "(%a)^2" Poly.pp (LinPoly.pol_of_linpol s)
    | MulC (p, pr) ->
      Printf.fprintf o "(%a) * (%a)" Poly.pp (LinPoly.pol_of_linpol p)
        output_prf_rule pr
    | MulPrf (p1, p2) ->
      Printf.fprintf o "(%a) * (%a)" output_prf_rule p1 output_prf_rule p2
    | AddPrf (p1, p2) ->
      Printf.fprintf o "%a + %a" output_prf_rule p1 output_prf_rule p2
    | CutPrf p -> Printf.fprintf o "[%a]" output_prf_rule p
    | Gcd (c, p) -> Printf.fprintf o "(%a)/%s" output_prf_rule p (Z.to_string c)

  let rec output_proof o = function
    | Done -> Printf.fprintf o "."
    | Step (i, p, pf) ->
      Printf.fprintf o "%i:= %a\n ; %a" i output_prf_rule p output_proof pf
    | Split (i, v, p1, p2) ->
      Printf.fprintf o "%i:=%a ; { %a } { %a }" i Vect.pp v output_proof p1
        output_proof p2
    | Enum (i, p1, v, p2, pl) ->
      Printf.fprintf o "%i{%a<=%a<=%a}%a" i output_prf_rule p1 Vect.pp v
        output_prf_rule p2 (pp_list ";" output_proof) pl
    | ExProof (i, j, k, x, z, t, pr) ->
      Printf.fprintf o "%i := %i = %i - %i ; %i := %i >= 0 ; %i := %i >= 0 ; %a"
        i x z t j z k t output_proof pr

  module OrdPrfRule = struct
    type t = prf_rule

    let id_of_constr = function
      | Annot _ -> 0
      | Hyp _ -> 1
      | Def _ -> 2
      | Cst _ -> 3
      | Zero -> 4
      | Square _ -> 5
      | MulC _ -> 6
      | Gcd _ -> 7
      | MulPrf _ -> 8
      | AddPrf _ -> 9
      | CutPrf _ -> 10

    let cmp_pair c1 c2 (x1, x2) (y1, y2) =
      match c1 x1 y1 with 0 -> c2 x2 y2 | i -> i

    let rec compare p1 p2 =
      match (p1, p2) with
      | Annot (s1, p1), Annot (s2, p2) ->
        if s1 = s2 then compare p1 p2 else String.compare s1 s2
      | Hyp i, Hyp j -> Int.compare i j
      | Def i, Def j -> Int.compare i j
      | Cst n, Cst m -> Q.compare n m
      | Zero, Zero -> 0
      | Square v1, Square v2 -> Vect.compare v1 v2
      | MulC (v1, p1), MulC (v2, p2) ->
        cmp_pair Vect.compare compare (v1, p1) (v2, p2)
      | Gcd (b1, p1), Gcd (b2, p2) ->
        cmp_pair Z.compare compare (b1, p1) (b2, p2)
      | MulPrf (p1, q1), MulPrf (p2, q2) ->
        cmp_pair compare compare (p1, q1) (p2, q2)
      | AddPrf (p1, q1), AddPrf (p2, q2) ->
        cmp_pair compare compare (p1, q1) (p2, q2)
      | CutPrf p, CutPrf p' -> compare p p'
      | _, _ -> Int.compare (id_of_constr p1) (id_of_constr p2)
  end

  module PrfRuleMap = Map.Make (OrdPrfRule)

  let rec pr_size = function
    | Annot (_, p) -> pr_size p
    | Zero | Square _ -> Q.zero
    | Hyp _ -> Q.one
    | Def _ -> Q.one
    | Cst n -> n
    | Gcd (i, p) -> pr_size p // Q.of_bigint i
    | MulPrf (p1, p2) | AddPrf (p1, p2) -> pr_size p1 +/ pr_size p2
    | CutPrf p -> pr_size p
    | MulC (v, p) -> pr_size p

  let rec pr_rule_max_hyp = function
    | Annot (_, p) -> pr_rule_max_hyp p
    | Hyp i -> i
    | Def i -> -1
    | Cst _ | Zero | Square _ -> -1
    | MulC (_, p) | CutPrf p | Gcd (_, p) -> pr_rule_max_hyp p
    | MulPrf (p1, p2) | AddPrf (p1, p2) ->
      max (pr_rule_max_hyp p1) (pr_rule_max_hyp p2)

  let rec pr_rule_max_def = function
    | Annot (_, p) -> pr_rule_max_hyp p
    | Hyp i -> -1
    | Def i -> i
    | Cst _ | Zero | Square _ -> -1
    | MulC (_, p) | CutPrf p | Gcd (_, p) -> pr_rule_max_def p
    | MulPrf (p1, p2) | AddPrf (p1, p2) ->
      max (pr_rule_max_def p1) (pr_rule_max_def p2)

  let rec proof_max_def = function
    | Done -> -1
    | Step (i, pr, prf) -> max i (max (pr_rule_max_def pr) (proof_max_def prf))
    | Split (i, _, p1, p2) -> max i (max (proof_max_def p1) (proof_max_def p2))
    | Enum (i, p1, _, p2, l) ->
      let m = max (pr_rule_max_def p1) (pr_rule_max_def p2) in
      List.fold_left (fun i prf -> max i (proof_max_def prf)) (max i m) l
    | ExProof (i, j, k, _, _, _, prf) ->
      max (max (max i j) k) (proof_max_def prf)

  (** [pr_rule_def_cut id pr] gives an explicit [id] to cut rules.
      This is because the Coq proof format only accept they as a proof-step *)
  let pr_rule_def_cut m id p =
    let rec pr_rule_def_cut m id = function
      | Annot (_, p) -> pr_rule_def_cut m id p
      | MulC (p, prf) ->
        let bds, m, id', prf' = pr_rule_def_cut m id prf in
        (bds, m, id', MulC (p, prf'))
      | MulPrf (p1, p2) ->
        let bds1, m, id, p1 = pr_rule_def_cut m id p1 in
        let bds2, m, id, p2 = pr_rule_def_cut m id p2 in
        (bds2 @ bds1, m, id, MulPrf (p1, p2))
      | AddPrf (p1, p2) ->
        let bds1, m, id, p1 = pr_rule_def_cut m id p1 in
        let bds2, m, id, p2 = pr_rule_def_cut m id p2 in
        (bds2 @ bds1, m, id, AddPrf (p1, p2))
      | CutPrf p | Gcd (_, p) -> (
        let bds, m, id, p = pr_rule_def_cut m id p in
        try
          let id' = PrfRuleMap.find p m in
          (bds, m, id, Def id')
        with Not_found ->
          let m = PrfRuleMap.add p id m in
          ((id, p) :: bds, m, id + 1, Def id) )
      | (Square _ | Cst _ | Def _ | Hyp _ | Zero) as x -> ([], m, id, x)
    in
    pr_rule_def_cut m id p

  (* Do not define top-level cuts *)
  let pr_rule_def_cut m id = function
    | CutPrf p ->
      let bds, m, ids, p' = pr_rule_def_cut m id p in
      (bds, m, ids, CutPrf p')
    | p -> pr_rule_def_cut m id p

  let rec implicit_cut p = match p with CutPrf p -> implicit_cut p | _ -> p

  let rec pr_rule_collect_defs pr =
    match pr with
    | Annot (_, pr) -> pr_rule_collect_defs pr
    | Def i -> ISet.add i ISet.empty
    | Hyp i -> ISet.empty
    | Cst _ | Zero | Square _ -> ISet.empty
    | MulC (_, pr) | Gcd (_, pr) | CutPrf pr -> pr_rule_collect_defs pr
    | MulPrf (p1, p2) | AddPrf (p1, p2) ->
      ISet.union (pr_rule_collect_defs p1) (pr_rule_collect_defs p2)

  let add_proof x y =
    match (x, y) with Zero, p | p, Zero -> p | _ -> AddPrf (x, y)

  let rec mul_cst_proof c p =
    match p with
    | Annot (s, p) -> Annot (s, mul_cst_proof c p)
    | MulC (v, p') -> MulC (Vect.mul c v, p')
    | _ -> (
      match Q.sign c with
      | 0 -> Zero (* This is likely to be a bug *)
      | -1 ->
        MulC (LinPoly.constant c, p) (* [p] should represent an equality *)
      | 1 -> if Q.one =/ c then p else MulPrf (Cst c, p)
      | _ -> assert false )

  let sMulC v p =
    let c, v' = Vect.decomp_cst v in
    if Vect.is_null v' then mul_cst_proof c p else MulC (v, p)

  let mul_proof p1 p2 =
    match (p1, p2) with
    | Zero, _ | _, Zero -> Zero
    | Cst c, p | p, Cst c -> mul_cst_proof c p
    | _, _ -> MulPrf (p1, p2)

  let prf_rule_of_map m =
    PrfRuleMap.fold (fun k v acc -> add_proof (sMulC v k) acc) m Zero

  let rec dev_prf_rule p =
    match p with
    | Annot (s, p) -> dev_prf_rule p
    | Hyp _ | Def _ | Cst _ | Zero | Square _ ->
      PrfRuleMap.singleton p (LinPoly.constant Q.one)
    | MulC (v, p) ->
      PrfRuleMap.map (fun v1 -> LinPoly.product v v1) (dev_prf_rule p)
    | AddPrf (p1, p2) ->
      PrfRuleMap.merge
        (fun k o1 o2 ->
          match (o1, o2) with
          | None, None -> None
          | None, Some v | Some v, None -> Some v
          | Some v1, Some v2 -> Some (LinPoly.addition v1 v2))
        (dev_prf_rule p1) (dev_prf_rule p2)
    | MulPrf (p1, p2) -> (
      let p1' = dev_prf_rule p1 in
      let p2' = dev_prf_rule p2 in
      let p1'' = prf_rule_of_map p1' in
      let p2'' = prf_rule_of_map p2' in
      match p1'' with
      | Cst c -> PrfRuleMap.map (fun v1 -> Vect.mul c v1) p2'
      | _ -> PrfRuleMap.singleton (MulPrf (p1'', p2'')) (LinPoly.constant Q.one)
      )
    | Gcd (c, p) ->
      PrfRuleMap.singleton
        (Gcd (c, prf_rule_of_map (dev_prf_rule p)))
        (LinPoly.constant Q.one)
    | CutPrf p ->
      PrfRuleMap.singleton
        (CutPrf (prf_rule_of_map (dev_prf_rule p)))
        (LinPoly.constant Q.one)

  let simplify_prf_rule p = prf_rule_of_map (dev_prf_rule p)

  (** [simplify_proof p] removes proof steps that are never re-used. *)
  let rec simplify_proof p =
    match p with
    | Done -> (Done, ISet.empty)
    | Step (i, pr, Done) -> (p, ISet.add i (pr_rule_collect_defs pr))
    | Step (i, pr, prf) ->
      let prf', hyps = simplify_proof prf in
      if not (ISet.mem i hyps) then (prf', hyps)
      else
        ( Step (i, pr, prf')
        , ISet.add i (ISet.union (pr_rule_collect_defs pr) hyps) )
    | Split (i, v, p1, p2) ->
      let p1, h1 = simplify_proof p1 in
      let p2, h2 = simplify_proof p2 in
      if not (ISet.mem i h1) then (p1, h1) (* Should not have computed p2 *)
      else if not (ISet.mem i h2) then (p2, h2)
      else (Split (i, v, p1, p2), ISet.add i (ISet.union h1 h2))
    | Enum (i, p1, v, p2, pl) ->
      let pl, hl = List.split (List.map simplify_proof pl) in
      let hyps = List.fold_left ISet.union ISet.empty hl in
      ( Enum (i, p1, v, p2, pl)
      , ISet.add i
          (ISet.union
             (ISet.union (pr_rule_collect_defs p1) (pr_rule_collect_defs p2))
             hyps) )
    | ExProof (i, j, k, x, z, t, prf) ->
      let prf', hyps = simplify_proof prf in
      if
        (not (ISet.mem i hyps))
        && (not (ISet.mem j hyps))
        && not (ISet.mem k hyps)
      then (prf', hyps)
      else
        ( ExProof (i, j, k, x, z, t, prf')
        , ISet.add i (ISet.add j (ISet.add k hyps)) )

  let rec normalise_proof id prf =
    match prf with
    | Done -> (id, Done)
    | Step (i, Gcd (c, p), Done) -> normalise_proof id (Step (i, p, Done))
    | Step (i, p, prf) ->
      let bds, m, id, p' =
        pr_rule_def_cut PrfRuleMap.empty id (simplify_prf_rule p)
      in
      let id, prf = normalise_proof id prf in
      let prf =
        List.fold_left
          (fun acc (i, p) -> Step (i, CutPrf p, acc))
          (Step (i, p', prf))
          bds
      in
      (id, prf)
    | Split (i, v, p1, p2) ->
      let id, p1 = normalise_proof id p1 in
      let id, p2 = normalise_proof id p2 in
      (id, Split (i, v, p1, p2))
    | ExProof (i, j, k, x, z, t, prf) ->
      let id, prf = normalise_proof id prf in
      (id, ExProof (i, j, k, x, z, t, prf))
    | Enum (i, p1, v, p2, pl) ->
      (* Why do I have  top-level cuts ? *)
      (* let p1 = implicit_cut p1 in
         let p2 = implicit_cut p2 in
         let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in
           (List.fold_left max  0 ids ,
              Enum(i,p1,v,p2,prfs))
      *)
      let bds1, m, id, p1' =
        pr_rule_def_cut PrfRuleMap.empty id (implicit_cut p1)
      in
      let bds2, m, id, p2' = pr_rule_def_cut m id (implicit_cut p2) in
      let ids, prfs = List.split (List.map (normalise_proof id) pl) in
      ( List.fold_left max 0 ids
      , List.fold_left
          (fun acc (i, p) -> Step (i, CutPrf p, acc))
          (Enum (i, p1', v, p2', prfs))
          (bds2 @ bds1) )

  let normalise_proof id prf =
    let prf = fst (simplify_proof prf) in
    let res = normalise_proof id prf in
    if debug then
      Printf.printf "normalise_proof %a -> %a" output_proof prf output_proof
        (snd res);
    res

  (*
  let mul_proof p1 p2 =
    let res = mul_proof p1 p2 in
    Printf.printf "mul_proof %a %a = %a\n"
    output_prf_rule p1 output_prf_rule p2 output_prf_rule res; res

  let add_proof p1 p2 =
    let res = add_proof p1 p2 in
    Printf.printf "add_proof %a %a = %a\n"
    output_prf_rule p1 output_prf_rule p2 output_prf_rule res; res


  let sMulC v p =
    let res = sMulC v p in
    Printf.printf "sMulC %a %a = %a\n" Vect.pp v output_prf_rule p output_prf_rule res ;
    res

  let mul_cst_proof c p  =
    let res = mul_cst_proof c p in
    Printf.printf "mul_cst_proof %s %a = %a\n" (Num.string_of_num c) output_prf_rule p output_prf_rule res ;
    res
 *)

  let proof_of_farkas env vect =
    Vect.fold
      (fun prf x n -> add_proof (mul_cst_proof n (IMap.find x env)) prf)
      Zero vect

  module Env = struct
    let output_hyp_or_def o = function
      | Hyp i -> Printf.fprintf o "Hyp %i" i
      | Def i -> Printf.fprintf o "Def %i" i
      | _ -> ()

    let rec output_hyps o l =
      match l with
      | [] -> ()
      | i :: l -> Printf.fprintf o "%a,%a" output_hyp_or_def i output_hyps l

    let id_of_hyp hyp l =
      let rec xid_of_hyp i l' =
        match l' with
        | [] ->
          Printf.fprintf stdout "\nid_of_hyp: %a notin [%a]\n" output_hyp_or_def
            hyp output_hyps l;
          failwith "Cannot find hyp or def"
        | hyp' :: l' -> if hyp = hyp' then i else xid_of_hyp (i + 1) l'
      in
      xid_of_hyp 0 l
  end

  let cmpl_prf_rule norm (cst : Q.t -> 'a) env prf =
    let rec cmpl = function
      | Annot (s, p) -> cmpl p
      | (Hyp _ | Def _) as h -> Mc.PsatzIn (CamlToCoq.nat (Env.id_of_hyp h env))
      | Cst i -> Mc.PsatzC (cst i)
      | Zero -> Mc.PsatzZ
      | MulPrf (p1, p2) -> Mc.PsatzMulE (cmpl p1, cmpl p2)
      | AddPrf (p1, p2) -> Mc.PsatzAdd (cmpl p1, cmpl p2)
      | MulC (lp, p) ->
        let lp = norm (LinPoly.coq_poly_of_linpol cst lp) in
        Mc.PsatzMulC (lp, cmpl p)
      | Square lp -> Mc.PsatzSquare (norm (LinPoly.coq_poly_of_linpol cst lp))
      | _ -> failwith "Cuts should already be compiled"
    in
    cmpl prf

  let cmpl_prf_rule_z env r =
    cmpl_prf_rule Mc.normZ (fun x -> CamlToCoq.bigint (Q.num x)) env r

  let cmpl_pol_z lp =
    try
      let cst x = CamlToCoq.bigint (Q.num x) in
      Mc.normZ (LinPoly.coq_poly_of_linpol cst lp)
    with x ->
      Printf.printf "cmpl_pol_z %s %a\n" (Printexc.to_string x) LinPoly.pp lp;
      raise x

  let rec cmpl_proof env prf =
    match prf with
    | Done -> Mc.DoneProof
    | Step (i, p, prf) -> (
      match p with
      | CutPrf p' ->
        Mc.CutProof (cmpl_prf_rule_z env p', cmpl_proof (Def i :: env) prf)
      | _ -> Mc.RatProof (cmpl_prf_rule_z env p, cmpl_proof (Def i :: env) prf)
      )
    | Split (i, v, p1, p2) ->
      Mc.SplitProof
        ( cmpl_pol_z v
        , cmpl_proof (Def i :: env) p1
        , cmpl_proof (Def i :: env) p2 )
    | Enum (i, p1, _, p2, l) ->
      Mc.EnumProof
        ( cmpl_prf_rule_z env p1
        , cmpl_prf_rule_z env p2
        , List.map (cmpl_proof (Def i :: env)) l )
    | ExProof (i, j, k, x, _, _, prf) ->
      Mc.ExProof
        (CamlToCoq.positive x, cmpl_proof (Def i :: Def j :: Def k :: env) prf)

  let compile_proof env prf =
    let id = 1 + proof_max_def prf in
    let _, prf = normalise_proof id prf in
    cmpl_proof (List.map (fun i -> Hyp i) env) prf

  let rec eval_prf_rule env = function
    | Annot (s, p) -> eval_prf_rule env p
    | Hyp i | Def i -> env i
    | Cst n -> (
      ( Vect.set 0 n Vect.null
      , match Q.compare n Q.zero with
        | 0 -> Ge
        | 1 -> Gt
        | _ -> failwith "eval_prf_rule : negative constant" ) )
    | Zero -> (Vect.null, Ge)
    | Square v -> (LinPoly.product v v, Ge)
    | MulC (v, p) -> (
      let p1, o = eval_prf_rule env p in
      match o with
      | Eq -> (LinPoly.product v p1, Eq)
      | _ ->
        Printf.fprintf stdout "MulC(%a,%a) invalid 2d arg %a %s" Vect.pp v
          output_prf_rule p Vect.pp p1 (string_of_op o);
        failwith "eval_prf_rule : not an equality" )
    | Gcd (g, p) ->
      let v, op = eval_prf_rule env p in
      (Vect.div (Q.of_bigint g) v, op)
    | MulPrf (p1, p2) ->
      let v1, o1 = eval_prf_rule env p1 in
      let v2, o2 = eval_prf_rule env p2 in
      (LinPoly.product v1 v2, opMult o1 o2)
    | AddPrf (p1, p2) ->
      let v1, o1 = eval_prf_rule env p1 in
      let v2, o2 = eval_prf_rule env p2 in
      (LinPoly.addition v1 v2, opAdd o1 o2)
    | CutPrf p -> eval_prf_rule env p

  let is_unsat (p, o) =
    let c, r = Vect.decomp_cst p in
    if Vect.is_null r then not (eval_op o c Q.zero) else false

  let rec eval_proof env p =
    match p with
    | Done -> failwith "Proof is not finished"
    | Step (i, prf, rst) ->
      let p, o = eval_prf_rule (fun i -> IMap.find i env) prf in
      if is_unsat (p, o) then true
      else if rst = Done then begin
        Printf.fprintf stdout "Last inference %a %s\n" LinPoly.pp p
          (string_of_op o);
        false
      end
      else eval_proof (IMap.add i (p, o) env) rst
    | Split (i, v, p1, p2) -> failwith "Not implemented"
    | Enum (i, r1, v, r2, l) ->
      let _ = eval_prf_rule (fun i -> IMap.find i env) r1 in
      let _ = eval_prf_rule (fun i -> IMap.find i env) r2 in
      (* Should check bounds *)
      failwith "Not implemented"
    | ExProof _ -> failwith "Not implemented"
end

module WithProof = struct
  type t = (LinPoly.t * op) * ProofFormat.prf_rule

  (* The comparison ignores proofs on purpose *)
  let compare : t -> t -> int =
   fun ((lp1, o1), _) ((lp2, o2), _) ->
    let c = Vect.compare lp1 lp2 in
    if c = 0 then compare_op o1 o2 else c

  let annot s (p, prf) = (p, ProofFormat.Annot (s, prf))

  let output o ((lp, op), prf) =
    Printf.fprintf o "%a %s 0 by %a\n" LinPoly.pp lp (string_of_op op)
      ProofFormat.output_prf_rule prf

  let output_sys o l = List.iter (Printf.fprintf o "%a\n" output) l

  exception InvalidProof

  let zero = ((Vect.null, Eq), ProofFormat.Zero)
  let const n = ((LinPoly.constant n, Ge), ProofFormat.Cst n)
  let of_cstr (c, prf) = ((Vect.set 0 (Q.neg c.cst) c.coeffs, c.op), prf)

  let product : t -> t -> t =
   fun ((p1, o1), prf1) ((p2, o2), prf2) ->
    ((LinPoly.product p1 p2, opMult o1 o2), ProofFormat.mul_proof prf1 prf2)

  let addition : t -> t -> t =
   fun ((p1, o1), prf1) ((p2, o2), prf2) ->
    ((Vect.add p1 p2, opAdd o1 o2), ProofFormat.add_proof prf1 prf2)

  let neg : t -> t =
   fun ((p1, o1), prf1) ->
    match o1 with
    | Eq ->
      ((Vect.mul Q.minus_one p1, o1), ProofFormat.mul_cst_proof Q.minus_one prf1)
    | _ -> failwith "neg: invalid proof"

  let mult p ((p1, o1), prf1) =
    match o1 with
    | Eq -> ((LinPoly.product p p1, o1), ProofFormat.sMulC p prf1)
    | Gt | Ge ->
      let n, r = Vect.decomp_cst p in
      if Vect.is_null r && n >/ Q.zero then
        ((LinPoly.product p p1, o1), ProofFormat.mul_cst_proof n prf1)
      else (
        if debug then
          Printf.printf "mult_error %a [*] %a\n" LinPoly.pp p output
            ((p1, o1), prf1);
        raise InvalidProof )

  let cutting_plane ((p, o), prf) =
    let c, p' = Vect.decomp_cst p in
    let g = Vect.gcd p' in
    if Z.equal Z.one g || c =/ Q.zero || not (Z.equal (Q.den c) Z.one) then None
      (* Nothing to do *)
    else
      let c1 = c // Q.of_bigint g in
      let c1' = Q.floor c1 in
      if c1 =/ c1' then None
      else
        match o with
        | Eq ->
          Some ((Vect.set 0 Q.minus_one Vect.null, Eq), ProofFormat.CutPrf prf)
        | Gt -> failwith "cutting_plane ignore strict constraints"
        | Ge ->
          (* This is a non-trivial common divisor *)
          Some
            ( (Vect.set 0 c1' (Vect.div (Q.of_bigint g) p), o)
            , ProofFormat.CutPrf prf )

  let construct_sign p =
    let c, p' = Vect.decomp_cst p in
    if Vect.is_null p' then
      Some
        ( match Q.sign c with
        | 0 -> (true, Eq, ProofFormat.Zero)
        | 1 -> (true, Gt, ProofFormat.Cst c)
        | _ (*-1*) -> (false, Gt, ProofFormat.Cst (Q.neg c)) )
    else None

  let get_sign l p =
    match construct_sign p with
    | None -> (
      try
        let (p', o), prf =
          List.find (fun ((p', o), prf) -> Vect.equal p p') l
        in
        Some (true, o, prf)
      with Not_found -> (
        let p = Vect.uminus p in
        try
          let (p', o), prf =
            List.find (fun ((p', o), prf) -> Vect.equal p p') l
          in
          Some (false, o, prf)
        with Not_found -> None ) )
    | Some s -> Some s

  let mult_sign : bool -> t -> t =
   fun b ((p, o), prf) -> if b then ((p, o), prf) else ((Vect.uminus p, o), prf)

  let rec linear_pivot sys ((lp1, op1), prf1) x ((lp2, op2), prf2) =
    (* lp1 = a1.x + b1 *)
    let a1, b1 = LinPoly.factorise x lp1 in
    (* lp2 = a2.x + b2 *)
    let a2, b2 = LinPoly.factorise x lp2 in
    if Vect.is_null a2 then (* We are done *)
      Some ((lp2, op2), prf2)
    else
      match (op1, op2) with
      | Eq, (Ge | Gt) -> (
        match get_sign sys a1 with
        | None -> None (* Impossible to pivot without sign information *)
        | Some (b, o, prf) ->
          let sa1 = mult_sign b ((a1, o), prf) in
          let sa2 = if b then Vect.uminus a2 else a2 in
          let (lp2, op2), prf2 =
            addition
              (product sa1 ((lp2, op2), prf2))
              (mult sa2 ((lp1, op1), prf1))
          in
          linear_pivot sys ((lp1, op1), prf1) x ((lp2, op2), prf2) )
      | Eq, Eq ->
        let (lp2, op2), prf2 =
          addition
            (mult a1 ((lp2, op2), prf2))
            (mult (Vect.uminus a2) ((lp1, op1), prf1))
        in
        linear_pivot sys ((lp1, op1), prf1) x ((lp2, op2), prf2)
      | (Ge | Gt), (Ge | Gt) -> (
        match (get_sign sys a1, get_sign sys a2) with
        | Some (b1, o1, p1), Some (b2, o2, p2) ->
          if b1 <> b2 then
            let (lp2, op2), prf2 =
              addition
                (product (mult_sign b1 ((a1, o1), p1)) ((lp2, op2), prf2))
                (product (mult_sign b2 ((a2, o2), p2)) ((lp1, op1), prf1))
            in
            linear_pivot sys ((lp1, op1), prf1) x ((lp2, op2), prf2)
          else None
        | _ -> None )
      | (Ge | Gt), Eq -> failwith "pivot: equality as second argument"

  let linear_pivot sys ((lp1, op1), prf1) x ((lp2, op2), prf2) =
    match linear_pivot sys ((lp1, op1), prf1) x ((lp2, op2), prf2) with
    | None -> None
    | Some (c, p) -> Some (c, ProofFormat.simplify_prf_rule p)

  let is_substitution strict ((p, o), prf) =
    let pred v = if strict then v =/ Q.one || v =/ Q.minus_one else true in
    match o with Eq -> LinPoly.search_linear pred p | _ -> None

  let subst1 sys0 =
    let oeq, sys' = extract (is_substitution true) sys0 in
    match oeq with
    | None -> sys0
    | Some (v, pc) -> (
      match simplify (linear_pivot sys0 pc v) sys' with
      | None -> sys0
      | Some sys' -> sys' )

  let sort (sys : t list) =
    let size ((p, o), prf) =
      let _, p' = Vect.decomp_cst p in
      let (x, q), p' = Vect.decomp_fst p' in
      Vect.fold
        (fun (l, (q, x)) x' q' ->
          let q' = Q.abs q' in
          (l + 1, if q </ q then (q, x) else (q', x')))
        (1, (Q.abs q, x))
        p
    in
    let cmp ((l1, (q1, _)), ((_, o), _)) ((l2, (q2, _)), ((_, o'), _)) =
      if l1 < l2 then -1 else if l1 = l2 then Q.compare q1 q2 else 1
    in
    List.sort cmp (List.rev_map (fun wp -> (size wp, wp)) sys)

  let iterate_pivot p sys0 =
    let elim sys =
      let oeq, sys' = extract p sys in
      match oeq with
      | None -> None
      | Some (v, pc) -> simplify (linear_pivot sys0 pc v) sys'
    in
    iterate_until_stable elim (List.map snd (sort sys0))

  let subst_constant is_int sys =
    let is_integer q = Q.(q =/ floor q) in
    let is_constant ((c, o), p) =
      match o with
      | Ge | Gt -> None
      | Eq -> (
        Vect.Bound.(
          match of_vect c with
          | None -> None
          | Some b ->
            if (not is_int) || is_integer (b.cst // b.coeff) then
              Monomial.get_var (LinPoly.MonT.retrieve b.var)
            else None) )
    in
    iterate_pivot is_constant sys

  let subst sys0 = iterate_pivot (is_substitution true) sys0

  let saturate_subst b sys0 =
    let select = is_substitution b in
    let gen (v, pc) ((c, op), prf) =
      if ISet.mem v (LinPoly.variables c) then
        linear_pivot sys0 pc v ((c, op), prf)
      else None
    in
    saturate select gen sys0

  let simple_pivot (q1, x) ((v1, o1), prf1) ((v2, o2), prf2) =
    let q2 = Vect.get x v2 in
    if q2 =/ Q.zero then None
    else
      let cv1, cv2 =
        if Q.sign q1 <> Q.sign q2 then (Q.abs q2, Q.abs q1)
        else
          match (o1, o2) with
          | Eq, _ -> (q2, Q.abs q1)
          | _, Eq -> (Q.abs q2, q2)
          | _, _ -> (Q.zero, Q.zero)
      in
      if cv2 =/ Q.zero then None
      else
        Some
          ( (Vect.mul_add cv1 v1 cv2 v2, opAdd o1 o2)
          , ProofFormat.add_proof
              (ProofFormat.mul_cst_proof cv1 prf1)
              (ProofFormat.mul_cst_proof cv2 prf2) )

  open Vect.Bound

  let mul_bound w1 w2 =
    let (p1, o1), prf1 = w1 in
    let (p2, o2), prf2 = w2 in
    match (LinPoly.get_bound p1, LinPoly.get_bound p2) with
    | None, _ | _, None -> None
    | ( Some {cst = c1; var = v1; coeff = c1'}
      , Some {cst = c2; var = v2; coeff = c2'} ) -> (
      let good_coeff b o =
        match o with
        | Eq -> Some (Q.neg b)
        | _ -> if b <=/ Q.zero then Some (Q.neg b) else None
      in
      match (good_coeff c1 o2, good_coeff c2 o1) with
      | None, _ | _, None -> None
      | Some c1, Some c2 ->
        let ext_mult c w =
          if c =/ Q.zero then zero else mult (LinPoly.constant c) w
        in
        Some
          (addition
             (addition (product w1 w2) (ext_mult c1 w2))
             (ext_mult c2 w1)) )
end

(* Local Variables: *)
(* coding: utf-8 *)
(* End: *)