Refactoring of Micromega code using a Simplex linear solver

- Simplex based linear prover
  Unset Simplex to get Fourier elimination
  For lia and nia, do not enumerate but generate cutting planes.
- Better non-linear support
  Factorisation of the non-linear pre-processing
  Careful handling of equation x=e, x is only eliminated if x is used linearly
- More opaque interfaces
  (Linear solvers Simplex and Mfourier are independent)
- Set Dump Arith "file" so that lia,nia calls generate Coq goals
  in filexxx.v. Used to collect benchmarks and regressions.
- Rationalise the test-suite
  example.v only tests psatz Z
  example_nia.v only tests lia, nia
  In both files, the tests are in essence the same.
  In particular, if a test is solved by psatz but not by nia,
  we finish the goal by an explicit Abort.
  There are additional tests in example_nia.v which require specific
  integer reasoning out of scope of psatz.

1 files changed, 269 insertions, 51 deletions

diff --git a/plugins/micromega/polynomial.mli b/plugins/micromega/polynomial.mli
index 4c095202ab..6f26f7a959 100644
--- a/plugins/micromega/polynomial.mli
+++ b/plugins/micromega/polynomial.mli
@@ -8,111 +8,329 @@
 (*         *     (see LICENSE file for the text of the license)         *)
 (************************************************************************)
 
+open Mutils
+
+module Mc = Micromega
+
+val max_nb_cstr : int ref
+
 type var = int
 
 module Monomial : sig
-
+  (** A monomial is represented by a multiset of variables  *)
   type t
+
+   (** [fold f m acc]
+       folds over the variables with multiplicities *)
   val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a
+
+  (** [const]
+      @return the empty monomial i.e. without any variable *)
   val const : t
+
+  (** [var x]
+      @return the monomial x^1 *)
+  val var : var -> t
+
+  (** [sqrt m]
+      @return [Some r] iff r^2 = m *)
   val sqrt : t -> t option
+
+  (** [is_var m]
+      @return [true] iff m = x^1 for some variable x *)
   val is_var : t -> bool
+
+  (** [div m1 m2]
+      @return a pair [mr,n] such that mr * (m2)^n = m1 where n is maximum *)
   val div : t -> t -> t * int
 
+  (** [compare m1 m2] provides a total order over monomials*)
   val compare : t -> t -> int
 
