[micromega] fix of bug #11191

- Add an instance to ZifyInst to instruct zify that
 0 < x -> 0 < y -> 0 < Z.pow x y

- More aggressive interval analysis to bound non-linear monomials.

6 files changed, 103 insertions, 43 deletions

diff --git a/plugins/micromega/ZifyInst.v b/plugins/micromega/ZifyInst.v
index 97f6fe0613..edfb5a2a94 100644
--- a/plugins/micromega/ZifyInst.v
+++ b/plugins/micromega/ZifyInst.v
@@ -523,3 +523,22 @@ Instance SatProdPos : Saturate Z.mul :=
     SatOk := Z.mul_pos_pos
   |}.
 Add Saturate SatProdPos.
+
+Lemma pow_pos_strict :
+  forall a b,
+    0 < a -> 0 < b -> 0 < a ^ b.
+Proof.
+  intros.
+  apply Z.pow_pos_nonneg; auto.
+  apply Z.lt_le_incl;auto.
+Qed.
+
+
+Instance SatPowPos : Saturate Z.pow :=
+  {|
+    PArg1 := fun x => 0 < x;
+    PArg2 := fun y => 0 < y;
+    PRes  := fun r => 0 < r;
+    SatOk := pow_pos_strict
+  |}.
+Add Saturate SatPowPos.
diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml
index cb15274736..61234145e1 100644
--- a/plugins/micromega/certificate.ml
+++ b/plugins/micromega/certificate.ml
@@ -395,50 +395,40 @@ let saturate_by_linear_equalities sys =
       output_sys sys output_sys sys';
   sys'
 
-(* let saturate_linear_equality_non_linear sys0 =
-  let (l,_) = extract_all (is_substitution false) sys0 in
-  let rec elim l acc =
-    match l with
-    | [] -> acc
-    | (v,pc)::l' ->
-       let nc = saturate (non_linear_pivot sys0 pc v) (sys0@acc) in
-       elim l' (nc@acc) in
-  elim l []
- *)
-
-let bounded_vars (sys : WithProof.t list) =
-  let l = fst (extract_all (fun ((p, o), prf) -> LinPoly.is_variable p) sys) in
-  List.fold_left (fun acc (i, wp) -> IMap.add i wp acc) IMap.empty l
-
-let rec power n p = if n = 1 then p else WithProof.product p (power (n - 1) p)
-
-let bound_monomial mp m =
-  if Monomial.is_var m || Monomial.is_const m then None
-  else
-    try
-      Some
-        (Monomial.fold
-           (fun v i acc ->
-             let wp = IMap.find v mp in
-             WithProof.product (power i wp) acc)
-           m (WithProof.const (Int 1)))
-    with Not_found -> None
-
 let bound_monomials (sys : WithProof.t list) =
-  let mp = bounded_vars sys in
-  let m =
+  let l =
+    extract_all
+      (fun ((p, o), _) ->
+        match LinPoly.get_bound p with
+        | None -> None
+        | Some Vect.Bound.{cst; var; coeff} ->
+          Some (Monomial.degree (LinPoly.MonT.retrieve var)))
+      sys
+  in
+  let deg =
+    List.fold_left (fun acc ((p, o), _) -> max acc (LinPoly.degree p)) 0 sys
+  in
+  let vars =
     List.fold_left
-      (fun acc ((p, _), _) ->
-        Vect.fold
-          (fun acc v _ ->
-            let m = LinPoly.MonT.retrieve v in
-            match bound_monomial mp m with
-            | None -> acc
-            | Some r -> IMap.add v r acc)
-          acc p)
-      IMap.empty sys
+      (fun acc ((p, o), _) -> ISet.union (LinPoly.monomials p) acc)
+      ISet.empty sys
+  in
+  let bounds =
+    saturate_bin
+      (fun (i1, w1) (i2, w2) ->
+        if i1 + i2 > deg then None
+        else
+          match WithProof.mul_bound w1 w2 with
+          | None -> None
+          | Some b -> Some (i1 + i2, b))
+      (fst l)
+  in
+  let has_mon (_, ((p, o), _)) =
+    match LinPoly.get_bound p with
+    | None -> false
+    | Some Vect.Bound.{cst; var; coeff} -> ISet.mem var vars
   in