+  (** [variables m]
+      @return the set of variables with (strictly) positive multiplicities *)
+  val variables : t -> ISet.t
+end
+
+module MonMap : sig
+  include Map.S with type key = Monomial.t
+
+  val union : (Monomial.t -> 'a -> 'a -> 'a option) -> 'a t -> 'a t -> 'a t
 end
 
 module Poly : sig
+  (** Representation of polonomial with rational coefficient.
+      a1.m1 + ... + c where
+      - ai are rational constants (num type)
+      - mi are monomials
+      - c is a rational constant
+
+   *)
 
   type t
 
+  (** [constant c]
+      @return the constant polynomial c *)
   val constant : Num.num -> t
+
+  (** [variable x]
+      @return the polynomial 1.x^1 *)
   val variable : var -> t
+
+  (** [addition p1 p2]
+      @return the polynomial p1+p2 *)
   val addition : t -> t -> t
+
+  (** [product p1 p2]
+      @return the polynomial p1*p2 *)
   val product : t -> t -> t
+
+  (** [uminus p]
+      @return the polynomial -p i.e product by -1 *)
   val uminus : t -> t
+
+  (** [get mi p]
+      @return the coefficient ai of the  monomial mi. *)
   val get : Monomial.t -> t -> Num.num
+
+
+  (** [fold f p a] folds f over the monomials of p with non-zero coefficient *)
   val fold : (Monomial.t -> Num.num -> 'a -> 'a) -> t -> 'a -> 'a
 
+  (** [is_linear p]
+      @return true if the polynomial is of the form a1.x1 +...+ an.xn + c
+      i.e every monomial is made of at most a variable *)
   val is_linear : t -> bool
 
+
+  (** [add m n p]
+      @return the polynomial n*m + p *)
   val add : Monomial.t -> Num.num -> t -> t
 
+  (** [variables p]
+      @return the set of variables of the polynomial p *)
+  val variables : t -> ISet.t
+
 end
 
-module Vect : sig
+type cstr = {coeffs : Vect.t ; op : op ; cst : Num.num} (** Representation of linear constraints *)
+and op = Eq | Ge | Gt
 
-  type var = int
-  type t = (var * Num.num) list
-  val hash : t -> int
-  val equal : t -> t -> bool
-  val compare : t -> t -> int
-  val pp_vect : 'a -> t -> unit
+val eval_op : op -> Num.num -> Num.num -> bool
+
+(*val opMult : op -> op -> op*)
+
+val opAdd  : op -> op -> op
+
+(** [is_strict c]
+    @return whether the constraint is strict i.e. c.op = Gt *)
+val is_strict : cstr -> bool
+
+exception Strict
+
+module LinPoly : sig
+  (** Linear(ised) polynomials represented as a [Vect.t]
+      i.e a sorted association list.
+      The constant is the coefficient of the variable 0
+
+      Each linear polynomial can be interpreted as a multi-variate polynomial.
+      There is a bijection mapping between a linear variable and a monomial
+      (see module [MonT])
+   *)
+
+  type t = Vect.t
+
+  (** Each variable of a linear polynomial is mapped to a monomial.
+      This is done using the monomial tables of the module MonT. *)
+
+  module MonT : sig
+    (** [clear ()] clears the mapping. *)
+    val clear : unit -> unit
+
+    (** [retrieve x]
+        @return the monomial corresponding to the variable [x] *)
+    val retrieve : int -> Monomial.t
+
+  end
+
+  (** [linpol_of_pol p] linearise the polynomial p *)
+  val linpol_of_pol : Poly.t -> t
+
+  (** [var x]
+      @return 1.y where y is the variable index of the monomial x^1.
+   *)
+  val var : var -> t
 
-  val get : var -> t -> Num.num option
-  val set : var -> Num.num -> t -> t
-  val fresh : (int * 'a) list -> int
-  val update : Int.t -> (Num.num -> Num.num) ->
-    (Int.t * Num.num) list -> (Int.t * Num.num) list
-  val null : t
+  (** [coq_poly_of_linpol c p]
+      @param p is a multi-variate polynomial.
+      @param c maps a rational to a Coq polynomial coefficient.
+      @return the coq expression corresponding to polynomial [p].*)
+  val coq_poly_of_linpol : (Num.num -> 'a) -> t -> 'a Mc.pExpr
 
-  val from_list : Num.num list -> t
-  val to_list : t -> Num.num list
+  (** [of_monomial m]
+      @returns 1.x where x is the variable (index) for monomial m *)
+  val of_monomial : Monomial.t -> t
 
-  val add : t -> t -> t
-  val mul : Num.num -> t -> t
+  (** [variables p]
+      @return the set of variables of the polynomial p
+      interpreted as a multi-variate polynomial *)
+  val variables : t -> ISet.t
+
+  (** [is_linear p]
+      @return whether the multi-variate polynomial is linear. *)
+  val is_linear : t -> bool
+
+  (** [is_linear_for x p]
+      @return true if the polynomial is linear in x
+      i.e can be written c*x+r where c is a constant and r is independent from x *)
+  val is_linear_for : var -> t -> bool
+
+  (** [constant c]
+      @return the constant polynomial c
+   *)
+  val constant : Num.num -> t
+
+  (** [search_linear pred p]
+      @return a variable x such p = a.x + b such that
+      p is linear in x i.e x does not occur in b and
+      a is a constant such that [pred a] *)
+
+  val search_linear : (Num.num -> bool) -> t -> var option
+
+  (** [search_all_linear pred p]
+      @return all the variables x such p = a.x + b such that
+      p is linear in x i.e x does not occur in b and
+      a is a constant such that [pred a] *)
+  val search_all_linear : (Num.num -> bool) -> t -> var list
+
+ (** [product p q]
+     @return the product of the polynomial [p*q] *)
+  val product : t -> t -> t
+
+  (** [factorise x p]
+      @return [a,b] such that [p = a.x + b]
+      and [x] does not occur in [b] *)
+  val factorise : var -> t -> t * t
+
+  (** [collect_square p]
+      @return a mapping m such that m[s] = s^2
+      for every s^2 that is a monomial of [p] *)
+  val collect_square : t -> Monomial.t MonMap.t
+
+
+  (** [pp_var o v] pretty-prints a monomial indexed by v. *)
+  val pp_var : out_channel -> var -> unit
+
+  (** [pp o p] pretty-prints a polynomial. *)
+  val pp : out_channel -> t -> unit
+
+  (** [pp_goal typ o l] pretty-prints the list of constraints as a Coq goal. *)
+  val pp_goal : string -> out_channel -> (t * op) list -> unit
 
 end
 
-type cstr_compat = {coeffs : Vect.t ; op : op ; cst : Num.num}
-and op = Eq | Ge
+module ProofFormat : sig
+  (** Proof format used by the proof-generating procedures.
+      It is fairly close to Coq format but a bit more liberal.
 
-type prf_rule =
-  | Hyp of int
-  | Def of int
-  | Cst  of Big_int.big_int
-  | Zero
-  | Square of (Vect.t * Num.num)
-  | MulC of (Vect.t * Num.num) * prf_rule
-  | Gcd of Big_int.big_int * prf_rule
-  | MulPrf of prf_rule * prf_rule
-  | AddPrf of prf_rule * prf_rule
-  | CutPrf of prf_rule
+      It is used for proofs over Z, Q, R.
+      However, certain constructions e.g. [CutPrf] are only relevant for Z.
+   *)
 
-type proof =
-  | Done
-  | Step of int * prf_rule * proof
-  | Enum of int * prf_rule * Vect.t * prf_rule * proof list
+  type prf_rule =
+    | Annot of string * prf_rule
+    | Hyp of int
+    | Def of int
+    | Cst  of Num.num
+    | Zero
+    | Square of Vect.t
+    | MulC of Vect.t * prf_rule
+    | Gcd of Big_int.big_int * prf_rule
+    | MulPrf of prf_rule * prf_rule
+    | AddPrf of prf_rule * prf_rule
+    | CutPrf of prf_rule
 
-val proof_max_id : proof -> int
+  type proof =
+    | Done
+    | Step of int * prf_rule * proof
+    | Enum of int * prf_rule * Vect.t * prf_rule * proof list
 
-val normalise_proof : int -> proof -> int * proof
+  val pr_rule_max_id : prf_rule -> int
 
-val output_proof : out_channel -> proof -> unit
+  val proof_max_id : proof -> int
 
-val add_proof : prf_rule -> prf_rule -> prf_rule
-val mul_proof : Big_int.big_int -> prf_rule -> prf_rule
+  val normalise_proof : int -> proof -> int * proof
 
-module LinPoly : sig
+  val output_prf_rule : out_channel -> prf_rule -> unit
 
-  type t = Vect.t * Num.num
+  val output_proof : out_channel -> proof -> unit
 
-  module MonT : sig
+  val add_proof : prf_rule -> prf_rule -> prf_rule
 
-    val clear : unit -> unit
-    val retrieve : int -> Monomial.t
+  val mul_cst_proof : Num.num -> prf_rule -> prf_rule
 
-  end
+  val mul_proof : prf_rule -> prf_rule -> prf_rule
 
-  val pivot_eq : Vect.var ->
-           cstr_compat * prf_rule ->
-           cstr_compat * prf_rule -> (cstr_compat * prf_rule) option
+  val compile_proof : int list -> proof -> Micromega.zArithProof
 
-  val linpol_of_pol : Poly.t -> t
+  val cmpl_prf_rule : ('a Micromega.pExpr -> 'a Micromega.pol) ->
+                      (Num.num -> 'a) -> (int list) -> prf_rule -> 'a Micromega.psatz
+
+  val proof_of_farkas : prf_rule IMap.t -> Vect.t -> prf_rule
+
+  val eval_prf_rule : (int -> LinPoly.t * op) -> prf_rule -> LinPoly.t * op
+
+  val eval_proof : (LinPoly.t * op) IMap.t -> proof -> bool
 
 end
 
-val output_cstr : out_channel -> cstr_compat -> unit
+val output_cstr : out_channel -> cstr -> unit
+
+val output_nlin_cstr : out_channel -> cstr -> unit
 
 val opMult : op -> op -> op
+
+(** [module WithProof] constructs polynomials packed with the proof that their sign is correct. *)
+module WithProof :
+sig
+
+  type t = (LinPoly.t * op) * ProofFormat.prf_rule
+
+  (** [InvalidProof] is raised if the operation is invalid. *)
+  exception InvalidProof
+
+  val annot : string -> t -> t
+
+  val of_cstr : cstr * ProofFormat.prf_rule -> t
+
+  (** [out_channel chan c] pretty-prints the constraint [c] over the channel [chan] *)
+  val output : out_channel -> t -> unit
+
+  (** [zero] represents the tautology (0=0) *)
+  val zero : t
+
+  (** [product p q]
+      @return the polynomial p*q with its sign and proof *)
+  val product : t -> t -> t
+
+  (** [addition p q]
+      @return the polynomial p+q with its sign and proof *)
+  val addition : t -> t -> t
+
+  (** [mult p q]
+      @return the polynomial p*q with its sign and proof.
+      @raise InvalidProof if p is not a constant and p  is not an equality *)
+  val mult : LinPoly.t -> t -> t
+
+  (** [cutting_plane p] does integer reasoning and adjust the constant to be integral *)
+  val cutting_plane : t -> t option
+
+  (** [linear_pivot sys p x q]
+      @return the polynomial [q] where [x] is eliminated using the polynomial [p]
+      The pivoting operation is only defined if
+      - p is linear in x i.e p = a.x+b and x neither occurs in a and b
+      - The pivoting also requires some sign conditions for [a]
+   *)
+  val linear_pivot : t list -> t -> Vect.var -> t -> t option
+
+end