-  IMap.fold (fun _ e acc -> e :: acc) m []
+  List.map snd (List.filter has_mon bounds) @ snd l
 
 let develop_constraints prfdepth n_spec sys =
   LinPoly.MonT.clear ();
diff --git a/plugins/micromega/mutils.ml b/plugins/micromega/mutils.ml
index 03f042647c..160b492d3d 100644
--- a/plugins/micromega/mutils.ml
+++ b/plugins/micromega/mutils.ml
@@ -140,6 +140,25 @@ let saturate p f sys =
     Printexc.print_backtrace stdout;
     raise x
 
+let saturate_bin (f : 'a -> 'a -> 'a option) (l : 'a list) =
+  let rec map_with acc e l =
+    match l with
+    | [] -> acc
+    | e' :: l' -> (
+      match f e e' with
+      | None -> map_with acc e l'
+      | Some r -> map_with (r :: acc) e l' )
+  in
+  let rec map2_with acc l' =
+    match l' with [] -> acc | e' :: l' -> map2_with (map_with acc e' l) l'
+  in
+  let rec iterate acc l' =
+    match map2_with [] l' with
+    | [] -> List.rev_append l' acc
+    | res -> iterate (List.rev_append l' acc) res
+  in
+  iterate [] l
+
 open Num
 open Big_int
 
diff --git a/plugins/micromega/mutils.mli b/plugins/micromega/mutils.mli
index ef8d154b13..5dcaf3be44 100644
--- a/plugins/micromega/mutils.mli
+++ b/plugins/micromega/mutils.mli
@@ -116,6 +116,7 @@ val simplify : ('a -> 'a option) -> 'a list -> 'a list option
 val saturate :
   ('a -> 'b option) -> ('b * 'a -> 'a -> 'a option) -> 'a list -> 'a list
 
+val saturate_bin : ('a -> 'a -> 'a option) -> 'a list -> 'a list
 val generate : ('a -> 'b option) -> 'a list -> 'b list
 val app_funs : ('a -> 'b option) list -> 'a -> 'b option
 val command : string -> string array -> 'a -> 'b
diff --git a/plugins/micromega/polynomial.ml b/plugins/micromega/polynomial.ml
index a4f9b60b14..b20213979b 100644
--- a/plugins/micromega/polynomial.ml
+++ b/plugins/micromega/polynomial.ml
@@ -379,6 +379,8 @@ module LinPoly = struct
         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)
 
@@ -892,8 +894,9 @@ module WithProof = struct
       if Vect.is_null r && n >/ Int 0 then
         ((LinPoly.product p p1, o1), ProofFormat.mul_cst_proof n prf1)
       else (
-        Printf.printf "mult_error %a [*] %a\n" LinPoly.pp p output
-          ((p1, o1), prf1);
+        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) =
@@ -1027,6 +1030,31 @@ module WithProof = struct
       else None
     in
     saturate select gen sys0
+
+  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 (minus_num b)
+        | _ -> if b <=/ Int 0 then Some (minus_num 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 =/ Int 0 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: *)
diff --git a/plugins/micromega/polynomial.mli b/plugins/micromega/polynomial.mli
index 7e905ac69b..4b56b037e0 100644
--- a/plugins/micromega/polynomial.mli
+++ b/plugins/micromega/polynomial.mli
@@ -224,6 +224,8 @@ module LinPoly : sig
       p is linear in x i.e x does not occur in b and
       a is a constant such that [pred a] *)
 
+  val get_bound : t -> Vect.Bound.t option
+
   val product : t -> t -> t
   (** [product p q]
      @return the product of the polynomial [p*q] *)
@@ -372,4 +374,5 @@ module WithProof : sig
 
   val saturate_subst : bool -> t list -> t list
   val is_substitution : bool -> t -> var option
+  val mul_bound : t -> t -> t option
 end