aboutsummaryrefslogtreecommitdiff
path: root/plugins
diff options
context:
space:
mode:
Diffstat (limited to 'plugins')
-rw-r--r--plugins/micromega/Lia.v2
-rw-r--r--plugins/micromega/MExtraction.v11
-rw-r--r--plugins/micromega/QMicromega.v2
-rw-r--r--plugins/micromega/ZMicromega.v14
-rw-r--r--plugins/micromega/certificate.ml1718
-rw-r--r--plugins/micromega/certificate.mli28
-rw-r--r--plugins/micromega/coq_micromega.ml45
-rw-r--r--plugins/micromega/itv.ml80
-rw-r--r--plugins/micromega/itv.mli18
-rw-r--r--plugins/micromega/mfourier.ml224
-rw-r--r--plugins/micromega/mfourier.mli25
-rw-r--r--plugins/micromega/micromega.ml18
-rw-r--r--plugins/micromega/micromega.mli140
-rw-r--r--plugins/micromega/micromega_plugin.mlpack5
-rw-r--r--plugins/micromega/mutils.ml47
-rw-r--r--plugins/micromega/mutils.mli25
-rw-r--r--plugins/micromega/polynomial.ml1336
-rw-r--r--plugins/micromega/polynomial.mli320
-rw-r--r--plugins/micromega/simplex.ml621
-rw-r--r--plugins/micromega/simplex.mli18
-rw-r--r--plugins/micromega/vect.ml295
-rw-r--r--plugins/micromega/vect.mli156
22 files changed, 3349 insertions, 1799 deletions
diff --git a/plugins/micromega/Lia.v b/plugins/micromega/Lia.v
index ae05cf5459..dd6319d5c4 100644
--- a/plugins/micromega/Lia.v
+++ b/plugins/micromega/Lia.v
@@ -32,7 +32,7 @@ Ltac zchange :=
Ltac zchecker_no_abstract := zchange ; vm_compute ; reflexivity.
-Ltac zchecker_abstract := zchange ; vm_cast_no_check (eq_refl true).
+Ltac zchecker_abstract := abstract (zchange ; vm_cast_no_check (eq_refl true)).
Ltac zchecker := zchecker_no_abstract.
diff --git a/plugins/micromega/MExtraction.v b/plugins/micromega/MExtraction.v
index 158ddb589b..5f01f981ef 100644
--- a/plugins/micromega/MExtraction.v
+++ b/plugins/micromega/MExtraction.v
@@ -53,12 +53,11 @@ Extract Constant Rinv => "fun x -> 1 / x".
(** In order to avoid annoying build dependencies the actual
extraction is only performed as a test in the test suite. *)
-(* Extraction "plugins/micromega/micromega.ml" *)
-(* Recursive Extraction *)
-(* List.map simpl_cone (*map_cone indexes*) *)
-(* denorm Qpower vm_add *)
-(* n_of_Z N.of_nat ZTautoChecker ZWeakChecker QTautoChecker RTautoChecker find. *)
-
+(*Extraction "micromega.ml"
+(*Recursive Extraction*) List.map simpl_cone (*map_cone indexes*)
+ denorm Qpower vm_add
+ normZ normQ normQ n_of_Z N.of_nat ZTautoChecker ZWeakChecker QTautoChecker RTautoChecker find.
+*)
(* Local Variables: *)
(* coding: utf-8 *)
(* End: *)
diff --git a/plugins/micromega/QMicromega.v b/plugins/micromega/QMicromega.v
index ddf4064a03..2880a05d8d 100644
--- a/plugins/micromega/QMicromega.v
+++ b/plugins/micromega/QMicromega.v
@@ -179,6 +179,8 @@ Definition qunsat := check_inconsistent 0 Qeq_bool Qle_bool.
Definition qdeduce := nformula_plus_nformula 0 Qplus Qeq_bool.
+Definition normQ := norm 0 1 Qplus Qmult Qminus Qopp Qeq_bool.
+Declare Equivalent Keys normQ RingMicromega.norm.
Definition QTautoChecker (f : BFormula (Formula Q)) (w: list QWitness) : bool :=
diff --git a/plugins/micromega/ZMicromega.v b/plugins/micromega/ZMicromega.v
index 892858e63f..f341a04e03 100644
--- a/plugins/micromega/ZMicromega.v
+++ b/plugins/micromega/ZMicromega.v
@@ -162,8 +162,8 @@ Declare Equivalent Keys psub RingMicromega.psub.
Definition padd := padd Z0 Z.add Zeq_bool.
Declare Equivalent Keys padd RingMicromega.padd.
-Definition norm := norm 0 1 Z.add Z.mul Z.sub Z.opp Zeq_bool.
-Declare Equivalent Keys norm RingMicromega.norm.
+Definition normZ := norm 0 1 Z.add Z.mul Z.sub Z.opp Zeq_bool.
+Declare Equivalent Keys normZ RingMicromega.norm.
Definition eval_pol := eval_pol Z.add Z.mul (fun x => x).
Declare Equivalent Keys eval_pol RingMicromega.eval_pol.
@@ -180,7 +180,7 @@ Proof.
apply (eval_pol_add Zsor ZSORaddon).
Qed.
-Lemma eval_pol_norm : forall env e, eval_expr env e = eval_pol env (norm e) .
+Lemma eval_pol_norm : forall env e, eval_expr env e = eval_pol env (normZ e) .
Proof.
intros.
apply (eval_pol_norm Zsor ZSORaddon).
@@ -188,8 +188,8 @@ Qed.
Definition xnormalise (t:Formula Z) : list (NFormula Z) :=
let (lhs,o,rhs) := t in
- let lhs := norm lhs in
- let rhs := norm rhs in
+ let lhs := normZ lhs in
+ let rhs := normZ rhs in
match o with
| OpEq =>
((psub lhs (padd rhs (Pc 1))),NonStrict)::((psub rhs (padd lhs (Pc 1))),NonStrict)::nil
@@ -225,8 +225,8 @@ Qed.
Definition xnegate (t:RingMicromega.Formula Z) : list (NFormula Z) :=
let (lhs,o,rhs) := t in
- let lhs := norm lhs in
- let rhs := norm rhs in
+ let lhs := normZ lhs in
+ let rhs := normZ rhs in
match o with
| OpEq => (psub lhs rhs,Equal) :: nil
| OpNEq => ((psub lhs (padd rhs (Pc 1))),NonStrict)::((psub rhs (padd lhs (Pc 1))),NonStrict)::nil
diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml
index 3a9709b6ce..e6edd50878 100644
--- a/plugins/micromega/certificate.ml
+++ b/plugins/micromega/certificate.ml
@@ -28,109 +28,80 @@ module Mc = Micromega
module Ml2C = Mutils.CamlToCoq
module C2Ml = Mutils.CoqToCaml
+let use_simplex = ref true
+
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 -> bool
-}
+ bigint_to_number : big_int -> 'a;
+ number_to_num : 'a -> num;
+ zero : 'a;
+ unit : 'a;
+ mult : 'a -> 'a -> 'a;
+ eqb : 'a -> 'a -> 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.Z.mul;
- eqb = Mc.zeq_bool
-}
-
+ 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.Z.mul;
+ 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
-}
+ 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 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 (C2Ml.positive 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 Int.equal 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 rec dev_form p =
+ match p with
+ | Mc.PEc z -> Poly.constant (n_spec.number_to_num z)
+ | Mc.PEX v -> Poly.variable (C2Ml.positive 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 Int.equal 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 rec fixpoint f x =
- let y' = f x in
- if Pervasives.(=) y' x then y'
- else fixpoint f y'
+ let y' = f x in
+ if Pervasives.(=) 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.PsatzMulE(t1, t2) ->
- simpl_cone (Mc.PsatzMulE (rec_simpl_cone t1, rec_simpl_cone t2))
- | Mc.PsatzAdd(t1,t2) ->
- simpl_cone (Mc.PsatzAdd (rec_simpl_cone t1, rec_simpl_cone t2))
- | x -> simpl_cone x in
- rec_simpl_cone 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.PsatzMulE(t1, t2) ->
+ simpl_cone (Mc.PsatzMulE (rec_simpl_cone t1, rec_simpl_cone t2))
+ | Mc.PsatzAdd(t1,t2) ->
+ simpl_cone (Mc.PsatzAdd (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
-let factorise_linear_cone c =
-
- let rec cone_list c l =
- match c with
- | Mc.PsatzAdd (x,r) -> cone_list r (x::l)
- | _ -> c :: l in
-
- let factorise c1 c2 =
- match c1 , c2 with
- | Mc.PsatzMulC(x,y) , Mc.PsatzMulC(x',y') ->
- if Pervasives.(=) x x' then Some (Mc.PsatzMulC(x, Mc.PsatzAdd(y,y'))) else None
- | Mc.PsatzMulE(x,y) , Mc.PsatzMulE(x',y') ->
- if Pervasives.(=) x x' then Some (Mc.PsatzMulE(x, Mc.PsatzAdd(y,y'))) else None
- | _ -> None in
-
- let rec rebuild_cone l pending =
- match l with
- | [] -> (match pending with
- | None -> Mc.PsatzZ
- | Some p -> p
- )
- | e::l ->
- (match pending with
- | None -> rebuild_cone l (Some e)
- | Some p -> (match factorise p e with
- | None -> Mc.PsatzAdd(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
@@ -147,956 +118,921 @@ let factorise_linear_cone c =
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
+ The variable c is at index 1
+ *)
(* 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 Pervasives.(=) mn Monomial.const
- then
- { coeffs = Vect.from_list ((Big_int unit_big_int):: (List.rev coeffs)) ;
+let constrain_variable v l =
+ let coeffs = List.fold_left (fun acc p -> (Vect.get v p.coeffs)::acc) [] l in
+ { coeffs = Vect.from_list ((Big_int zero_big_int):: (Big_int zero_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)) ;
+
+
+
+let constrain_constant l =
+ let coeffs = List.fold_left (fun acc p -> minus_num p.cst ::acc) [] l in
+ { coeffs = Vect.from_list ((Big_int zero_big_int):: (Big_int unit_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 rec xpositivity i l =
+ match l with
+ | [] -> []
+ | c::l -> match c.op with
+ | Eq -> xpositivity (i+1) l
+ | _ ->
+ {coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ;
+ op = Ge ;
+ cst = Int 0 } :: (xpositivity (i+1) l)
+ in
+ xpositivity 1 l
+
+
+let cstr_of_poly (p,o) =
+ let (c,l) = Vect.decomp_cst p in
+ {coeffs = l; op = o ; cst = minus_num c}
+
+
+
+let variables_of_cstr c = Vect.variables c.coeffs
-module MonSet = Set.Make(Monomial)
(* 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 ==> This does not work anymore *)
- let l' = List.map fst l in
-
- let monomials =
- List.fold_left (fun acc p ->
- Poly.fold (fun m _ acc -> MonSet.add m acc) p acc)
- (MonSet.singleton Monomial.const) l'
- in (* For each monomial, compute a constraint *)
- let s0 =
- MonSet.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
+
+let build_dual_linear_system l =
+
+ let variables =
+ List.fold_left (fun acc p -> ISet.union acc (variables_of_cstr p)) ISet.empty l in
+ (* For each monomial, compute a constraint *)
+ let s0 =
+ ISet.fold (fun mn res -> (constrain_variable mn l)::res) variables [] in
+ let c = constrain_constant l in
+
+ (* I need at least something strictly positive *)
+ let strict = {
+ coeffs = Vect.from_list ((Big_int zero_big_int) :: (Big_int unit_big_int)::
+ (List.map (fun c -> if is_strict c then Big_int unit_big_int else 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)
-
-(* 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
- | [] -> failwith "empty_certificate"
- | e::cert' ->
- (* let cst = match compare_big_int e zero_big_int with
- | 0 -> Mc.PsatzZ
- | 1 -> Mc.PsatzC (bint_to_cst e)
- | _ -> failwith "positivity error"
- in *)
- let rec scalar_product cert l =
- match cert with
- | [] -> Mc.PsatzZ
- | 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.PsatzAdd (
- Mc.PsatzMulC (Mc.Pc ( bint_to_cst c), Mc.PsatzIn (Ml2C.nat i)),
- r)
- | 0 -> r
- | _ -> Mc.PsatzAdd (
- Mc.PsatzMulE (Mc.PsatzC (bint_to_cst c), Mc.PsatzIn (Ml2C.nat i)),
- r) in
- (factorise_linear_cone
- (simplify_cone n_spec (scalar_product cert' li)))
-
-
-exception Strict
-
-module MonMap = Map.Make(Monomial)
-
-let primal l =
- let vr = ref 0 in
-
- let vect_of_poly map p =
- Poly.fold (fun mn vl (map,vect) ->
- if Pervasives.(=) mn Monomial.const
- then (map,vect)
- else
- let (mn,m) = try (MonMap.find mn map,map) with Not_found -> let res = (!vr, MonMap.add mn !vr map) in incr vr ; res in
- (m,if Int.equal (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 = Int.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 (MonMap.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
- match Fourier.find_point sys with
- | Inr _ -> None
- | Inl cert -> Some (rats_to_ints (Vect.to_list cert))
- (* should not use rats_to_ints *)
- with x when CErrors.noncritical x ->
- if debug
- then (Printf.printf "raw certificate %s" (Printexc.to_string x);
- flush stdout) ;
- 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))
+ {coeffs = Vect.from_list ([Big_int zero_big_int ;Big_int unit_big_int]) ;
+ op = Ge ;
+ cst = Big_int zero_big_int}::(strict::(positivity l)@c::s0)
+
+
+(** [direct_linear_prover l] does not handle strict inegalities *)
+let fourier_linear_prover l =
+ match Mfourier.Fourier.find_point l with
+ | Inr prf ->
+ if debug then Printf.printf "AProof : %a\n" Mfourier.pp_proof prf ;
+ let cert = (*List.map (fun (x,n) -> x+1,n)*) (fst (List.hd (Mfourier.Proof.mk_proof l prf))) in
+ if debug then Printf.printf "CProof : %a" Vect.pp cert ;
+ (*Some (rats_to_ints (Vect.to_list cert))*)
+ Some (Vect.normalise cert)
| Inl _ -> None
- with Strict ->
+
+
+let direct_linear_prover l =
+ if !use_simplex
+ then Simplex.find_unsat_certificate l
+ else fourier_linear_prover l
+
+let find_point l =
+ if !use_simplex
+ then Simplex.find_point l
+ else match Mfourier.Fourier.find_point l with
+ | Inr _ -> None
+ | Inl cert -> Some cert
+
+let optimise v l =
+ if !use_simplex
+ then Simplex.optimise v l
+ else Mfourier.Fourier.optimise v l
+
+
+
+let dual_raw_certificate l =
+ if debug
+ then begin
+ Printf.printf "dual_raw_certificate\n";
+ List.iter (fun c -> Printf.fprintf stdout "%a\n" output_cstr c) l
+ end;
+
+ let sys = build_dual_linear_system l in
+
+ if debug then begin
+ Printf.printf "dual_system\n";
+ List.iter (fun c -> Printf.fprintf stdout "%a\n" output_cstr c) sys
+ end;
+
+ try
+ match find_point sys with
+ | None -> None
+ | Some cert ->
+ match Vect.choose cert with
+ | None -> failwith "dual_raw_certificate: empty_certificate"
+ | Some _ ->
+ (*Some (rats_to_ints (Vect.to_list (Vect.decr_var 2 (Vect.set 1 (Int 0) cert))))*)
+ Some (Vect.normalise (Vect.decr_var 2 (Vect.set 1 (Int 0) cert)))
+ (* should not use rats_to_ints *)
+ with x when CErrors.noncritical x ->
+ if debug
+ then (Printf.printf "dual raw certificate %s" (Printexc.to_string x);
+ flush stdout) ;
+ None
+
+
+
+let simple_linear_prover l =
+ try
+ direct_linear_prover l
+ with Strict ->
(* Fourier elimination should handle > *)
- dual_raw_certificate l
+ dual_raw_certificate l
+open ProofFormat
+
+
+let env_of_list l =
+ snd (List.fold_left (fun (i,m) p -> (i+1,IMap.add i p m)) (0,IMap.empty) l)
-let simple_linear_prover l =
- let (lc,li) = List.split l in
- match raw_certificate lc with
- | None -> None (* No certificate *)
- | Some cert -> Some (cert,li)
-
+let linear_prover_cstr sys =
+ let (sysi,prfi) = List.split sys in
-let linear_prover n_spec l =
- let build_system n_spec l =
- let li = List.combine l (CList.interval 0 (List.length l -1)) in
- let (l1,l') = List.partition
- (fun (x,_) -> if Pervasives.(=) (snd x) Mc.NonEqual then true else false) li in
- List.map
- (fun ((x,y),i) -> match y with
- Mc.NonEqual -> failwith "cannot happen"
- | y -> ((dev_form n_spec x, y),i)) l' in
- let l' = build_system n_spec l in
- simple_linear_prover (*n_spec*) l'
+
+ match simple_linear_prover sysi with
+ | None -> None
+ | Some cert -> Some (proof_of_farkas (env_of_list prfi) cert)
+
+let linear_prover_cstr =
+ if debug
+ then
+ fun sys ->
+ Printf.printf "<linear_prover"; flush stdout ;
+ let res = linear_prover_cstr sys in
+ Printf.printf ">"; flush stdout ;
+ res
+ else linear_prover_cstr
-let linear_prover n_spec l =
- try linear_prover n_spec l
- with x when CErrors.noncritical x ->
- (print_string (Printexc.to_string x); None)
let compute_max_nb_cstr l d =
- let len = List.length l in
- max len (max d (len * d))
+ let len = List.length l in
+ max len (max d (len * d))
-let linear_prover_with_cert prfdepth spec l =
- max_nb_cstr := compute_max_nb_cstr l prfdepth ;
- match linear_prover spec l with
- | None -> None
- | Some cert -> Some (make_certificate spec cert)
-let nlinear_prover prfdepth (sys: (Mc.q Mc.pExpr * Mc.op1) list) =
- LinPoly.MonT.clear ();
- max_nb_cstr := compute_max_nb_cstr sys prfdepth ;
- (* Assign a proof to the initial hypotheses *)
- let sys = List.mapi (fun i c -> (c,Mc.PsatzIn (Ml2C.nat i))) sys in
+let develop_constraint z_spec (e,k) =
+ (dev_form z_spec e,
+ match k with
+ | Mc.NonStrict -> Ge
+ | Mc.Equal -> Eq
+ | Mc.Strict -> Gt
+ | _ -> assert false
+ )
+
+(** A single constraint can be unsat for the following reasons:
+ - 0 >= c for c a negative constant
+ - 0 = c for c a non-zero constant
+ - e = c when the coeffs of e are all integers and c is rational
+ *)
+open ProofFormat
+type checksat =
+ | Tauto (* Tautology *)
+ | Unsat of prf_rule (* Unsatisfiable *)
+ | Cut of cstr * prf_rule (* Cutting plane *)
+ | Normalise of cstr * prf_rule (* Coefficients may be normalised i.e relatively prime *)
- (* Add all the product of hypotheses *)
- let prod = all_pairs (fun ((c,o),p) ((c',o'),p') ->
- ((Mc.PEmul(c,c') , Mc.opMult o o') , Mc.PsatzMulE(p,p'))) sys in
+exception FoundProof of prf_rule
+
+
+(** [check_sat]
+ - detects constraints that are not satisfiable;
+ - normalises constraints and generate cuts.
+ *)
+
+let check_int_sat (cstr,prf) =
+ let {coeffs=coeffs ; op=op ; cst=cst} = cstr in
+ match Vect.choose coeffs with
+ | None ->
+ if eval_op op (Int 0) cst then Tauto else Unsat prf
+ | _ ->
+ let gcdi = Vect.gcd coeffs in
+ let gcd = Big_int gcdi in
+ if eq_num gcd (Int 1)
+ then Normalise(cstr,prf)
+ else
+ if Int.equal (sign_num (mod_num cst gcd)) 0
+ then (* We can really normalise *)
+ begin
+ assert (sign_num gcd >=1 ) ;
+ let cstr = {
+ coeffs = Vect.div gcd coeffs;
+ op = op ; cst = cst // gcd
+ } in
+ Normalise(cstr,Gcd(gcdi,prf))
+ (* Normalise(cstr,CutPrf prf)*)
+ end
+ else
+ match op with
+ | Eq -> Unsat (CutPrf prf)
+ | Ge ->
+ let cstr = {
+ coeffs = Vect.div gcd coeffs;
+ op = op ; cst = ceiling_num (cst // gcd)
+ } in Cut(cstr,CutPrf prf)
+ | Gt -> failwith "check_sat : Unexpected operator"
+
+
+let apply_and_normalise check f psys =
+ List.fold_left (fun acc pc' ->
+ match f pc' with
+ | None -> pc'::acc
+ | Some pc' ->
+ match check pc' with
+ | Tauto -> acc
+ | Unsat prf -> raise (FoundProof prf)
+ | Cut(c,p) -> (c,p)::acc
+ | Normalise (c,p) -> (c,p)::acc
+ ) [] psys
+
+
+let simplify f sys =
+ let (sys',b) =
+ List.fold_left (fun (sys',b) c ->
+ match f c with
+ | None -> (c::sys',b)
+ | Some c' ->
+ (c'::sys',true)
+ ) ([],false) sys in
+ if b then Some sys' else None
+
+let saturate f sys =
+ List.fold_left (fun sys' c -> match f c with
+ | None -> sys'
+ | Some c' -> c'::sys'
+ ) [] sys
+
+let is_substitution strict ((p,o),prf) =
+ let pred v = if strict then v =/ Int 1 || v =/ Int (-1) else true in
- (* Only filter those have a meaning *)
- let prod = List.fold_left (fun l ((c,o),p) ->
match o with
- | None -> l
- | Some o -> ((c,o),p) :: l) [] prod in
-
- let sys = sys @ prod in
-
- let square =
- (* Collect the squares and state that they are positive *)
- let pols = List.map (fun ((p,_),_) -> dev_form q_spec p) sys in
- let square =
- List.fold_left (fun acc p ->
- Poly.fold
- (fun m _ acc ->
- match Monomial.sqrt m with
- | None -> acc
- | Some s -> MonMap.add s m acc) p acc) MonMap.empty pols in
-
- let pol_of_mon m =
- Monomial.fold (fun x v p -> Mc.PEmul(Mc.PEpow(Mc.PEX(Ml2C.positive x),Ml2C.n v),p)) m (Mc.PEc q_spec.unit) in
-
- let norm0 =
- Mc.norm q_spec.zero q_spec.unit Mc.qplus Mc.qmult Mc.qminus Mc.qopp Mc.qeq_bool in
+ | Eq -> LinPoly.search_linear pred p
+ | _ -> None
+
+
+let is_linear_for v pc =
+ LinPoly.is_linear (fst (fst pc)) || LinPoly.is_linear_for v (fst (fst pc))
+
+
+
+
+let non_linear_pivot sys pc v pc' =
+ if LinPoly.is_linear (fst (fst pc'))
+ then None (* There are other ways to deal with those *)
+ else WithProof.linear_pivot sys pc v pc'
+
+
+let is_linear_substitution sys ((p,o),prf) =
+ let pred v = v =/ Int 1 || v =/ Int (-1) in
+ match o with
+ | Eq -> begin
+ match
+ List.filter (fun v -> List.for_all (is_linear_for v) sys) (LinPoly.search_all_linear pred p)
+ with
+ | [] -> None
+ | v::_ -> Some v (* make a choice *)
+ end
+ | _ -> None
+
+
+let elim_simple_linear_equality sys0 =
+
+ let elim sys =
+ let (oeq,sys') = extract (is_linear_substitution sys) sys in
+ match oeq with
+ | None -> None
+ | Some(v,pc) -> simplify (WithProof.linear_pivot sys0 pc v) sys' in
+
+ iterate_until_stable elim sys0
+
+
+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 develop_constraints prfdepth n_spec sys =
+ LinPoly.MonT.clear ();
+ max_nb_cstr := compute_max_nb_cstr sys prfdepth ;
+ let sys = List.map (develop_constraint n_spec) sys in
+ List.mapi (fun i (p,o) -> ((LinPoly.linpol_of_pol p,o),Hyp i)) sys
+
+let square_of_var i =
+ let x = LinPoly.var i in
+ ((LinPoly.product x x,Ge),(Square x))
+
+(** [nlinear_preprocess sys] augments the system [sys] by performing some limited non-linear reasoning.
+ For instance, it asserts that the x² ≥0 but also that if c₁ ≥ 0 ∈ sys and c₂ ≥ 0 ∈ sys then c₁ × c₂ ≥ 0.
+ The resulting system is linearised.
+ *)
+
+let nlinear_preprocess (sys:WithProof.t list) =
+
+ let is_linear = List.for_all (fun ((p,_),_) -> LinPoly.is_linear p) sys in
+
+ if is_linear then sys
+ else
+ let collect_square =
+ List.fold_left (fun acc ((p,_),_) -> MonMap.union (fun k e1 e2 -> Some e1) acc (LinPoly.collect_square p)) MonMap.empty sys in
+ let sys = MonMap.fold (fun s m acc ->
+ let s = LinPoly.of_monomial s in
+ let m = LinPoly.of_monomial m in
+ ((m, Ge), (Square s))::acc) collect_square sys in
+
+ let collect_vars = List.fold_left (fun acc p -> ISet.union acc (LinPoly.variables (fst (fst p)))) ISet.empty sys in
+
+ let sys = ISet.fold (fun i acc -> square_of_var i :: acc) collect_vars sys in
+
+ let sys = sys @ (all_pairs WithProof.product sys) in
- MonMap.fold (fun s m acc -> ((pol_of_mon m , Mc.NonStrict), Mc.PsatzSquare(norm0 (pol_of_mon s)))::acc) square [] in
+ if debug then begin
+ Printf.fprintf stdout "Preprocessed\n";
+ List.iter (fun s -> Printf.fprintf stdout "%a\n" WithProof.output s) sys;
+ end ;
+
+ List.map (WithProof.annot "P") sys
+
- let sys = sys @ square in
+let nlinear_prover prfdepth sys =
+ let sys = develop_constraints prfdepth q_spec sys in
+ let sys1 = elim_simple_linear_equality sys in
+ let sys2 = saturate_linear_equality_non_linear sys1 in
+ let sys = nlinear_preprocess sys1@sys2 in
+ let sys = List.map (fun ((p,o),prf) -> (cstr_of_poly (p,o), prf)) sys in
+ let id = (List.fold_left
+ (fun acc (_,r) -> max acc (ProofFormat.pr_rule_max_id r)) 0 sys) in
+ let env = CList.interval 0 id in
+ match linear_prover_cstr sys with
+ | None -> None
+ | Some cert ->
+ Some (cmpl_prf_rule Mc.normQ CamlToCoq.q env cert)
- (* Call the linear prover without the proofs *)
- let sys_no_prf = List.map fst sys in
- match linear_prover q_spec sys_no_prf with
- | None -> None
- | Some cert ->
- let cert = make_certificate q_spec cert in
- let rec map_psatz = function
- | Mc.PsatzIn n -> snd (List.nth sys (C2Ml.nat n))
- | Mc.PsatzSquare c -> Mc.PsatzSquare c
- | Mc.PsatzMulC(c,p) -> Mc.PsatzMulC(c, map_psatz p)
- | Mc.PsatzMulE(p1,p2) -> Mc.PsatzMulE(map_psatz p1,map_psatz p2)
- | Mc.PsatzAdd(p1,p2) -> Mc.PsatzAdd(map_psatz p1,map_psatz p2)
- | Mc.PsatzC c -> Mc.PsatzC c
- | Mc.PsatzZ -> Mc.PsatzZ in
- Some (map_psatz cert)
+let linear_prover_with_cert prfdepth sys =
+ let sys = develop_constraints prfdepth q_spec sys in
+ (* let sys = nlinear_preprocess sys in *)
+ let sys = List.map (fun (c,p) -> cstr_of_poly c,p) sys in
+
+ match linear_prover_cstr sys with
+ | None -> None
+ | Some cert ->
+ Some (cmpl_prf_rule Mc.normQ CamlToCoq.q (List.mapi (fun i e -> i) sys) cert)
(* The prover is (probably) incomplete --
only searching for naive cutting planes *)
-let develop_constraint z_spec (e,k) =
- match k with
- | Mc.NonStrict -> (dev_form z_spec e , Ge)
- | Mc.Equal -> (dev_form z_spec e , Eq)
- | _ -> assert false
-
open Sos_types
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 Int.equal (compare_big_int e unit_big_int) 0
- then (unit_big_int, Add (y1,y2))
- else e, Add (Mul(Const (Big_int s2'), y1),
+ 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 Int.equal (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"
+ | 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 (s,t') = scale_term t in
+ s,t'
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)
+ | 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"
+ | 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 term_to_q_pol e = Mc.norm_aux (Ml2C.q (Int 0)) (Ml2C.q (Int 1)) Mc.qplus Mc.qmult Mc.qminus Mc.qopp Mc.qeq_bool (term_to_q_expr e)
let rec product l =
- match l with
- | [] -> Mc.PsatzZ
- | [i] -> Mc.PsatzIn (Ml2C.nat i)
- | i ::l -> Mc.PsatzMulE(Mc.PsatzIn (Ml2C.nat i), product l)
+ match l with
+ | [] -> Mc.PsatzZ
+ | [i] -> Mc.PsatzIn (Ml2C.nat i)
+ | i ::l -> Mc.PsatzMulE(Mc.PsatzIn (Ml2C.nat i), product l)
let q_cert_of_pos pos =
- let rec _cert_of_pos = function
- Axiom_eq i -> Mc.PsatzIn (Ml2C.nat i)
- | Axiom_le i -> Mc.PsatzIn (Ml2C.nat i)
- | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i)
- | Monoid l -> product l
- | Rational_eq n | Rational_le n | Rational_lt n ->
- if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ else
- Mc.PsatzC (Ml2C.q n)
- | Square t -> Mc.PsatzSquare (term_to_q_pol t)
- | Eqmul (t, y) -> Mc.PsatzMulC(term_to_q_pol t, _cert_of_pos y)
- | Sum (y, z) -> Mc.PsatzAdd (_cert_of_pos y, _cert_of_pos z)
- | Product (y, z) -> Mc.PsatzMulE (_cert_of_pos y, _cert_of_pos z) in
- simplify_cone q_spec (_cert_of_pos pos)
+ let rec _cert_of_pos = function
+ Axiom_eq i -> Mc.PsatzIn (Ml2C.nat i)
+ | Axiom_le i -> Mc.PsatzIn (Ml2C.nat i)
+ | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i)
+ | Monoid l -> product l
+ | Rational_eq n | Rational_le n | Rational_lt n ->
+ if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ else
+ Mc.PsatzC (Ml2C.q n)
+ | Square t -> Mc.PsatzSquare (term_to_q_pol t)
+ | Eqmul (t, y) -> Mc.PsatzMulC(term_to_q_pol t, _cert_of_pos y)
+ | Sum (y, z) -> Mc.PsatzAdd (_cert_of_pos y, _cert_of_pos z)
+ | Product (y, z) -> Mc.PsatzMulE (_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"
+ | 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 term_to_z_pol e = Mc.norm_aux (Ml2C.z 0) (Ml2C.z 1) Mc.Z.add Mc.Z.mul Mc.Z.sub Mc.Z.opp Mc.zeq_bool (term_to_z_expr e)
let z_cert_of_pos pos =
- let s,pos = (scale_certificate pos) in
- let rec _cert_of_pos = function
- Axiom_eq i -> Mc.PsatzIn (Ml2C.nat i)
- | Axiom_le i -> Mc.PsatzIn (Ml2C.nat i)
- | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i)
- | Monoid l -> product l
- | Rational_eq n | Rational_le n | Rational_lt n ->
- if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ else
- Mc.PsatzC (Ml2C.bigint (big_int_of_num n))
- | Square t -> Mc.PsatzSquare (term_to_z_pol t)
- | Eqmul (t, y) ->
- let is_unit =
- match t with
- | Const n -> n =/ Int 1
- | _ -> false in
- if is_unit
- then _cert_of_pos y
- else Mc.PsatzMulC(term_to_z_pol t, _cert_of_pos y)
- | Sum (y, z) -> Mc.PsatzAdd (_cert_of_pos y, _cert_of_pos z)
- | Product (y, z) -> Mc.PsatzMulE (_cert_of_pos y, _cert_of_pos z) in
- simplify_cone z_spec (_cert_of_pos pos)
+ let s,pos = (scale_certificate pos) in
+ let rec _cert_of_pos = function
+ Axiom_eq i -> Mc.PsatzIn (Ml2C.nat i)
+ | Axiom_le i -> Mc.PsatzIn (Ml2C.nat i)
+ | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i)
+ | Monoid l -> product l
+ | Rational_eq n | Rational_le n | Rational_lt n ->
+ if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ else
+ Mc.PsatzC (Ml2C.bigint (big_int_of_num n))
+ | Square t -> Mc.PsatzSquare (term_to_z_pol t)
+ | Eqmul (t, y) ->
+ let is_unit =
+ match t with
+ | Const n -> n =/ Int 1
+ | _ -> false in
+ if is_unit
+ then _cert_of_pos y
+ else Mc.PsatzMulC(term_to_z_pol t, _cert_of_pos y)
+ | Sum (y, z) -> Mc.PsatzAdd (_cert_of_pos y, _cert_of_pos z)
+ | Product (y, z) -> Mc.PsatzMulE (_cert_of_pos y, _cert_of_pos z) in
+ simplify_cone z_spec (_cert_of_pos pos)
(** All constraints (initial or derived) have an index and have a justification i.e., proof.
Given a constraint, all the coefficients are always integers.
-*)
+ *)
open Mutils
-open Mfourier
open Num
open Big_int
open Polynomial
-module Env =
-struct
-
- let id_of_hyp hyp l =
- let rec xid_of_hyp i l =
- match l with
- | [] -> failwith "id_of_hyp"
- | hyp'::l -> if Pervasives.(=) hyp hyp' then i else xid_of_hyp (i+1) l in
- xid_of_hyp 0 l
-
-end
+type prf_sys = (cstr * prf_rule) list
-let coq_poly_of_linpol (p,c) =
-
- let pol_of_mon m =
- Monomial.fold (fun x v p -> Mc.PEmul(Mc.PEpow(Mc.PEX(Ml2C.positive x),Ml2C.n v),p)) m (Mc.PEc (Mc.Zpos Mc.XH)) in
-
- List.fold_left (fun acc (x,v) ->
- let mn = LinPoly.MonT.retrieve x in
- Mc.PEadd(Mc.PEmul(Mc.PEc (Ml2C.bigint (numerator v)), pol_of_mon mn),acc)) (Mc.PEc (Ml2C.bigint (numerator c))) p
-
-
-
-
-let rec cmpl_prf_rule env = function
- | Hyp i | Def i -> Mc.PsatzIn (Ml2C.nat (Env.id_of_hyp i env))
- | Cst i -> Mc.PsatzC (Ml2C.bigint i)
- | Zero -> Mc.PsatzZ
- | MulPrf(p1,p2) -> Mc.PsatzMulE(cmpl_prf_rule env p1, cmpl_prf_rule env p2)
- | AddPrf(p1,p2) -> Mc.PsatzAdd(cmpl_prf_rule env p1 , cmpl_prf_rule env p2)
- | MulC(lp,p) -> let lp = Mc.norm0 (coq_poly_of_linpol lp) in
- Mc.PsatzMulC(lp,cmpl_prf_rule env p)
- | Square lp -> Mc.PsatzSquare (Mc.norm0 (coq_poly_of_linpol lp))
- | _ -> failwith "Cuts should already be compiled"
-
-
-let rec cmpl_proof env = function
- | Done -> Mc.DoneProof
- | Step(i,p,prf) ->
- begin
- match p with
- | CutPrf p' ->
- Mc.CutProof(cmpl_prf_rule env p', cmpl_proof (i::env) prf)
- | _ -> Mc.RatProof(cmpl_prf_rule env p,cmpl_proof (i::env) prf)
- end
- | Enum(i,p1,_,p2,l) ->
- Mc.EnumProof(cmpl_prf_rule env p1,cmpl_prf_rule env p2,List.map (cmpl_proof (i::env)) l)
-
-
-let compile_proof env prf =
- let id = 1 + proof_max_id prf in
- let _,prf = normalise_proof id prf in
- if debug then Printf.fprintf stdout "compiled proof %a\n" output_proof prf;
- cmpl_proof env prf
-
-type prf_sys = (cstr_compat * prf_rule) list
-
-
-let xlinear_prover sys =
- match Fourier.find_point sys 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 sys prf))) in
- if debug then Printf.printf "CProof : %a" Vect.pp_vect cert ;
- Some (rats_to_ints (Vect.to_list cert))
- | Inl _ -> None
-
-
-let proof_of_farkas prf cert =
- (* Printf.printf "\nproof_of_farkas %a , %a \n" (pp_list output_prf_rule) prf (pp_list output_bigint) cert ; *)
- let rec mk_farkas acc prf cert =
- match prf, cert with
- | _ , [] -> acc
- | [] , _ -> failwith "proof_of_farkas : not enough hyps"
- | p::prf,c::cert ->
- mk_farkas (add_proof (mul_proof c p) acc) prf cert in
- let res = mk_farkas Zero prf cert in
- (*Printf.printf "==> %a" output_prf_rule res ; *)
- res
-
-
-let linear_prover sys =
- let (sysi,prfi) = List.split sys in
- match xlinear_prover sysi with
- | None -> None
- | Some cert -> Some (proof_of_farkas prfi cert)
-
-let linear_prover =
- if debug
- then
- fun sys ->
- Printf.printf "<linear_prover"; flush stdout ;
- let res = linear_prover sys in
- Printf.printf ">"; flush stdout ;
- res
- else linear_prover
-
-
-
-
-(** A single constraint can be unsat for the following reasons:
- - 0 >= c for c a negative constant
- - 0 = c for c a non-zero constant
- - e = c when the coeffs of e are all integers and c is rational
-*)
-
-type checksat =
-| Tauto (* Tautology *)
-| Unsat of prf_rule (* Unsatisfiable *)
-| Cut of cstr_compat * prf_rule (* Cutting plane *)
-| Normalise of cstr_compat * prf_rule (* coefficients are relatively prime *)
-
-
-(** [check_sat]
- - detects constraints that are not satisfiable;
- - normalises constraints and generate cuts.
-*)
-
-let check_sat (cstr,prf) =
- let {coeffs=coeffs ; op=op ; cst=cst} = cstr in
- match coeffs with
- | [] ->
- if eval_op op (Int 0) cst then Tauto else Unsat prf
- | _ ->
- let gcdi = (gcd_list (List.map snd coeffs)) in
- let gcd = Big_int gcdi in
- if eq_num gcd (Int 1)
- then Normalise(cstr,prf)
- else
- if Int.equal (sign_num (mod_num cst gcd)) 0
- then (* We can really normalise *)
- begin
- assert (sign_num gcd >=1 ) ;
- let cstr = {
- coeffs = List.map (fun (x,v) -> (x, v // gcd)) coeffs;
- op = op ; cst = cst // gcd
- } in
- Normalise(cstr,Gcd(gcdi,prf))
- (* Normalise(cstr,CutPrf prf)*)
- end
- else
- match op with
- | Eq -> Unsat (CutPrf prf)
- | Ge ->
- let cstr = {
- coeffs = List.map (fun (x,v) -> (x, v // gcd)) coeffs;
- op = op ; cst = ceiling_num (cst // gcd)
- } in Cut(cstr,CutPrf prf)
(** Proof generating pivoting over variable v *)
let pivot v (c1,p1) (c2,p2) =
- let {coeffs = v1 ; op = op1 ; cst = n1} = c1
- and {coeffs = v2 ; op = op2 ; cst = n2} = c2 in
+ let {coeffs = v1 ; op = op1 ; cst = n1} = c1
+ and {coeffs = v2 ; op = op2 ; cst = n2} = c2 in
(* Could factorise gcd... *)
- let xpivot cv1 cv2 =
- (
- {coeffs = Vect.add (Vect.mul cv1 v1) (Vect.mul cv2 v2) ;
- op = Proof.add_op op1 op2 ;
- cst = n1 */ cv1 +/ n2 */ cv2 },
-
- AddPrf(mul_proof (numerator cv1) p1,mul_proof (numerator cv2) p2)) in
-
- match Vect.get v v1 , Vect.get v v2 with
- | None , _ | _ , None -> None
- | Some a , Some b ->
- if Int.equal ((sign_num a) * (sign_num b)) (-1)
- then
- let cv1 = abs_num b
- and cv2 = abs_num a in
- Some (xpivot cv1 cv2)
- else
- if op1 == Eq
- then
- let cv1 = minus_num (b */ (Int (sign_num a)))
- and cv2 = abs_num a in
- Some (xpivot cv1 cv2)
- else if op2 == Eq
- then
- let cv1 = abs_num b
- and cv2 = minus_num (a */ (Int (sign_num b))) in
- Some (xpivot cv1 cv2)
- else None (* op2 could be Eq ... this might happen *)
-
-exception FoundProof of prf_rule
+ let xpivot cv1 cv2 =
+ (
+ {coeffs = Vect.add (Vect.mul cv1 v1) (Vect.mul cv2 v2) ;
+ op = opAdd op1 op2 ;
+ cst = n1 */ cv1 +/ n2 */ cv2 },
+
+ AddPrf(mul_cst_proof cv1 p1,mul_cst_proof cv2 p2)) in
+
+ match Vect.get v v1 , Vect.get v v2 with
+ | Int 0 , _ | _ , Int 0 -> None
+ | a , b ->
+ if Int.equal ((sign_num a) * (sign_num b)) (-1)
+ then
+ let cv1 = abs_num b
+ and cv2 = abs_num a in
+ Some (xpivot cv1 cv2)
+ else
+ if op1 == Eq
+ then
+ let cv1 = minus_num (b */ (Int (sign_num a)))
+ and cv2 = abs_num a in
+ Some (xpivot cv1 cv2)
+ else if op2 == Eq
+ then
+ let cv1 = abs_num b
+ and cv2 = minus_num (a */ (Int (sign_num b))) in
+ Some (xpivot cv1 cv2)
+ else None (* op2 could be Eq ... this might happen *)
+
let simpl_sys sys =
- List.fold_left (fun acc (c,p) ->
- match check_sat (c,p) with
- | Tauto -> acc
- | Unsat prf -> raise (FoundProof prf)
- | Cut(c,p) -> (c,p)::acc
- | Normalise (c,p) -> (c,p)::acc) [] sys
+ List.fold_left (fun acc (c,p) ->
+ match check_int_sat (c,p) with
+ | Tauto -> acc
+ | Unsat prf -> raise (FoundProof prf)
+ | Cut(c,p) -> (c,p)::acc
+ | Normalise (c,p) -> (c,p)::acc) [] sys
(** [ext_gcd a b] is the extended Euclid algorithm.
[ext_gcd a b = (x,y,g)] iff [ax+by=g]
Source: http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
-*)
+ *)
let rec ext_gcd a b =
- if Int.equal (sign_big_int b) 0
- then (unit_big_int,zero_big_int)
- else
- let (q,r) = quomod_big_int a b in
- let (s,t) = ext_gcd b r in
- (t, sub_big_int s (mult_big_int q t))
+ if Int.equal (sign_big_int b) 0
+ then (unit_big_int,zero_big_int)
+ else
+ let (q,r) = quomod_big_int a b in
+ let (s,t) = ext_gcd b r in
+ (t, sub_big_int s (mult_big_int q t))
let extract_coprime (c1,p1) (c2,p2) =
- let rec exist2 vect1 vect2 =
- match vect1 , vect2 with
- | _ , [] | [], _ -> None
- | (v1,n1)::vect1' , (v2, n2) :: vect2' ->
- if Pervasives.(=) v1 v2
- then
- if Int.equal (compare_big_int (gcd_big_int (numerator n1) (numerator n2)) unit_big_int) 0
- then Some (v1,n1,n2)
- else
- exist2 vect1' vect2'
- else
- if v1 < v2
- then exist2 vect1' vect2
- else exist2 vect1 vect2' in
-
- if c1.op == Eq && c2.op == Eq
- then exist2 c1.coeffs c2.coeffs
- else None
+ if c1.op == Eq && c2.op == Eq
+ then Vect.exists2 (fun n1 n2 ->
+ Int.equal (compare_big_int (gcd_big_int (numerator n1) (numerator n2)) unit_big_int) 0)
+ c1.coeffs c2.coeffs
+ else None
let extract2 pred l =
- let rec xextract2 rl l =
- match l with
- | [] -> (None,rl) (* Did not find *)
- | e::l ->
- match extract (pred e) l with
- | None,_ -> xextract2 (e::rl) l
- | Some (r,e'),l' -> Some (r,e,e'), List.rev_append rl l' in
-
- xextract2 [] l
+ let rec xextract2 rl l =
+ match l with
+ | [] -> (None,rl) (* Did not find *)
+ | e::l ->
+ match extract (pred e) l with
+ | None,_ -> xextract2 (e::rl) l
+ | Some (r,e'),l' -> Some (r,e,e'), List.rev_append rl l' in
+ xextract2 [] l
-let extract_coprime_equation psys =
- extract2 extract_coprime psys
+let extract_coprime_equation psys =
+ extract2 extract_coprime psys
-let apply_and_normalise f psys =
- List.fold_left (fun acc pc' ->
- match f pc' with
- | None -> pc'::acc
- | Some pc' ->
- match check_sat pc' with
- | Tauto -> acc
- | Unsat prf -> raise (FoundProof prf)
- | Cut(c,p) -> (c,p)::acc
- | Normalise (c,p) -> (c,p)::acc
- ) [] psys
-let pivot_sys v pc psys = apply_and_normalise (pivot v pc) psys
+let pivot_sys v pc psys = apply_and_normalise check_int_sat (pivot v pc) psys
let reduce_coprime psys =
- let oeq,sys = extract_coprime_equation psys in
- match oeq with
- | None -> None (* Nothing to do *)
- | Some((v,n1,n2),(c1,p1),(c2,p2) ) ->
- let (l1,l2) = ext_gcd (numerator n1) (numerator n2) in
- let l1' = Big_int l1 and l2' = Big_int l2 in
- let cstr =
- {coeffs = Vect.add (Vect.mul l1' c1.coeffs) (Vect.mul l2' c2.coeffs);
- op = Eq ;
- cst = (l1' */ c1.cst) +/ (l2' */ c2.cst)
- } in
- let prf = add_proof (mul_proof (numerator l1') p1) (mul_proof (numerator l2') p2) in
-
- Some (pivot_sys v (cstr,prf) ((c1,p1)::sys))
+ let oeq,sys = extract_coprime_equation psys in
+ match oeq with
+ | None -> None (* Nothing to do *)
+ | Some((v,n1,n2),(c1,p1),(c2,p2) ) ->
+ let (l1,l2) = ext_gcd (numerator n1) (numerator n2) in
+ let l1' = Big_int l1 and l2' = Big_int l2 in
+ let cstr =
+ {coeffs = Vect.add (Vect.mul l1' c1.coeffs) (Vect.mul l2' c2.coeffs);
+ op = Eq ;
+ cst = (l1' */ c1.cst) +/ (l2' */ c2.cst)
+ } in
+ let prf = add_proof (mul_cst_proof l1' p1) (mul_cst_proof l2' p2) in
+
+ Some (pivot_sys v (cstr,prf) ((c1,p1)::sys))
(** If there is an equation [eq] of the form 1.x + e = c, do a pivot over x with equation [eq] *)
let reduce_unary psys =
- let is_unary_equation (cstr,prf) =
- if cstr.op == Eq
- then
- try
- Some (fst (List.find (fun (_,n) -> n =/ (Int 1) || n=/ (Int (-1))) cstr.coeffs))
- with Not_found -> None
- else None in
-
- let (oeq,sys) = extract is_unary_equation psys in
- match oeq with
- | None -> None (* Nothing to do *)
- | Some(v,pc) ->
- Some(pivot_sys v pc sys)
-
-let reduce_non_lin_unary psys =
-
- let is_unary_equation (cstr,prf) =
- if cstr.op == Eq
- then
- try
- let x = fst (List.find (fun (x,n) -> (n =/ (Int 1) || n=/ (Int (-1))) && Monomial.is_var (LinPoly.MonT.retrieve x) ) cstr.coeffs) in
- let x' = LinPoly.MonT.retrieve x in
- if List.for_all (fun (y,_) -> Pervasives.(=) y x || Int.equal (snd (Monomial.div (LinPoly.MonT.retrieve y) x')) 0) cstr.coeffs
- then Some x
- else None
- with Not_found -> None
- else None in
-
-
- let (oeq,sys) = extract is_unary_equation psys in
- match oeq with
- | None -> None (* Nothing to do *)
- | Some(v,pc) ->
- Some(apply_and_normalise (LinPoly.pivot_eq v pc) sys)
+ let is_unary_equation (cstr,prf) =
+ if cstr.op == Eq
+ then
+ Vect.find (fun v n -> if n =/ (Int 1) || n=/ (Int (-1)) then Some v else None) cstr.coeffs
+ else None in
+
+ let (oeq,sys) = extract is_unary_equation psys in
+ match oeq with
+ | None -> None (* Nothing to do *)
+ | Some(v,pc) ->
+ Some(pivot_sys v pc sys)
+
let reduce_var_change psys =
- let rec rel_prime vect =
- match vect with
- | [] -> None
- | (x,v)::vect ->
- let v = numerator v in
- try
- let (x',v') = List.find (fun (_,v') ->
- let v' = numerator v' in
- eq_big_int (gcd_big_int v v') unit_big_int) vect in
- Some ((x,v),(x',numerator v'))
- with Not_found -> rel_prime vect in
-
- let rel_prime (cstr,prf) = if cstr.op == Eq then rel_prime cstr.coeffs else None in
-
- let (oeq,sys) = extract rel_prime psys in
-
- match oeq with
- | None -> None
- | Some(((x,v),(x',v')),(c,p)) ->
- let (l1,l2) = ext_gcd v v' in
- let l1,l2 = Big_int l1 , Big_int l2 in
-
- let get v vect =
- match Vect.get v vect with
- | None -> Int 0
- | Some n -> n in
-
- let pivot_eq (c',p') =
- let {coeffs = coeffs ; op = op ; cst = cst} = c' in
- let vx = get x coeffs in
- let vx' = get x' coeffs in
- let m = minus_num (vx */ l1 +/ vx' */ l2) in
- Some ({coeffs =
- Vect.add (Vect.mul m c.coeffs) coeffs ; op = op ; cst = m */ c.cst +/ cst} ,
- AddPrf(MulC(([], m),p),p')) in
-
- Some (apply_and_normalise pivot_eq sys)
-
-let iterate_until_stable f x =
- let rec iter x =
- match f x with
- | None -> x
- | Some x' -> iter x' in
- iter x
-
-let rec app_funs l x =
- match l with
- | [] -> None
- | f::fl ->
- match f x with
- | None -> app_funs fl x
- | Some x' -> Some x'
+ let rec rel_prime vect =
+ match Vect.choose vect with
+ | None -> None
+ | Some(x,v,vect) ->
+ let v = numerator v in
+ match Vect.find (fun x' v' ->
+ let v' = numerator v' in
+ if eq_big_int (gcd_big_int v v') unit_big_int
+ then Some(x',v') else None) vect with
+ | Some(x',v') -> Some ((x,v),(x', v'))
+ | None -> rel_prime vect in
+
+ let rel_prime (cstr,prf) = if cstr.op == Eq then rel_prime cstr.coeffs else None in
+
+ let (oeq,sys) = extract rel_prime psys in
+
+ match oeq with
+ | None -> None
+ | Some(((x,v),(x',v')),(c,p)) ->
+ let (l1,l2) = ext_gcd v v' in
+ let l1,l2 = Big_int l1 , Big_int l2 in
+
+
+ let pivot_eq (c',p') =
+ let {coeffs = coeffs ; op = op ; cst = cst} = c' in
+ let vx = Vect.get x coeffs in
+ let vx' = Vect.get x' coeffs in
+ let m = minus_num (vx */ l1 +/ vx' */ l2) in
+ Some ({coeffs =
+ Vect.add (Vect.mul m c.coeffs) coeffs ; op = op ; cst = m */ c.cst +/ cst} ,
+ AddPrf(MulC((LinPoly.constant m),p),p')) in
+
+ Some (apply_and_normalise check_int_sat pivot_eq sys)
+
let reduction_equations psys =
- iterate_until_stable (app_funs
- [reduce_unary ; reduce_coprime ;
- reduce_var_change (*; reduce_pivot*)]) psys
+ iterate_until_stable (app_funs
+ [reduce_unary ; reduce_coprime ;
+ reduce_var_change (*; reduce_pivot*)]) psys
-let reduction_non_lin_equations psys =
- iterate_until_stable (app_funs
- [reduce_non_lin_unary (*; reduce_coprime ;
- reduce_var_change ; reduce_pivot *)]) psys
- (** [get_bound sys] returns upon success an interval (lb,e,ub) with proofs *)
+(** [get_bound sys] returns upon success an interval (lb,e,ub) with proofs *)
let get_bound sys =
- let is_small (v,i) =
- match Itv.range i with
- | None -> false
- | Some i -> i <=/ (Int 1) in
-
- let select_best (x1,i1) (x2,i2) =
- if Itv.smaller_itv i1 i2
- then (x1,i1) else (x2,i2) in
-
- (* For lia, there are no equations => these precautions are not needed *)
- (* For nlia, there are equations => do not enumerate over equations! *)
- let all_planes sys =
- let (eq,ineq) = List.partition (fun c -> c.op == Eq) sys in
- match eq with
- | [] -> List.rev_map (fun c -> c.coeffs) ineq
- | _ ->
- List.fold_left (fun acc c ->
- if List.exists (fun c' -> Vect.equal c.coeffs c'.coeffs) eq
- then acc else c.coeffs ::acc) [] ineq in
-
- let smallest_interval =
- List.fold_left
- (fun acc vect ->
- if is_small acc
- then acc
- else
- match Fourier.optimise vect sys with
- | None -> acc
- | Some i ->
- if debug then Printf.printf "Found a new bound %a" Vect.pp_vect vect ;
- select_best (vect,i) acc) (Vect.null, (None,None)) (all_planes sys) in
- let smallest_interval =
- match smallest_interval
- with
- | (x,(Some i, Some j)) -> Some(i,x,j)
- | x -> None (* This should not be possible *)
- in
- match smallest_interval with
- | Some (lb,e,ub) ->
- let (lbn,lbd) = (sub_big_int (numerator lb) unit_big_int, denominator lb) in
- let (ubn,ubd) = (add_big_int unit_big_int (numerator ub) , denominator ub) in
- (match
- (* x <= ub -> x > ub *)
- xlinear_prover ({coeffs = Vect.mul (Big_int ubd) e ; op = Ge ; cst = Big_int ubn} :: sys),
- (* lb <= x -> lb > x *)
- xlinear_prover
- ({coeffs = Vect.mul (minus_num (Big_int lbd)) e ; op = Ge ; cst = minus_num (Big_int lbn)} :: sys)
- with
- | Some cub , Some clb -> Some(List.tl clb,(lb,e,ub), List.tl cub)
- | _ -> failwith "Interval without proof"
- )
- | None -> None
+ let is_small (v,i) =
+ match Itv.range i with
+ | None -> false
+ | Some i -> i <=/ (Int 1) in
+
+ let select_best (x1,i1) (x2,i2) =
+ if Itv.smaller_itv i1 i2
+ then (x1,i1) else (x2,i2) in
+
+ (* For lia, there are no equations => these precautions are not needed *)
+ (* For nlia, there are equations => do not enumerate over equations! *)
+ let all_planes sys =
+ let (eq,ineq) = List.partition (fun c -> c.op == Eq) sys in
+ match eq with
+ | [] -> List.rev_map (fun c -> c.coeffs) ineq
+ | _ ->
+ List.fold_left (fun acc c ->
+ if List.exists (fun c' -> Vect.equal c.coeffs c'.coeffs) eq
+ then acc else c.coeffs ::acc) [] ineq in
+
+ let smallest_interval =
+ List.fold_left
+ (fun acc vect ->
+ if is_small acc
+ then acc
+ else
+ match optimise vect sys with
+ | None -> acc
+ | Some i ->
+ if debug then Printf.printf "Found a new bound %a in %a" Vect.pp vect Itv.pp i;
+ select_best (vect,i) acc) (Vect.null, (None,None)) (all_planes sys) in
+ let smallest_interval =
+ match smallest_interval
+ with
+ | (x,(Some i, Some j)) -> Some(i,x,j)
+ | x -> None (* This should not be possible *)
+ in
+ match smallest_interval with
+ | Some (lb,e,ub) ->
+ let (lbn,lbd) = (sub_big_int (numerator lb) unit_big_int, denominator lb) in
+ let (ubn,ubd) = (add_big_int unit_big_int (numerator ub) , denominator ub) in
+ (match
+ (* x <= ub -> x > ub *)
+ direct_linear_prover ({coeffs = Vect.mul (Big_int ubd) e ; op = Ge ; cst = Big_int ubn} :: sys),
+ (* lb <= x -> lb > x *)
+ direct_linear_prover
+ ({coeffs = Vect.mul (minus_num (Big_int lbd)) e ; op = Ge ; cst = minus_num (Big_int lbn)} :: sys)
+ with
+ | Some cub , Some clb -> Some(List.tl (Vect.to_list clb),(lb,e,ub), List.tl (Vect.to_list cub))
+ | _ -> failwith "Interval without proof"
+ )
+ | None -> None
let check_sys sys =
- List.for_all (fun (c,p) -> List.for_all (fun (_,n) -> sign_num n <> 0) c.coeffs) sys
+ List.for_all (fun (c,p) -> Vect.for_all (fun _ n -> sign_num n <> 0) c.coeffs) sys
let xlia (can_enum:bool) reduction_equations sys =
-
- let rec enum_proof (id:int) (sys:prf_sys) : proof option =
- if debug then (Printf.printf "enum_proof\n" ; flush stdout) ;
- assert (check_sys sys) ;
-
- let nsys,prf = List.split sys in
- match get_bound nsys with
- | None -> None (* Is the systeme really unbounded ? *)
- | Some(prf1,(lb,e,ub),prf2) ->
- if debug then Printf.printf "Found interval: %a in [%s;%s] -> " Vect.pp_vect e (string_of_num lb) (string_of_num ub) ;
- (match start_enum id e (ceiling_num lb) (floor_num ub) sys
- with
- | Some prfl ->
- Some(Enum(id,proof_of_farkas prf prf1,e, proof_of_farkas prf prf2,prfl))
- | None -> None
- )
- and start_enum id e clb cub sys =
- if clb >/ cub
- then Some []
- else
- let eq = {coeffs = e ; op = Eq ; cst = clb} in
- match aux_lia (id+1) ((eq, Def id) :: sys) with
- | None -> None
- | Some prf ->
- match start_enum id e (clb +/ (Int 1)) cub sys with
- | None -> None
- | Some l -> Some (prf::l)
-
- and aux_lia (id:int) (sys:prf_sys) : proof option =
- assert (check_sys sys) ;
- if debug then Printf.printf "xlia: %a \n" (pp_list (fun o (c,_) -> output_cstr o c)) sys ;
- try
- let sys = reduction_equations sys in
- if debug then
- Printf.printf "after reduction: %a \n" (pp_list (fun o (c,_) -> output_cstr o c)) sys ;
- match linear_prover sys with
- | Some prf -> Some (Step(id,prf,Done))
- | None -> if can_enum then enum_proof id sys else None
- with FoundProof prf ->
+ let rec enum_proof (id:int) (sys:prf_sys) : ProofFormat.proof option =
+ if debug then (Printf.printf "enum_proof\n" ; flush stdout) ;
+ assert (check_sys sys) ;
+
+ let nsys,prf = List.split sys in
+ match get_bound nsys with
+ | None -> None (* Is the systeme really unbounded ? *)
+ | Some(prf1,(lb,e,ub),prf2) ->
+ if debug then Printf.printf "Found interval: %a in [%s;%s] -> " Vect.pp e (string_of_num lb) (string_of_num ub) ;
+ (match start_enum id e (ceiling_num lb) (floor_num ub) sys
+ with
+ | Some prfl ->
+ Some(Enum(id,proof_of_farkas (env_of_list prf) (Vect.from_list prf1),e,
+ proof_of_farkas (env_of_list prf) (Vect.from_list prf2),prfl))
+ | None -> None
+ )
+
+ and start_enum id e clb cub sys =
+ if clb >/ cub
+ then Some []
+ else
+ let eq = {coeffs = e ; op = Eq ; cst = clb} in
+ match aux_lia (id+1) ((eq, Def id) :: sys) with
+ | None -> None
+ | Some prf ->
+ match start_enum id e (clb +/ (Int 1)) cub sys with
+ | None -> None
+ | Some l -> Some (prf::l)
+
+ and aux_lia (id:int) (sys:prf_sys) : ProofFormat.proof option =
+ assert (check_sys sys) ;
+ if debug then Printf.printf "xlia: %a \n" (pp_list ";" (fun o (c,_) -> output_cstr o c)) sys ;
+ try
+ let sys = reduction_equations sys in
+ if debug then
+ Printf.printf "after reduction: %a \n" (pp_list ";" (fun o (c,_) -> output_cstr o c)) sys ;
+ match linear_prover_cstr sys with
+ | Some prf -> Some (Step(id,prf,Done))
+ | None -> if can_enum then enum_proof id sys else None
+ with FoundProof prf ->
(* [reduction_equations] can find a proof *)
- Some(Step(id,prf,Done)) in
+ Some(Step(id,prf,Done)) in
(* let sys' = List.map (fun (p,o) -> Mc.norm0 p , o) sys in*)
- let id = List.length sys in
- let orpf =
- try
- let sys = simpl_sys sys in
- aux_lia id sys
- with FoundProof pr -> Some(Step(id,pr,Done)) in
- match orpf with
- | None -> None
- | Some prf ->
- (*Printf.printf "direct proof %a\n" output_proof prf ; *)
- let env = List.mapi (fun i _ -> i) sys in
- let prf = compile_proof env prf in
- (*try
+ let id = 1 + (List.fold_left
+ (fun acc (_,r) -> max acc (ProofFormat.pr_rule_max_id r)) 0 sys) in
+ let orpf =
+ try
+ let sys = simpl_sys sys in
+ aux_lia id sys
+ with FoundProof pr -> Some(Step(id,pr,Done)) in
+ match orpf with
+ | None -> None
+ | Some prf ->
+ let env = CList.interval 0 (id - 1) in
+ if debug then begin
+ Printf.fprintf stdout "direct proof %a\n" output_proof prf;
+ flush stdout;
+ end;
+ let prf = compile_proof env prf in
+ (*try
if Mc.zChecker sys' prf then Some prf else
raise Certificate.BadCertificate
with Failure s -> (Printf.printf "%s" s ; Some prf)
- *) Some prf
-
-
-let cstr_compat_of_poly (p,o) =
- let (v,c) = LinPoly.linpol_of_pol p in
- {coeffs = v ; op = o ; cst = minus_num c }
-
+ *) Some prf
+
+let xlia_simplex env sys =
+ match Simplex.integer_solver sys with
+ | None -> None
+ | Some prf ->
+ (*let _ = ProofFormat.eval_proof (env_of_list env) prf in *)
+
+ let id = 1 + (List.fold_left
+ (fun acc (_,r) -> max acc (ProofFormat.pr_rule_max_id r)) 0 sys) in
+ let env = CList.interval 0 (id - 1) in
+ Some (compile_proof env prf)
+
+let xlia env0 en red sys =
+ if !use_simplex then xlia_simplex env0 sys
+ else xlia en red sys
+
+
+let dump_file = ref None
+
+let gen_bench (tac, prover) can_enum prfdepth sys =
+ let res = prover can_enum prfdepth sys in
+ (match !dump_file with
+ | None -> ()
+ | Some file ->
+ begin
+ let o = open_out (Filename.temp_file ~temp_dir:(Sys.getcwd ()) file ".v") in
+ let sys = develop_constraints prfdepth z_spec sys in
+ Printf.fprintf o "Require Import ZArith Lia. Open Scope Z_scope.\n";
+ Printf.fprintf o "Goal %a.\n" (LinPoly.pp_goal "Z") (List.map fst sys) ;
+ begin
+ match res with
+ | None ->
+ Printf.fprintf o "Proof.\n intros. Fail %s.\nAbort.\n" tac
+ | Some res ->
+ Printf.fprintf o "Proof.\n intros. %s.\nQed.\n" tac
+ end
+ ;
+ flush o ;
+ close_out o ;
+ end);
+ res
let lia (can_enum:bool) (prfdepth:int) sys =
- LinPoly.MonT.clear ();
- max_nb_cstr := compute_max_nb_cstr sys prfdepth ;
- let sys = List.map (develop_constraint z_spec) sys in
- let (sys:cstr_compat list) = List.map cstr_compat_of_poly sys in
- let sys = List.mapi (fun i c -> (c,Hyp i)) sys in
- xlia can_enum reduction_equations sys
+ let sys = develop_constraints prfdepth z_spec sys in
+ if debug then begin
+ Printf.fprintf stdout "Input problem\n";
+ List.iter (fun s -> Printf.fprintf stdout "%a\n" WithProof.output s) sys;
+ end;
+
+ let sys' = List.map (fun ((p,o),prf) -> (cstr_of_poly (p,o), prf)) sys in
+ xlia (List.map fst sys) can_enum reduction_equations sys'
+let make_cstr_system sys =
+ List.map (fun ((p,o),prf) -> (cstr_of_poly (p,o), prf)) sys
let nlia enum prfdepth sys =
- LinPoly.MonT.clear ();
- max_nb_cstr := compute_max_nb_cstr sys prfdepth;
- let sys = List.map (develop_constraint z_spec) sys in
- let sys = List.mapi (fun i c -> (c,Hyp i)) sys in
-
- let is_linear = List.for_all (fun ((p,_),_) -> Poly.is_linear p) sys in
-
- let collect_square =
- List.fold_left (fun acc ((p,_),_) -> Poly.fold
- (fun m _ acc ->
- match Monomial.sqrt m with
- | None -> acc
- | Some s -> MonMap.add s m acc) p acc) MonMap.empty sys in
- let sys = MonMap.fold (fun s m acc ->
- let s = LinPoly.linpol_of_pol (Poly.add s (Int 1) (Poly.constant (Int 0))) in
- let m = Poly.add m (Int 1) (Poly.constant (Int 0)) in
- ((m, Ge), (Square s))::acc) collect_square sys in
-
- (* List.iter (fun ((p,_),_) -> Printf.printf "square %a\n" Poly.pp p) gen_square*)
-
- let sys =
- if is_linear then sys
- else sys @ (all_sym_pairs (fun ((c,o),p) ((c',o'),p') ->
- ((Poly.product c c',opMult o o'), MulPrf(p,p'))) sys) in
+ let sys = develop_constraints prfdepth z_spec sys in
+ let is_linear = List.for_all (fun ((p,_),_) -> LinPoly.is_linear p) sys in
+
+ if debug then begin
+ Printf.fprintf stdout "Input problem\n";
+ List.iter (fun s -> Printf.fprintf stdout "%a\n" WithProof.output s) sys;
+ end;
+
+ if is_linear
+ then xlia (List.map fst sys) enum reduction_equations (make_cstr_system sys)
+ else
+ (*
+ let sys1 = elim_every_substitution sys in
+ No: if a wrong equation is chosen, the proof may fail.
+ It would only be safe if the variable is linear...
+ *)
+ let sys1 = elim_simple_linear_equality sys in
+ let sys2 = saturate_linear_equality_non_linear sys1 in
+ let sys3 = nlinear_preprocess (sys1@sys2) in
+
+ let sys4 = make_cstr_system ((*sys2@*)sys3) in
+ (* [reduction_equations] is too brutal - there should be some non-linear reasoning *)
+ xlia (List.map fst sys) enum reduction_equations sys4
+
+(* For regression testing, if bench = true generate a Coq goal *)
+
+let lia can_enum prfdepth sys =
+ gen_bench ("lia",lia) can_enum prfdepth sys
+
+let nlia enum prfdepth sys =
+ gen_bench ("nia",nlia) enum prfdepth sys
+
- let sys = List.map (fun (c,p) -> cstr_compat_of_poly c,p) sys in
- assert (check_sys sys) ;
- xlia enum (if is_linear then reduction_equations else reduction_non_lin_equations) sys
diff --git a/plugins/micromega/certificate.mli b/plugins/micromega/certificate.mli
index 13d50d1eee..e925f1bc5e 100644
--- a/plugins/micromega/certificate.mli
+++ b/plugins/micromega/certificate.mli
@@ -10,13 +10,33 @@
module Mc = Micromega
-type 'a number_spec
+(** [use_simplex] is bound to the Coq option Simplex.
+ If set, use the Simplex method, otherwise use Fourier *)
+val use_simplex : bool ref
+
+(** [dump_file] is bound to the Coq option Dump Arith.
+ If set to some [file], arithmetic goals are dumped in filexxx.v *)
+val dump_file : string option ref
+
+(** [q_cert_of_pos prf] converts a Sos proof into a rational Coq proof *)
val q_cert_of_pos : Sos_types.positivstellensatz -> Mc.q Mc.psatz
+
+(** [z_cert_of_pos prf] converts a Sos proof into an integer Coq proof *)
val z_cert_of_pos : Sos_types.positivstellensatz -> Mc.z Mc.psatz
+
+(** [lia enum depth sys] generates an unsat proof for the linear constraints in [sys].
+ If the Simplex option is set, any failure to find a proof should be considered as a bug. *)
val lia : bool -> int -> (Mc.z Mc.pExpr * Mc.op1) list -> Mc.zArithProof option
+
+(** [nlia enum depth sys] generates an unsat proof for the non-linear constraints in [sys].
+ The solver is incomplete -- the problem is undecidable *)
val nlia : bool -> int -> (Mc.z Mc.pExpr * Mc.op1) list -> Mc.zArithProof option
+
+(** [linear_prover_with_cert depth sys] generates an unsat proof for the linear constraints in [sys].
+ Over the rationals, the solver is complete. *)
+val linear_prover_with_cert : int -> (Mc.q Mc.pExpr * Mc.op1) list -> Mc.q Micromega.psatz option
+
+(** [nlinear depth sys] generates an unsat proof for the non-linear constraints in [sys].
+ The solver is incompete -- the problem is decidable. *)
val nlinear_prover : int -> (Mc.q Mc.pExpr * Mc.op1) list -> Mc.q Mc.psatz option
-val linear_prover_with_cert : int -> 'a number_spec ->
- ('a Mc.pExpr * Mc.op1) list -> 'a Mc.psatz option
-val q_spec : Mc.q number_spec
diff --git a/plugins/micromega/coq_micromega.ml b/plugins/micromega/coq_micromega.ml
index e0a369ca5f..ec080edbdb 100644
--- a/plugins/micromega/coq_micromega.ml
+++ b/plugins/micromega/coq_micromega.ml
@@ -12,7 +12,7 @@
(* *)
(* ** Toplevel definition of tactics ** *)
(* *)
-(* - Modules ISet, M, Mc, Env, Cache, CacheZ *)
+(* - Modules M, Mc, Env, Cache, CacheZ *)
(* *)
(* Frédéric Besson (Irisa/Inria) 2006-20011 *)
(* *)
@@ -44,7 +44,7 @@ let lia_enum = ref true
let lia_proof_depth = ref max_depth
let get_lia_option () =
- (!lia_enum,!lia_proof_depth)
+ (!Certificate.use_simplex,!lia_enum,!lia_proof_depth)
let get_lra_option () =
!lra_proof_depth
@@ -70,10 +70,32 @@ let _ =
optread = (fun () -> !lia_enum);
optwrite = (fun x -> lia_enum := x)
} in
+
+ let solver_opt =
+ {
+ optdepr = false;
+ optname = "Use the Simplex instead of Fourier elimination";
+ optkey = ["Simplex"];
+ optread = (fun () -> !Certificate.use_simplex);
+ optwrite = (fun x -> Certificate.use_simplex := x)
+ } in
+
+ let dump_file_opt =
+ {
+ optdepr = false;
+ optname = "Generate Coq goals in file from calls to 'lia' 'nia'";
+ optkey = ["Dump"; "Arith"];
+ optread = (fun () -> !Certificate.dump_file);
+ optwrite = (fun x -> Certificate.dump_file := x)
+ } in
+
+ let _ = declare_bool_option solver_opt in
+ let _ = declare_stringopt_option dump_file_opt in
let _ = declare_int_option (int_opt ["Lra"; "Depth"] lra_proof_depth) in
let _ = declare_int_option (int_opt ["Lia"; "Depth"] lia_proof_depth) in
let _ = declare_bool_option lia_enum_opt in
()
+
(**
* Initialize a tag type to the Tag module declaration (see Mutils).
@@ -288,11 +310,6 @@ let rec add_term t0 = function
xcnf true f
-(**
- * MODULE: Ordered set of integers.
- *)
-
-module ISet = Set.Make(Int)
(**
* Given a set of integers s=\{i0,...,iN\} and a list m, return the list of
@@ -1937,7 +1954,9 @@ let really_call_csdpcert : provername -> micromega_polys -> Sos_types.positivste
["plugins"; "micromega"; "csdpcert" ^ Coq_config.exec_extension] in
match ((command cmdname [|cmdname|] (provername,poly)) : csdp_certificate) with
- | F str -> failwith str
+ | F str ->
+ if debug then Printf.fprintf stdout "really_call_csdpcert : %s\n" str;
+ raise (failwith str)
| S res -> res
(**
@@ -2047,7 +2066,7 @@ let compact_pt pt f =
let lift_pexpr_prover p l = p (List.map (fun (e,o) -> Mc.denorm e , o) l)
module CacheZ = PHashtable(struct
- type prover_option = bool * int
+ type prover_option = bool * bool* int
type t = prover_option * ((Mc.z Mc.pol * Mc.op1) list)
let equal = (=)
@@ -2060,8 +2079,8 @@ module CacheQ = PHashtable(struct
let hash = Hashtbl.hash
end)
-let memo_zlinear_prover = CacheZ.memo ".lia.cache" (fun ((ce,b),s) -> lift_pexpr_prover (Certificate.lia ce b) s)
-let memo_nlia = CacheZ.memo ".nia.cache" (fun ((ce,b),s) -> lift_pexpr_prover (Certificate.nlia ce b) s)
+let memo_zlinear_prover = CacheZ.memo ".lia.cache" (fun ((_,ce,b),s) -> lift_pexpr_prover (Certificate.lia ce b) s)
+let memo_nlia = CacheZ.memo ".nia.cache" (fun ((_,ce,b),s) -> lift_pexpr_prover (Certificate.nlia ce b) s)
let memo_nra = CacheQ.memo ".nra.cache" (fun (o,s) -> lift_pexpr_prover (Certificate.nlinear_prover o) s)
@@ -2069,7 +2088,7 @@ let memo_nra = CacheQ.memo ".nra.cache" (fun (o,s) -> lift_pexpr_prover (Certifi
let linear_prover_Q = {
name = "linear prover";
get_option = get_lra_option ;
- prover = (fun (o,l) -> lift_pexpr_prover (Certificate.linear_prover_with_cert o Certificate.q_spec) l) ;
+ prover = (fun (o,l) -> lift_pexpr_prover (Certificate.linear_prover_with_cert o ) l) ;
hyps = hyps_of_cone ;
compact = compact_cone ;
pp_prf = pp_psatz pp_q ;
@@ -2080,7 +2099,7 @@ let linear_prover_Q = {
let linear_prover_R = {
name = "linear prover";
get_option = get_lra_option ;
- prover = (fun (o,l) -> lift_pexpr_prover (Certificate.linear_prover_with_cert o Certificate.q_spec) l) ;
+ prover = (fun (o,l) -> lift_pexpr_prover (Certificate.linear_prover_with_cert o ) l) ;
hyps = hyps_of_cone ;
compact = compact_cone ;
pp_prf = pp_psatz pp_q ;
diff --git a/plugins/micromega/itv.ml b/plugins/micromega/itv.ml
new file mode 100644
index 0000000000..dc1df7ec9f
--- /dev/null
+++ b/plugins/micromega/itv.ml
@@ -0,0 +1,80 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+
+(** Intervals (extracted from mfourier.ml) *)
+
+open Num
+ (** 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 explicitly normalised.
+ *)
+
+ let pp o (n1,n2) =
+ (match n1 with
+ | None -> output_string o "]-oo"
+ | Some n -> Printf.fprintf o "[%s" (string_of_num n)
+ );
+ output_string o ",";
+ (match n2 with
+ | None -> output_string o "+oo["
+ | Some n -> Printf.fprintf o "%s]" (string_of_num n)
+ )
+
+
+
+ (** 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
+
+(** [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 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
diff --git a/plugins/micromega/itv.mli b/plugins/micromega/itv.mli
new file mode 100644
index 0000000000..31f6a89fe2
--- /dev/null
+++ b/plugins/micromega/itv.mli
@@ -0,0 +1,18 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+open Num
+
+type interval = num option * num option
+val pp : out_channel -> interval -> unit
+val inter : interval -> interval -> interval option
+val range : interval -> num option
+val smaller_itv : interval -> interval -> bool
+val in_bound : interval -> num -> bool
+val norm_itv : interval -> interval option
diff --git a/plugins/micromega/mfourier.ml b/plugins/micromega/mfourier.ml
index 3328abdab7..baf8c82355 100644
--- a/plugins/micromega/mfourier.ml
+++ b/plugins/micromega/mfourier.ml
@@ -1,3 +1,13 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+
open Util
open Num
open Polynomial
@@ -8,66 +18,6 @@ let debug = false
let compare_float (p : float) q = Pervasives.compare p q
(** Implementation of intervals *)
-module Itv =
-struct
-
- (** 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 explicitly normalised.
- *)
-
- (** 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
-
-(** [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 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 Itv
type vector = Vect.t
@@ -84,8 +34,6 @@ type proof =
| Elim of var * proof * proof
| And of proof * proof
-let max_nb_cstr = ref max_int
-
type system = {
sys : cstr_info ref System.t ;
vars : ISet.t
@@ -126,7 +74,7 @@ let pp_cstr o (vect,bnd) =
| None -> ()
| Some n -> Printf.fprintf o "%s <= " (string_of_num n))
;
- pp_vect o vect ;
+ Vect.pp o vect ;
(match r with
| None -> output_string o"\n"
| Some n -> Printf.fprintf o "<=%s\n" (string_of_num n))
@@ -185,30 +133,23 @@ 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),
+ match Vect.choose vect with
+ | None -> if Itv.in_bound (l,r) (Int 0) then Redundant else Contradiction
+ | Some (_,n,_) -> Cstr(Vect.div n vect,
let divn x = x // n in
if Int.equal (sign_num n) 1
then{cinfo with bound = (Option.map divn l , Option.map divn r) }
else {cinfo with pos = cinfo.neg ; neg = cinfo.pos ; bound = (Option.map divn r , Option.map divn l)})
-(** For compatibility, there is an external representation of constraints *)
+(** For compatibility, there is an external representation of constraints *)
-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 Int.equal sg 1 then count n (p+1) v else count (n+1) p v in
- count 0 0 v
+ Vect.fold (fun (n,p) _ vl ->
+ let sg = sign_num vl in
+ assert (sg <> 0) ;
+ if Int.equal sg 1 then (n,p+1)else (n+1, p)) (0,0) v
let norm_cstr {coeffs = v ; op = o ; cst = c} idx =
@@ -217,7 +158,9 @@ let norm_cstr {coeffs = v ; op = o ; cst = c} idx =
normalise_cstr v {pos = p ; neg = n ; bound =
(match o with
| Eq -> Some c , Some c
- | Ge -> Some c , None) ;
+ | Ge -> Some c , None
+ | Gt -> raise Polynomial.Strict
+ ) ;
prf = Assum idx }
@@ -237,7 +180,7 @@ let load_system l =
| 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
+ Vect.fold (fun s v _ -> ISet.add v s) vrs cstr.coeffs) ISet.empty li in
{sys = sys ;vars = vars}
@@ -255,27 +198,7 @@ let system_list sys =
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 Int.equal 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 res = xadd v1 v2 in
+ let res = mul_add (Int 1 // c1) v1 (Int 1 // c2) v2 in
(res, count res)
let add (v1,c1) (v2,c2) =
@@ -294,9 +217,9 @@ let add (v1,c1) (v2,c2) =
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 *)
+ | Int 0 -> (* The constraint does not mention [x], store it in m *)
(l,(vect,info)::m,r)
- | Some vl -> (* otherwise *)
+ | vl -> (* otherwise *)
let cons_bound lst bd =
match bd with
@@ -352,7 +275,8 @@ let project vr sys =
let project_using_eq vr c vect bound prf (vect',info') =
match get vr vect' with
- | Some c2 ->
+ | Int 0 -> (vect',info')
+ | c2 ->
let c1 = if c2 >=/ Int 0 then minus_num c else c in
let c2 = abs_num c2 in
@@ -367,10 +291,10 @@ let project_using_eq vr c vect bound prf (vect',info') =
(Option.map f l , Option.map 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 = Option.get (get vr vect) in
+ let c = get vr vect in
let elim_var = project_using_eq vr c vect cst prf in
@@ -397,16 +321,13 @@ module IMap = CMap.Make(Int)
The function returns as second argument, a sub-vector consisting in the variables that are not in [map]. *)
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) []
+ Vect.fold (fun (sum,rst) v vl ->
+ try
+ let val_v = IMap.find v map in
+ (sum +/ (val_v */ vl), rst)
+ with
+ Not_found -> (sum, Vect.set v vl rst)) (Int 0,Vect.null) vect
+
(** [restrict_bound n sum itv] returns the interval of [x]
@@ -427,11 +348,13 @@ let restrict_bound n sum (itv:interval) =
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
+ let vl = Vect.get v rst in
match inter bnd (restrict_bound vl sum (!iref).bound) with
- | None -> failwith "bound_of_variable: impossible"
+ | None ->
+ Printf.fprintf stdout "bound_of_variable: eval_vecr %a = %s,%a\n"
+ Vect.pp vect (Num.string_of_num sum) Vect.pp rst ;
+ Printf.fprintf stdout "current interval: %a\n" Itv.pp (!iref).bound;
+ failwith "bound_of_variable: impossible"
| Some itv -> itv) sys (None,None)
@@ -458,12 +381,13 @@ 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);
+ if debug then Printf.printf "solve_sys :\n %a" pp_system sys.sys ;
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 ;
+ (Printf.printf "\nE %a = %s variable %i\n" Vect.pp 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)
@@ -503,9 +427,9 @@ struct
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 ->
+ match Vect.choose l1 with
+ | None -> xpart rl ((Vect.null,info)::ltl) n (info.neg+info.pos+z) p
+ | Some(vr, vl, rl1) ->
if Int.equal v vr
then
let cons_bound lst bd =
@@ -557,24 +481,26 @@ struct
| _ -> false
let rec unroll_until v l =
- match l with
- | [] -> (false,[])
- | (i,_)::rl -> if Int.equal i v
+ match Vect.choose l with
+ | None -> (false,Vect.null)
+ | Some(i,_,rl) -> if Int.equal i v
then (true,rl)
else if i < v then unroll_until v rl else (false,l)
+
let rec choose_simple_equation eqs =
match eqs with
| [] -> None
| (vect,a,prf,ln)::eqs ->
- match vect with
- | [i,_] -> Some (i,vect,a,prf,ln)
- | _ -> choose_simple_equation eqs
+ match Vect.choose vect with
+ | Some(i,v,rst) -> if Vect.is_null rst
+ then Some (i,vect,a,prf,ln)
+ else choose_simple_equation eqs
+ | _ -> choose_simple_equation eqs
-
- let choose_primal_equation eqs sys_l =
+ let choose_primal_equation eqs (sys_l: (Vect.t *cstr_info) list) =
(* Counts the number of equations referring to variable [v] --
It looks like nb_cst is dead...
@@ -586,9 +512,9 @@ struct
else nb_eq) 0 sys_l in
let rec find_var vect =
- match vect with
- | [] -> None
- | (i,_)::vect ->
+ match Vect.choose vect with
+ | None -> None
+ | Some(i,_,vect) ->
let nb_eq = is_primal_equation_var i in
if Int.equal nb_eq 2
then Some i else find_var vect in
@@ -638,9 +564,9 @@ struct
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
+ match Vect.choose eqr with
+ | None -> costs
+ | Some(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
@@ -692,10 +618,10 @@ struct
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
+ let vect = IMap.fold (fun v i vect -> Vect.set v i vect) map Vect.null in
+ if debug then Printf.printf "SOLUTION %a" Vect.pp vect ;
+ let res = Inl vect in
+ res
end
@@ -735,8 +661,8 @@ struct
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 ->
+ | Int 0 , _ | _ , Int 0 -> None
+ | a , b ->
if Int.equal ((sign_num a) * (sign_num b)) (-1)
then
Some (add (p1,abs_num a) (p2,abs_num b) ,
@@ -768,7 +694,7 @@ struct
| Cstr(v,info) -> Inl ((prf,cstr,v,info)::acc)) (Inl []) l
- type oproof = (vector * cstr_compat * num) option
+ type oproof = (vector * cstr * num) option
let merge_proof (oleft:oproof) (prf,cstr,v,info) (oright:oproof) =
let (l,r) = info.bound in
@@ -789,9 +715,9 @@ struct
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 Vect.choose cstrr.coeffs with
+ | None -> Inr (add (prfl,Int 1) (prfr,Int 1), cstrr) (* this is wrong *)
+ | Some(v,_,_) ->
match pivot v (prfl,cstrl) (prfr,cstrr) with
| None -> failwith "merge_proof : pivot is not possible"
| Some x -> Inr x
@@ -804,7 +730,7 @@ let mk_proof hyps prf =
let rec mk_proof prf =
match prf with
- | Assum i -> [ ([i, Int 1] , List.nth hyps i) ]
+ | Assum i -> [ (Vect.set i (Int 1) Vect.null , List.nth hyps i) ]
| Elim(v,prf1,prf2) ->
let prfsl = mk_proof prf1
diff --git a/plugins/micromega/mfourier.mli b/plugins/micromega/mfourier.mli
index f1d8edeab6..45a81cc118 100644
--- a/plugins/micromega/mfourier.mli
+++ b/plugins/micromega/mfourier.mli
@@ -8,25 +8,18 @@
(* * (see LICENSE file for the text of the license) *)
(************************************************************************)
-module Itv : sig
-
- type interval = Num.num option * Num.num option
- val range : interval -> Num.num option
- val smaller_itv : interval -> interval -> bool
-
-end
-
module IMap : CSig.MapS with type key = int
type proof
module Fourier : sig
- val find_point : Polynomial.cstr_compat list ->
- ((IMap.key * Num.num) list, proof) Util.union
- val optimise : Polynomial.Vect.t ->
- Polynomial.cstr_compat list ->
+ val find_point : Polynomial.cstr list ->
+ (Vect.t, proof) Util.union
+
+ val optimise : Vect.t ->
+ Polynomial.cstr list ->
Itv.interval option
end
@@ -35,15 +28,11 @@ val pp_proof : out_channel -> proof -> unit
module Proof : sig
- val mk_proof : Polynomial.cstr_compat list ->
- proof -> (Polynomial.Vect.t * Polynomial.cstr_compat) list
+ val mk_proof : Polynomial.cstr list ->
+ proof -> (Vect.t * Polynomial.cstr) list
val add_op : Polynomial.op -> Polynomial.op -> Polynomial.op
end
-val max_nb_cstr : int ref
-
-val eval_op : Polynomial.op -> Num.num -> Num.num -> bool
-
exception TimeOut
diff --git a/plugins/micromega/micromega.ml b/plugins/micromega/micromega.ml
index 52c6ef983d..f67f1da146 100644
--- a/plugins/micromega/micromega.ml
+++ b/plugins/micromega/micromega.ml
@@ -1484,17 +1484,17 @@ let psub1 =
let padd1 =
padd0 Z0 Z.add zeq_bool
-(** val norm0 : z pExpr -> z pol **)
+(** val normZ : z pExpr -> z pol **)
-let norm0 =
+let normZ =
norm Z0 (Zpos XH) Z.add Z.mul Z.sub Z.opp zeq_bool
(** val xnormalise0 : z formula -> z nFormula list **)
let xnormalise0 t0 =
let { flhs = lhs; fop = o; frhs = rhs } = t0 in
- let lhs0 = norm0 lhs in
- let rhs0 = norm0 rhs in
+ let lhs0 = normZ lhs in
+ let rhs0 = normZ rhs in
(match o with
| OpEq ->
((psub1 lhs0 (padd1 rhs0 (Pc (Zpos XH)))),NonStrict)::(((psub1 rhs0
@@ -1516,8 +1516,8 @@ let normalise t0 =
let xnegate0 t0 =
let { flhs = lhs; fop = o; frhs = rhs } = t0 in
- let lhs0 = norm0 lhs in
- let rhs0 = norm0 rhs in
+ let lhs0 = normZ lhs in
+ let rhs0 = normZ rhs in
(match o with
| OpEq -> ((psub1 lhs0 rhs0),Equal)::[]
| OpNEq ->
@@ -1707,6 +1707,12 @@ let qunsat =
let qdeduce =
nformula_plus_nformula { qnum = Z0; qden = XH } qplus qeq_bool
+(** val normQ : q pExpr -> q pol **)
+
+let normQ =
+ norm { qnum = Z0; qden = XH } { qnum = (Zpos XH); qden = XH } qplus qmult
+ qminus qopp qeq_bool
+
(** val qTautoChecker : q formula bFormula -> qWitness list -> bool **)
let qTautoChecker f w =
diff --git a/plugins/micromega/micromega.mli b/plugins/micromega/micromega.mli
index 9619781786..72c2bf7da3 100644
--- a/plugins/micromega/micromega.mli
+++ b/plugins/micromega/micromega.mli
@@ -151,8 +151,7 @@ val mkPinj : positive -> 'a1 pol -> 'a1 pol
val mkPinj_pred : positive -> 'a1 pol -> 'a1 pol
-val mkPX :
- 'a1 -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
+val mkPX : 'a1 -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
val mkXi : 'a1 -> 'a1 -> positive -> 'a1 pol
@@ -164,49 +163,27 @@ val paddC : ('a1 -> 'a1 -> 'a1) -> 'a1 pol -> 'a1 -> 'a1 pol
val psubC : ('a1 -> 'a1 -> 'a1) -> 'a1 pol -> 'a1 -> 'a1 pol
-val paddI :
- ('a1 -> 'a1 -> 'a1) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol ->
- positive -> 'a1 pol -> 'a1 pol
+val paddI : ('a1 -> 'a1 -> 'a1) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
-val psubI :
- ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 pol -> 'a1 pol -> 'a1 pol) ->
- 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
+val psubI : ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
-val paddX :
- 'a1 -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol
- -> positive -> 'a1 pol -> 'a1 pol
+val paddX : 'a1 -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
-val psubX :
- 'a1 -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol -> 'a1
- pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
+val psubX : 'a1 -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
-val padd :
- 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol ->
- 'a1 pol
+val padd : 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
-val psub :
- 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1
- -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
+val psub : 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
-val pmulC_aux :
- 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 -> 'a1
- pol
+val pmulC_aux : 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 -> 'a1 pol
-val pmulC :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1
- -> 'a1 pol
+val pmulC : 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 -> 'a1 pol
-val pmulI :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 pol ->
- 'a1 pol -> 'a1 pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
+val pmulI : 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol -> 'a1 pol) -> 'a1 pol -> positive -> 'a1 pol -> 'a1 pol
-val pmul :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
+val pmul : 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
-val psquare :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> 'a1 pol -> 'a1 pol
+val psquare : 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol
type 'c pExpr =
| PEc of 'c
@@ -220,16 +197,12 @@ type 'c pExpr =
val mk_X : 'a1 -> 'a1 -> positive -> 'a1 pol
val ppow_pos :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> ('a1 pol -> 'a1 pol) -> 'a1 pol -> 'a1 pol -> positive -> 'a1 pol
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol) -> 'a1 pol -> 'a1 pol -> positive -> 'a1 pol
-val ppow_N :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> ('a1 pol -> 'a1 pol) -> 'a1 pol -> n -> 'a1 pol
+val ppow_N : 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 pol -> 'a1 pol) -> 'a1 pol -> n -> 'a1 pol
val norm_aux :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pExpr -> 'a1 pol
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pExpr -> 'a1 pol
type 'a bFormula =
| TT
@@ -251,32 +224,22 @@ val tt : 'a1 cnf
val ff : 'a1 cnf
-val add_term :
- ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 -> 'a1 clause -> 'a1
- clause option
+val add_term : ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 -> 'a1 clause -> 'a1 clause option
-val or_clause :
- ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 clause -> 'a1 clause ->
- 'a1 clause option
+val or_clause : ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 clause -> 'a1 clause -> 'a1 clause option
-val or_clause_cnf :
- ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 clause -> 'a1 cnf -> 'a1
- cnf
+val or_clause_cnf : ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 clause -> 'a1 cnf -> 'a1 cnf
-val or_cnf :
- ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 cnf -> 'a1 cnf -> 'a1 cnf
+val or_cnf : ('a1 -> bool) -> ('a1 -> 'a1 -> 'a1 option) -> 'a1 cnf -> 'a1 cnf -> 'a1 cnf
val and_cnf : 'a1 cnf -> 'a1 cnf -> 'a1 cnf
-val xcnf :
- ('a2 -> bool) -> ('a2 -> 'a2 -> 'a2 option) -> ('a1 -> 'a2 cnf) -> ('a1 ->
- 'a2 cnf) -> bool -> 'a1 bFormula -> 'a2 cnf
+val xcnf : ('a2 -> bool) -> ('a2 -> 'a2 -> 'a2 option) -> ('a1 -> 'a2 cnf) -> ('a1 -> 'a2 cnf) -> bool -> 'a1 bFormula -> 'a2 cnf
val cnf_checker : ('a1 list -> 'a2 -> bool) -> 'a1 cnf -> 'a2 list -> bool
val tauto_checker :
- ('a2 -> bool) -> ('a2 -> 'a2 -> 'a2 option) -> ('a1 -> 'a2 cnf) -> ('a1 ->
- 'a2 cnf) -> ('a2 list -> 'a3 -> bool) -> 'a1 bFormula -> 'a3 list -> bool
+ ('a2 -> bool) -> ('a2 -> 'a2 -> 'a2 option) -> ('a1 -> 'a2 cnf) -> ('a1 -> 'a2 cnf) -> ('a2 list -> 'a3 -> bool) -> 'a1 bFormula -> 'a3 list -> bool
val cneqb : ('a1 -> 'a1 -> bool) -> 'a1 -> 'a1 -> bool
@@ -307,32 +270,24 @@ type 'c psatz =
val map_option : ('a1 -> 'a2 option) -> 'a1 option -> 'a2 option
-val map_option2 :
- ('a1 -> 'a2 -> 'a3 option) -> 'a1 option -> 'a2 option -> 'a3 option
+val map_option2 : ('a1 -> 'a2 -> 'a3 option) -> 'a1 option -> 'a2 option -> 'a3 option
val pexpr_times_nformula :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> 'a1 polC -> 'a1 nFormula -> 'a1 nFormula option
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 polC -> 'a1 nFormula -> 'a1 nFormula option
val nformula_times_nformula :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> 'a1 nFormula -> 'a1 nFormula -> 'a1 nFormula option
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula -> 'a1 nFormula -> 'a1 nFormula option
-val nformula_plus_nformula :
- 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula -> 'a1
- nFormula -> 'a1 nFormula option
+val nformula_plus_nformula : 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula -> 'a1 nFormula -> 'a1 nFormula option
val eval_Psatz :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula list -> 'a1 psatz -> 'a1
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula list -> 'a1 psatz -> 'a1
nFormula option
-val check_inconsistent :
- 'a1 -> ('a1 -> 'a1 -> bool) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula -> bool
+val check_inconsistent : 'a1 -> ('a1 -> 'a1 -> bool) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula -> bool
val check_normalised_formulas :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- bool) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula list -> 'a1 psatz -> bool
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> ('a1 -> 'a1 -> bool) -> 'a1 nFormula list -> 'a1 psatz -> bool
type op2 =
| OpEq
@@ -345,36 +300,27 @@ type op2 =
type 't formula = { flhs : 't pExpr; fop : op2; frhs : 't pExpr }
val norm :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pExpr -> 'a1 pol
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pExpr -> 'a1 pol
-val psub0 :
- 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1
- -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
+val psub0 : 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
-val padd0 :
- 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol ->
- 'a1 pol
+val padd0 : 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 pol -> 'a1 pol -> 'a1 pol
val xnormalise :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1 nFormula
- list
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1
+ nFormula list
val cnf_normalise :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1 nFormula
- cnf
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1
+ nFormula cnf
val xnegate :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1 nFormula
- list
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1
+ nFormula list
val cnf_negate :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 ->
- 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1 nFormula
- cnf
+ 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 formula -> 'a1
+ nFormula cnf
val xdenorm : positive -> 'a1 pol -> 'a1 pExpr
@@ -384,9 +330,7 @@ val map_PExpr : ('a2 -> 'a1) -> 'a2 pExpr -> 'a1 pExpr
val map_Formula : ('a2 -> 'a1) -> 'a2 formula -> 'a1 formula
-val simpl_cone :
- 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 psatz ->
- 'a1 psatz
+val simpl_cone : 'a1 -> 'a1 -> ('a1 -> 'a1 -> 'a1) -> ('a1 -> 'a1 -> bool) -> 'a1 psatz -> 'a1 psatz
type q = { qnum : z; qden : positive }
@@ -431,7 +375,7 @@ val psub1 : z pol -> z pol -> z pol
val padd1 : z pol -> z pol -> z pol
-val norm0 : z pExpr -> z pol
+val normZ : z pExpr -> z pol
val xnormalise0 : z formula -> z nFormula list
@@ -487,6 +431,8 @@ val qunsat : q nFormula -> bool
val qdeduce : q nFormula -> q nFormula -> q nFormula option
+val normQ : q pExpr -> q pol
+
val qTautoChecker : q formula bFormula -> qWitness list -> bool
type rcst =
diff --git a/plugins/micromega/micromega_plugin.mlpack b/plugins/micromega/micromega_plugin.mlpack
index ed253da3fd..2baf6608a4 100644
--- a/plugins/micromega/micromega_plugin.mlpack
+++ b/plugins/micromega/micromega_plugin.mlpack
@@ -1,8 +1,11 @@
-Sos_types
Mutils
+Itv
+Vect
+Sos_types
Micromega
Polynomial
Mfourier
+Simplex
Certificate
Persistent_cache
Coq_micromega
diff --git a/plugins/micromega/mutils.ml b/plugins/micromega/mutils.ml
index 9d03560b71..40aeef3959 100644
--- a/plugins/micromega/mutils.ml
+++ b/plugins/micromega/mutils.ml
@@ -19,11 +19,31 @@
(* *)
(************************************************************************)
-let rec pp_list f o l =
+
+module ISet = Set.Make(Int)
+
+module IMap =
+ struct
+ include Map.Make(Int)
+
+ let from k m =
+ let (_,_,r) = split (k-1) m in
+ r
+ end
+
+(*let output_int o i = output_string o (string_of_int i)*)
+
+let iset_pp o s =
+ Printf.fprintf o "{ %a }"
+ (fun o s -> ISet.iter (fun i -> Printf.fprintf o "%i " i) s) s
+
+let rec pp_list s f o l =
match l with
| [] -> ()
- | e::l -> f o e ; output_string o ";" ; pp_list f o l
+ | [e] -> f o e
+ | e::l -> f o e ; output_string o s ; pp_list s f o l
+let output_bigint o bi = output_string o (Big_int.string_of_big_int bi)
let finally f rst =
try
@@ -79,6 +99,12 @@ let extract pred l =
| _ -> (fd, e::sys)
) (None,[]) l
+let extract_all pred l =
+ List.fold_left (fun (s1,s2) e ->
+ match pred e with
+ | None -> s1,e::s2
+ | Some v -> (v,e)::s1 , s2) ([],[]) l
+
open Num
open Big_int
@@ -117,7 +143,22 @@ let rats_to_ints l =
List.map (fun x -> (div_big_int (mult_big_int (numerator x) c)
(denominator x))) l
-(* assoc_pos j [a0...an] = [j,a0....an,j+n],j+n+1 *)
+let iterate_until_stable f x =
+ let rec iter x =
+ match f x with
+ | None -> x
+ | Some x' -> iter x' in
+ iter x
+
+let rec app_funs l x =
+ match l with
+ | [] -> None
+ | f::fl ->
+ match f x with
+ | None -> app_funs fl x
+ | Some x' -> Some x'
+
+
(**
* MODULE: Coq to Caml data-structure mappings
*)
diff --git a/plugins/micromega/mutils.mli b/plugins/micromega/mutils.mli
index 094429ea18..35ca1e5516 100644
--- a/plugins/micromega/mutils.mli
+++ b/plugins/micromega/mutils.mli
@@ -8,6 +8,22 @@
(* * (see LICENSE file for the text of the license) *)
(************************************************************************)
+
+module ISet : Set.S with type elt = int
+
+module IMap :
+sig
+ include Map.S with type key = int
+
+ (** [from k m] returns the submap of [m] with keys greater or equal k *)
+ val from : key -> 'elt t -> 'elt t
+
+end
+
+val iset_pp : out_channel -> ISet.t -> unit
+
+val output_bigint : out_channel -> Big_int.big_int -> unit
+
val numerator : Num.num -> Big_int.big_int
val denominator : Num.num -> Big_int.big_int
@@ -30,7 +46,7 @@ end
module TagSet : CSig.SetS with type elt = Tag.t
-val pp_list : (out_channel -> 'a -> unit) -> out_channel -> 'a list -> unit
+val pp_list : string -> (out_channel -> 'a -> unit) -> out_channel -> 'a list -> unit
module CamlToCoq : sig
@@ -56,6 +72,7 @@ module CoqToCaml : sig
end
+val ppcm : Big_int.big_int -> Big_int.big_int -> Big_int.big_int
val rats_to_ints : Num.num list -> Big_int.big_int list
val all_pairs : ('a -> 'a -> 'b) -> 'a list -> 'b list
@@ -67,4 +84,10 @@ val gcd_list : Num.num list -> Big_int.big_int
val extract : ('a -> 'b option) -> 'a list -> ('b * 'a) option * 'a list
+val extract_all : ('a -> 'b option) -> 'a list -> ('b * 'a) list * 'a list
+
+val iterate_until_stable : ('a -> 'a option) -> 'a -> 'a
+
+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 1d18a26f33..5f31b6f145 100644
--- a/plugins/micromega/polynomial.ml
+++ b/plugins/micromega/polynomial.ml
@@ -10,7 +10,7 @@
(* *)
(* Micromega: A reflexive tactic using the Positivstellensatz *)
(* *)
-(* Frédéric Besson (Irisa/Inria) 2006-20011 *)
+(* Frédéric Besson (Irisa/Inria) 2006-20018 *)
(* *)
(************************************************************************)
@@ -18,6 +18,10 @@ open Num
module Utils = Mutils
open Utils
+module Mc = Micromega
+
+let max_nb_cstr = ref max_int
+
type var = int
let debug = false
@@ -25,652 +29,882 @@ let debug = false
let (<+>) = add_num
let (<*>) = mult_num
-
module Monomial :
sig
- type t
- val const : t
- val is_const : t -> bool
- val var : var -> t
- val is_var : t -> bool
- 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
+ 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
end
- =
-struct
- (* A monomial is represented by a multiset of variables *)
- module Map = Map.Make(Int)
- open Map
-
- type t = int Map.t
-
- let pp o m = Map.iter
- (fun k v ->
- if v = 1 then Printf.fprintf o "x%i." k
- else Printf.fprintf o "x%i^%i." k v) 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 =
- try
- not (Map.fold (fun _ i fk ->
- if fk = true (* first key *)
- then
- if i = 1 then false
- else raise Not_found
- else raise Not_found) m true)
- with Not_found -> false
-
- 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' m
- 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 fold = fold
+ = struct
+ (* A monomial is represented by a multiset of variables *)
+ module Map = Map.Make(Int)
+ open Map
+
+ type t = int Map.t
+
+ 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 eq_num (Int 0) i then ()
+ else Printf.fprintf o "%s" (string_of_num i)
+ else
+ match i with
+ | Int 1 -> Monomial.pp o m
+ | Int -1 -> Printf.fprintf o "-%a" Monomial.pp m
+ | Int 0 -> ()
+ | _ -> Printf.fprintf o "%s*%a" (string_of_num i) Monomial.pp m
+
+
+
module Poly :
- (* A polynomial is a map of monomials *)
- (*
- This is probably a naive implementation
+(* 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 product : t -> t -> t
- val addition : t -> t -> t
- val uminus : t -> t
- val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a
- val is_linear : t -> bool
-end =
-struct
- (*normalisation bug : 0*x ... *)
- module P = Map.Make(Monomial)
- open P
-
- type t = num P.t
-
- (* 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 ->
+ type t
+ val pp : out_channel -> t -> unit
+ val get : Monomial.t -> t -> num
+ val variable : var -> t
+ val add : Monomial.t -> num -> t -> t
+ val constant : num -> 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 is_linear : t -> bool
+ val variables : t -> ISet.t
+ val factorise : var -> t -> t * t
+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 mn i -> Printf.fprintf o "%a + " pp_mon (mn, i)) 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 ->
if sign_num v = 0 then p
else
let vl = (get mn p) <+> v in
- if sign_num vl = 0 then
- remove mn p
- else add mn vl p
+ if sign_num vl = 0 then
+ remove mn p
+ else add mn vl p
- (** Design choice: empty is not a polynomial
- I do not remember why ....
- **)
+ (** 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 ->
- if sign_num v = 0
+ (* The product by a monomial *)
+ let mult : Monomial.t -> num -> t -> t =
+ fun mn v p ->
+ if sign_num v = 0
then constant (Int 0)
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 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 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 uminus : t -> t =
+ fun p -> map (fun v -> minus_num v) p
- let is_linear p = P.fold (fun m _ acc -> acc && (Monomial.is_const m || Monomial.is_var m)) p true
+ let fold = P.fold
-(* let is_linear p =
- let res = is_linear p in
- Printf.printf "is_linear %a = %b\n" pp p res ; res
-*)
-end
+ let is_linear p = P.fold (fun m _ acc -> acc && (Monomial.is_const m || Monomial.is_var m)) p true
+ let variables p = P.fold (fun m _ acc -> ISet.union (Monomial.variables m) acc) p ISet.empty
+
+ 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 (Int 0) , constant (Int 0))
+
+end
-module Vect =
- struct
- (** [t] is the type of vectors.
- A vector [(x1,v1) ; ... ; (xn,vn)] is such that:
- - variables indexes are ordered (x1 <c ... < xn
- - values are all non-zero
- *)
- type var = int
- type t = (var * num) list
-
-(** [equal v1 v2 = true] if the vectors are syntactically equal. *)
-
- let rec equal v1 v2 =
- match v1 , v2 with
- | [] , [] -> true
- | [] , _ -> false
- | _::_ , [] -> false
- | (i1,n1)::v1 , (i2,n2)::v2 ->
- (Int.equal 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 )
-
-
- 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
-
- xfrom_list 0 l
-
- let zero_num = Int 0
-
-
- 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 Int.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 Int.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 mul z t =
- match z with
- | Int 0 -> []
- | Int 1 -> t
- | _ -> List.map (fun (i,n) -> (i, mult_num z n)) t
-
-
- let rec add v1 v2 =
- match v1 , v2 with
- | (x1,n1)::v1' , (x2,n2)::v2' ->
- if x1 = x2
- then
- let n' = n1 +/ n2 in
- if n' =/ Int 0 then add v1' v2'
- else
- let res = add v1' v2' in
- (x1,n') ::res
- else if x1 < x2
- then let res = add v1' v2 in
- (x1, n1)::res
- else let res = add v1 v2' in
- (x2, n2)::res
- | [] , [] -> []
- | [] , _ -> v2
- | _ , [] -> v1
-
-
-
-
- let compare : t -> t -> int = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical
- [
- (fun () -> Int.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 Int.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
type vector = Vect.t
-type cstr_compat = {coeffs : vector ; op : op ; cst : num}
-and op = |Eq | Ge
+type cstr = {coeffs : vector ; op : op ; cst : num}
+and op = |Eq | Ge | Gt
-let string_of_op = function Eq -> "=" | Ge -> ">="
+exception Strict
-let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} =
- Printf.fprintf o "%a %s %s" Vect.pp_vect coeffs (string_of_op op) (string_of_num cst)
+let is_strict c = Pervasives.(=) c.op Gt
-let opMult o1 o2 =
- match o1, o2 with
- | Eq , Eq -> Eq
- | Eq , Ge | Ge , Eq -> Ge
- | Ge , Ge -> Ge
-
-open Big_int
-
-type prf_rule =
- | Hyp of int
- | Def of int
- | Cst of big_int
- | Zero
- | Square of (Vect.t * num)
- | MulC of (Vect.t * num) * prf_rule
- | Gcd of big_int * 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
- | Enum of int * prf_rule * Vect.t * prf_rule * proof list
-
-
-let rec output_prf_rule o = function
- | Hyp i -> Printf.fprintf o "Hyp %i" i
- | Def i -> Printf.fprintf o "Def %i" i
- | Cst c -> Printf.fprintf o "Cst %s" (string_of_big_int c)
- | Zero -> Printf.fprintf o "Zero"
- | Square _ -> Printf.fprintf o "( )^2"
- | MulC(p,pr) -> Printf.fprintf o "P * %a" 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 (string_of_big_int c)
-
-let rec output_proof o = function
- | Done -> Printf.fprintf o "."
- | Step(i,p,pf) -> Printf.fprintf o "%i:= %a ; %a" i output_prf_rule p output_proof pf
- | Enum(i,p1,v,p2,pl) -> Printf.fprintf o "%i{%a<=%a<=%a}%a" i
- output_prf_rule p1 Vect.pp_vect v output_prf_rule p2
- (pp_list output_proof) pl
-
-let rec pr_rule_max_id = function
- | Hyp i | Def i -> i
- | Cst _ | Zero | Square _ -> -1
- | MulC(_,p) | CutPrf p | Gcd(_,p) -> pr_rule_max_id p
- | MulPrf(p1,p2)| AddPrf(p1,p2) -> max (pr_rule_max_id p1) (pr_rule_max_id p2)
-
-let rec proof_max_id = function
- | Done -> -1
- | Step(i,pr,prf) -> max i (max (pr_rule_max_id pr) (proof_max_id prf))
- | Enum(i,p1,_,p2,l) ->
- let m = max (pr_rule_max_id p1) (pr_rule_max_id p2) in
- List.fold_left (fun i prf -> max i (proof_max_id prf)) (max i m) l
-
-let rec pr_rule_def_cut id = function
- | MulC(p,prf) ->
- let (bds,id',prf') = pr_rule_def_cut id prf in
- (bds, id', MulC(p,prf'))
- | MulPrf(p1,p2) ->
- let (bds1,id,p1) = pr_rule_def_cut id p1 in
- let (bds2,id,p2) = pr_rule_def_cut id p2 in
- (bds2@bds1,id,MulPrf(p1,p2))
- | AddPrf(p1,p2) ->
- let (bds1,id,p1) = pr_rule_def_cut id p1 in
- let (bds2,id,p2) = pr_rule_def_cut id p2 in
- (bds2@bds1,id,AddPrf(p1,p2))
- | CutPrf p ->
- let (bds,id,p) = pr_rule_def_cut id p in
- ((id,p)::bds,id+1,Def id)
- | Gcd(c,p) ->
- let (bds,id,p) = pr_rule_def_cut id p in
- ((id,p)::bds,id+1,Def id)
- | Square _|Cst _|Def _|Hyp _|Zero as x -> ([],id,x)
-
-
-(* Do not define top-level cuts *)
-let pr_rule_def_cut id = function
- | CutPrf p ->
- let (bds,ids,p') = pr_rule_def_cut id p in
- bds,ids, CutPrf p'
- | p -> pr_rule_def_cut id p
-
-
-let rec implicit_cut p =
- match p with
- | CutPrf p -> implicit_cut p
- | _ -> p
-
-
-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,id,p' = pr_rule_def_cut id 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)
- | 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 eval_op = function
+ | Eq -> (=/)
+ | Ge -> (>=/)
+ | Gt -> (>/)
- let bds1,id,p1' = pr_rule_def_cut id (implicit_cut p1) in
- let bds2,id,p2' = pr_rule_def_cut 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 string_of_op = function Eq -> "=" | Ge -> ">=" | Gt -> ">"
-let normalise_proof id prf =
- 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 output_cstr o {coeffs = coeffs ; op = op ; cst = cst} =
+ Printf.fprintf o "%a %s %s" Vect.pp coeffs (string_of_op op) (string_of_num cst)
+let opMult o1 o2 =
+ match o1, o2 with
+ | Eq , _ | _ , Eq -> Eq
+ | Ge , _ | _ , Ge -> Ge
+ | Gt , Gt -> Gt
-let add_proof x y =
- match x, y with
- | Zero , p | p , Zero -> p
- | _ -> AddPrf(x,y)
+let opAdd o1 o2 =
+ match o1, o2 with
+ | Eq , x | x , Eq -> x
+ | Gt , x | x , Gt -> Gt
+ | Ge , Ge -> Ge
-let mul_proof c p =
- match sign_big_int c with
- | 0 -> Zero (* This is likely to be a bug *)
- | -1 -> MulC(([],Big_int c),p) (* [p] should represent an equality *)
- | 1 ->
- if eq_big_int c unit_big_int
- then p
- else MulPrf(Cst c,p)
- | _ -> assert false
-let mul_proof_ext (p,c) prf =
- match p with
- | [] -> mul_proof (numerator c) prf
- | _ -> MulC((p,c),prf)
-
+module LinPoly = struct
+ (** A linear polynomial a0 + a1.x1 + ... + an.xn
+ By convention, the constant a0 is the coefficient of the variable 0.
+ *)
-module LinPoly =
-struct
- type t = Vect.t * num
+ type t = Vect.t
- module MonT =
- struct
+ 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 clear () =
+ let clear () =
index_of_monomial := MonoMap.empty;
- monomial_of_index := IntMap.empty ;
+ monomial_of_index := IntMap.empty ;
fresh := 0
- let register m =
+ let register m =
try
- MonoMap.find m !index_of_monomial
- with Not_found ->
- begin
- 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
- end
+ MonoMap.find m !index_of_monomial
+ with Not_found ->
+ begin
+ 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
+ end
let retrieve i = IntMap.find i !monomial_of_index
+ let _ = register Monomial.const
- end
+ end
- let normalise (v,c) =
- (List.sort (fun x y -> Int.compare (fst x) (fst y)) v , c)
+ let var v = Vect.set (MonT.register (Monomial.var v)) (Int 1) Vect.null
+ let of_monomial m =
+ let v = MonT.register m in
+ Vect.set v (Int 1) Vect.null
- let output_mon o (x,v) =
- Printf.fprintf o "%s.%a +" (string_of_num v) Monomial.pp (MonT.retrieve x)
+ 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 (Int 0)) v
+ let coq_poly_of_linpol cst p =
- let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} =
- Printf.fprintf o "%a %s %s" (pp_list output_mon) coeffs (string_of_op op) (string_of_num cst)
+ 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 (Int 1))) 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 (Int 0))) 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 sign_num 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 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_linear p l =
+
+ Vect.find (fun x v ->
+ if p v
+ then
+ let x' = MonT.retrieve x in
+ match Monomial.get_var x' with
+ | None -> None
+ | Some x -> if is_linear_for x l
+ then Some x
+ else None
+ else None) l
+
+
+ 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 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 variables p = Vect.fold
+ (fun acc v _ ->
+ ISet.union (Monomial.variables (MonT.retrieve v)) acc) ISet.empty p
+
+
+ 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 linpol_of_pol p =
- let (v,c) =
- Poly.fold
- (fun mon num (vct,cst) ->
- if Monomial.is_const mon then (vct,num)
- else
- let vr = MonT.register mon in
- ((vr,num)::vct,cst)) p ([], Int 0) in
- normalise (v,c)
- let mult v m (vect,c) =
- if Monomial.is_const m
- then
- (Vect.mul v vect, v <*> c)
- else
- if sign_num v <> 0
- then
- let hd =
- if sign_num c <> 0
- then [MonT.register m,v <*> c]
- else [] in
-
- let vect = hd @ (List.map (fun (x,n) ->
- let x = MonT.retrieve x in
- let x_m = MonT.register (Monomial.prod m x) in
- (x_m, v <*> n)) vect ) in
- normalise (vect , Int 0)
- else ([],Int 0)
- let mult v m (vect,c) =
- let (vect',c') = mult v m (vect,c) in
- if debug then
- Printf.printf "mult %s %a (%a,%s) -> (%a,%s)\n" (string_of_num v) Monomial.pp m
- (pp_list output_mon) vect (string_of_num c)
- (pp_list output_mon) vect' (string_of_num c') ;
- (vect',c')
+ 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
- let make_lin_pol v mon =
- if Monomial.is_const mon
- then [] , v
- else [MonT.register mon, v],Int 0
-
+end
+
+let output_nlin_cstr o {coeffs = coeffs ; op = op ; cst = cst} =
+ let p = LinPoly.pol_of_linpol coeffs in
+
+ Printf.fprintf o "%a %s %s" Poly.pp p (string_of_op op) (string_of_num cst)
+
+
+module ProofFormat = struct
+ open Big_int
+
+ 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
+
+ type proof =
+ | Done
+ | Step of int * prf_rule * proof
+ | Enum of int * prf_rule * Vect.t * prf_rule * proof list
+
+
+ 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" (string_of_num 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 (string_of_big_int c)
+
+ let rec output_proof o = function
+ | Done -> Printf.fprintf o "."
+ | Step(i,p,pf) -> Printf.fprintf o "%i:= %a ; %a" i output_prf_rule p output_proof pf
+ | 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
+
+ let rec pr_rule_max_id = function
+ | Annot(_,p) -> pr_rule_max_id p
+ | Hyp i | Def i -> i
+ | Cst _ | Zero | Square _ -> -1
+ | MulC(_,p) | CutPrf p | Gcd(_,p) -> pr_rule_max_id p
+ | MulPrf(p1,p2)| AddPrf(p1,p2) -> max (pr_rule_max_id p1) (pr_rule_max_id p2)
+
+ let rec proof_max_id = function
+ | Done -> -1
+ | Step(i,pr,prf) -> max i (max (pr_rule_max_id pr) (proof_max_id prf))
+ | Enum(i,p1,_,p2,l) ->
+ let m = max (pr_rule_max_id p1) (pr_rule_max_id p2) in
+ List.fold_left (fun i prf -> max i (proof_max_id prf)) (max i m) l
+
+
+ let rec pr_rule_def_cut id = function
+ | Annot(_,p) -> pr_rule_def_cut id p
+ | MulC(p,prf) ->
+ let (bds,id',prf') = pr_rule_def_cut id prf in
+ (bds, id', MulC(p,prf'))
+ | MulPrf(p1,p2) ->
+ let (bds1,id,p1) = pr_rule_def_cut id p1 in
+ let (bds2,id,p2) = pr_rule_def_cut id p2 in
+ (bds2@bds1,id,MulPrf(p1,p2))
+ | AddPrf(p1,p2) ->
+ let (bds1,id,p1) = pr_rule_def_cut id p1 in
+ let (bds2,id,p2) = pr_rule_def_cut id p2 in
+ (bds2@bds1,id,AddPrf(p1,p2))
+ | CutPrf p ->
+ let (bds,id,p) = pr_rule_def_cut id p in
+ ((id,p)::bds,id+1,Def id)
+ | Gcd(c,p) ->
+ let (bds,id,p) = pr_rule_def_cut id p in
+ ((id,p)::bds,id+1,Def id)
+ | Square _|Cst _|Def _|Hyp _|Zero as x -> ([],id,x)
+
+
+ (* Do not define top-level cuts *)
+ let pr_rule_def_cut id = function
+ | CutPrf p ->
+ let (bds,ids,p') = pr_rule_def_cut id p in
+ bds,ids, CutPrf p'
+ | p -> pr_rule_def_cut id p
+
+
+ let rec implicit_cut p =
+ match p with
+ | CutPrf p -> implicit_cut p
+ | _ -> p
+
+
+ let rec pr_rule_collect_hyps pr =
+ match pr with
+ | Annot(_,pr) -> pr_rule_collect_hyps pr
+ | Hyp i | Def i -> ISet.add i ISet.empty
+ | Cst _ | Zero | Square _ -> ISet.empty
+ | MulC(_,pr) | Gcd(_,pr)| CutPrf pr -> pr_rule_collect_hyps pr
+ | MulPrf(p1,p2) | AddPrf(p1,p2) -> ISet.union (pr_rule_collect_hyps p1) (pr_rule_collect_hyps p2)
+
+ let simplify_proof p =
+ let rec simplify_proof p =
+ match p with
+ | Done -> (Done, ISet.empty)
+ | Step(i,pr,Done) -> (p, ISet.add i (pr_rule_collect_hyps 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_hyps pr) hyps))
+ | 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_hyps p1) (pr_rule_collect_hyps p2)) hyps)) in
+ fst (simplify_proof p)
+
+
+ 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,id,p' = pr_rule_def_cut id 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)
+ | 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,id,p1' = pr_rule_def_cut id (implicit_cut p1) in
+ let bds2,id,p2' = pr_rule_def_cut 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 = 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 add_proof x y =
+ match x, y with
+ | Zero , p | p , Zero -> p
+ | _ -> AddPrf(x,y)
+
+
+ let mul_cst_proof c p =
+ match sign_num c with
+ | 0 -> Zero (* This is likely to be a bug *)
+ | -1 -> MulC(LinPoly.constant c,p) (* [p] should represent an equality *)
+ | 1 ->
+ if eq_num (Int 1) c
+ then p
+ else MulPrf(Cst c,p)
+ | _ -> assert false
+ let mul_proof p1 p2 =
+ match p1 , p2 with
+ | Zero , _ | _ , Zero -> Zero
+ | Cst (Int 1) , p | p , Cst (Int 1) -> p
+ | _ , _ -> MulPrf(p1,p2)
- let xpivot_eq (c,prf) x v (c',prf') =
- if debug then Printf.printf "xpivot_eq {%a} %a %s {%a}\n"
- output_cstr c
- Monomial.pp (MonT.retrieve x)
- (string_of_num v) output_cstr c' ;
+ 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
- let {coeffs = coeffs ; op = op ; cst = cst} = c' in
- let m = MonT.retrieve x in
- let apply_pivot (vqn,q,n) (c',prf') =
- (* Morally, we have (Vect.get (q*x^n) c'.coeffs) = vmn with n >=0 *)
+ module Env = struct
- let cc' = abs_num v in
- let cc_num = Int (- (sign_num v)) <*> vqn in
- let cc_mon = Monomial.prod q (Monomial.exp m (n-1)) in
+ let rec string_of_int_list l =
+ match l with
+ | [] -> ""
+ | i::l -> Printf.sprintf "%i,%s" i (string_of_int_list l)
- let (c_coeff,c_cst) = mult cc_num cc_mon (c.coeffs, minus_num c.cst) in
-
- let c' = {coeffs = Vect.add (Vect.mul cc' c'.coeffs) c_coeff ; op = op ; cst = (minus_num c_cst) <+> (cc' <*> c'.cst)} in
- let prf' = add_proof
- (mul_proof_ext (make_lin_pol cc_num cc_mon) prf)
- (mul_proof (numerator cc') prf') in
- if debug then Printf.printf "apply_pivot -> {%a}\n" output_cstr c' ;
- (c',prf') in
+ let id_of_hyp hyp l =
+ let rec xid_of_hyp i l' =
+ match l' with
+ | [] -> failwith (Printf.sprintf "id_of_hyp %i %s" hyp (string_of_int_list l))
+ | hyp'::l' -> if Pervasives.(=) hyp hyp' then i else xid_of_hyp (i+1) l' in
+ xid_of_hyp 0 l
+ end
+
+ let cmpl_prf_rule norm (cst:num-> 'a) env prf =
+ let rec cmpl =
+ function
+ | Annot(s,p) -> cmpl p
+ | Hyp i | Def i -> Mc.PsatzIn (CamlToCoq.nat (Env.id_of_hyp i 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 (numerator x)) env r
+
+ let rec cmpl_proof env = function
+ | Done -> Mc.DoneProof
+ | Step(i,p,prf) ->
+ begin
+ match p with
+ | CutPrf p' ->
+ Mc.CutProof(cmpl_prf_rule_z env p', cmpl_proof (i::env) prf)
+ | _ -> Mc.RatProof(cmpl_prf_rule_z env p,cmpl_proof (i::env) prf)
+ end
+ | Enum(i,p1,_,p2,l) ->
+ Mc.EnumProof(cmpl_prf_rule_z env p1,cmpl_prf_rule_z env p2,List.map (cmpl_proof (i::env)) l)
+
+
+ let compile_proof env prf =
+ let id = 1 + proof_max_id prf in
+ let _,prf = normalise_proof id prf in
+ cmpl_proof 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 Num.compare_num n (Int 0) 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
+ begin 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"
+ end
+ | Gcd(g,p) -> let (v,op) = eval_prf_rule env p in
+ (Vect.div (Big_int 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 (Int 0))
+ 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 Pervasives.(=) 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
+ | 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"
+
+end
+
+module WithProof = struct
+
+ type t = ((LinPoly.t * op) * ProofFormat.prf_rule)
+
+ 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 cmp (q,n) (q',n') =
- if n < n' then -1
- else if n = n' then Monomial.compare q q'
- else 1 in
+ exception InvalidProof
-
- let find_pivot (c',prf') =
- let (v,q,n) = List.fold_left
- (fun (v,q,n) (x,v') ->
- let x = MonT.retrieve x in
- let (q',n') = Monomial.div x m in
- if cmp (q,n) (q',n') = -1 then (v',q',n') else (v,q,n)) (Int 0, Monomial.const,0) c'.coeffs in
- if n > 0 then Some (v,q,n) else None in
+ let zero = ((Vect.null,Eq), ProofFormat.Zero)
- let rec pivot (q,n) (c',prf') =
- match find_pivot (c',prf') with
- | None -> (c',prf')
- | Some(v,q',n') ->
- if cmp (q',n') (q,n) = -1
- then pivot (q',n') (apply_pivot (v,q',n') (c',prf'))
- else (c',prf') in
- pivot (Monomial.const,max_int) (c',prf')
+ let of_cstr (c,prf) =
+ (Vect.set 0 (Num.minus_num (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 pivot_eq x (c,prf) =
- match Vect.get x c.coeffs with
- | None -> (fun x -> None)
- | Some v -> fun cp' -> Some (xpivot_eq (c,prf) x v cp')
+ 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 mult p ((p1,o1),prf1) =
+ match o1 with
+ | Eq -> ((LinPoly.product p p1,o1), ProofFormat.MulC(p, prf1))
+ | Gt| Ge -> let (n,r) = Vect.decomp_cst p in
+ if Vect.is_null r && n >/ Int 0
+ then ((LinPoly.product p p1, o1), ProofFormat.mul_cst_proof n prf1)
+ else raise InvalidProof
+
+
+ let cutting_plane ((p,o),prf) =
+ let (c,p') = Vect.decomp_cst p in
+ let g = (Vect.gcd p') in
+ if (Big_int.eq_big_int Big_int.unit_big_int g) || c =/ Int 0 ||
+ not (Big_int.eq_big_int (denominator c) Big_int.unit_big_int)
+ then None (* Nothing to do *)
+ else
+ let c1 = c // (Big_int g) in
+ let c1' = Num.floor_num c1 in
+ if c1 =/ c1'
+ then None
+ else
+ match o with
+ | Eq -> Some ((Vect.set 0 (Int (-1)) Vect.null,Eq), ProofFormat.Gcd(g,prf))
+ | Gt -> failwith "cutting_plane ignore strict constraints"
+ | Ge ->
+ (* This is a non-trivial common divisor *)
+ Some ((Vect.set 0 c1' (Vect.div (Big_int g) p),o),ProofFormat.Gcd(g, prf))
+
+
+ let construct_sign p =
+ let (c,p') = Vect.decomp_cst p in
+ if Vect.is_null p'
+ then
+ Some (begin match sign_num c with
+ | 0 -> (true, Eq, ProofFormat.Zero)
+ | 1 -> (true,Gt, ProofFormat.Cst c)
+ | _ (*-1*) -> (false,Gt, ProofFormat.Cst (minus_num c))
+ end)
+ else None
+
+
+ let get_sign l p =
+ match construct_sign p with
+ | None -> begin
+ 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
+ end
+ | 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) -> begin
+ 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)
+
+ end
+ | 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) -> begin
+ 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
+ end
+ | (Ge|Gt) , Eq -> failwith "pivot: equality as second argument"
end
+
+
+(* Local Variables: *)
+(* coding: utf-8 *)
+(* End: *)
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
diff --git a/plugins/micromega/simplex.ml b/plugins/micromega/simplex.ml
new file mode 100644
index 0000000000..8d8c6ea90b
--- /dev/null
+++ b/plugins/micromega/simplex.ml
@@ -0,0 +1,621 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+
+(** A naive simplex *)
+open Polynomial
+open Num
+open Util
+open Mutils
+
+let debug = false
+
+type iset = unit IMap.t
+
+type tableau = Vect.t IMap.t (** Mapping basic variables to their equation.
+ All variables >= than a threshold rst are restricted.*)
+module Restricted =
+ struct
+ type t =
+ {
+ base : int; (** All variables above [base] are restricted *)
+ exc : int option (** Except [exc] which is currently optimised *)
+ }
+
+ let pp o {base;exc} =
+ Printf.fprintf o ">= %a " LinPoly.pp_var base;
+ match exc with
+ | None ->Printf.fprintf o "-"
+ | Some x ->Printf.fprintf o "-%a" LinPoly.pp_var base
+
+ let is_exception (x:var) (r:t) =
+ match r.exc with
+ | None -> false
+ | Some x' -> x = x'
+
+ let restrict x rst =
+ if is_exception x rst
+ then
+ {base = rst.base;exc= None}
+ else failwith (Printf.sprintf "Cannot restrict %i" x)
+
+
+ let is_restricted x r0 =
+ x >= r0.base && not (is_exception x r0)
+
+ let make x = {base = x ; exc = None}
+
+ let set_exc x rst = {base = rst.base ; exc = Some x}
+
+ let fold rst f m acc =
+ IMap.fold (fun k v acc ->
+ if is_exception k rst then acc
+ else f k v acc) (IMap.from rst.base m) acc
+
+ end
+
+
+
+let pp_row o v = LinPoly.pp o v
+
+let output_tableau o t =
+ IMap.iter (fun k v ->
+ Printf.fprintf o "%a = %a\n" LinPoly.pp_var k pp_row v) t
+
+let output_vars o m =
+ IMap.iter (fun k _ -> Printf.fprintf o "%a " LinPoly.pp_var k) m
+
+
+(** A tableau is feasible iff for every basic restricted variable xi,
+ we have ci>=0.
+
+ When all the non-basic variables are set to 0, the value of a basic
+ variable xi is necessarily ci. If xi is restricted, it is feasible
+ if ci>=0.
+ *)
+
+
+let unfeasible (rst:Restricted.t) tbl =
+ Restricted.fold rst (fun k v m ->
+ if Vect.get_cst v >=/ Int 0 then m
+ else IMap.add k () m) tbl IMap.empty
+
+
+let is_feasible rst tb = IMap.is_empty (unfeasible rst tb)
+
+(** Let a1.x1+...+an.xn be a vector of non-basic variables.
+ It is maximised if all the xi are restricted
+ and the ai are negative.
+
+ If xi>= 0 (restricted) and ai is negative,
+ the maximum for ai.xi is obtained for xi = 0
+
+ Otherwise, it is possible to make ai.xi arbitrarily big:
+ - if xi is not restricted, take +/- oo depending on the sign of ai
+ - if ai is positive, take +oo
+ *)
+
+let is_maximised_vect rst v =
+ Vect.for_all (fun xi ai ->
+ if ai >/ Int 0
+ then false
+ else Restricted.is_restricted xi rst) v
+
+
+(** [is_maximised rst v]
+ @return None if the variable is not maximised
+ @return Some v where v is the maximal value
+ *)
+let is_maximised rst v =
+ try
+ let (vl,v) = Vect.decomp_cst v in
+ if is_maximised_vect rst v
+ then Some vl
+ else None
+ with Not_found -> None
+
+(** A variable xi is unbounded if for every
+ equation xj= ...ai.xi ...
+ if ai < 0 then xj is not restricted.
+ As a result, even if we
+ increase the value of xi, it is always
+ possible to adjust the value of xj without
+ violating a restriction.
+ *)
+
+(* let is_unbounded rst tbl vr =
+ IMap.for_all (fun x v -> if Vect.get vr v </ Int 0
+ then not (IMap.mem vr rst)
+ else true
+ ) tbl
+ *)
+
+type result =
+ | Max of num (** Maximum is reached *)
+ | Ubnd of var (** Problem is unbounded *)
+ | Feas (** Problem is feasible *)
+
+type pivot =
+ | Done of result
+ | Pivot of int * int * num
+
+
+
+
+type simplex =
+ | Opt of tableau * result
+
+(** For a row, x = ao.xo+...+ai.xi
+ a valid pivot variable is such that it can improve the value of xi.
+ it is the case, if xi is unrestricted (increase if ai> 0, decrease if ai < 0)
+ xi is restricted but ai > 0
+
+This is the entering variable.
+ *)
+
+let rec find_pivot_column (rst:Restricted.t) (r:Vect.t) =
+ match Vect.choose r with
+ | None -> failwith "find_pivot_column"
+ | Some(xi,ai,r') -> if ai </ Int 0
+ then if Restricted.is_restricted xi rst
+ then find_pivot_column rst r' (* ai.xi cannot be improved *)
+ else (xi, -1) (* r is not restricted, sign of ai does not matter *)
+ else (* ai is positive, xi can be increased *)
+ (xi,1)
+
+(** Finding the variable leaving the basis is more subtle because we need to:
+ - increase the objective function
+ - make sure that the entering variable has a feasible value
+ - but also that after pivoting all the other basic variables are still feasible.
+ This explains why we choose the pivot with the smallest score
+ *)
+
+let min_score s (i1,sc1) =
+ match s with
+ | None -> Some (i1,sc1)
+ | Some(i0,sc0) ->
+ if sc0 </ sc1 then s
+ else if sc1 </ sc0 then Some (i1,sc1)
+ else if i0 < i1 then s else Some(i1,sc1)
+
+let find_pivot_row rst tbl j sgn =
+ Restricted.fold rst
+ (fun i' v res ->
+ let aij = Vect.get j v in
+ if (Int sgn) */ aij </ Int 0
+ then (* This would improve *)
+ let score' = Num.abs_num ((Vect.get_cst v) // aij) in
+ min_score res (i',score')
+ else res) tbl None
+
+let safe_find err x t =
+ try
+ IMap.find x t
+ with Not_found ->
+ if debug
+ then Printf.fprintf stdout "safe_find %s x%i %a\n" err x output_tableau t;
+ failwith err
+
+
+(** [find_pivot vr t] aims at improving the objective function of the basic variable vr *)
+let find_pivot vr (rst:Restricted.t) tbl =
+ (* Get the objective of the basic variable vr *)
+ let v = safe_find "find_pivot" vr tbl in
+ match is_maximised rst v with
+ | Some mx -> Done (Max mx) (* Maximum is reached; we are done *)
+ | None ->
+ (* Extract the vector *)
+ let (_,v) = Vect.decomp_cst v in
+ let (j',sgn) = find_pivot_column rst v in
+ match find_pivot_row rst (IMap.remove vr tbl) j' sgn with
+ | None -> Done (Ubnd j')
+ | Some (i',sc) -> Pivot(i', j', sc)
+
+(** [solve_column c r e]
+ @param c is a non-basic variable
+ @param r is a basic variable
+ @param e is a vector such that r = e
+ and e is of the form ai.c+e'
+ @return the vector (-r + e').-1/ai i.e
+ c = (r - e')/ai
+ *)
+
+let solve_column (c : var) (r : var) (e : Vect.t) : Vect.t =
+ let a = Vect.get c e in
+ if a =/ Int 0
+ then failwith "Cannot solve column"
+ else
+ let a' = (Int (-1) // a) in
+ Vect.mul a' (Vect.set r (Int (-1)) (Vect.set c (Int 0) e))
+
+(** [pivot_row r c e]
+ @param c is such that c = e
+ @param r is a vector r = g.c + r'
+ @return g.e+r' *)
+
+let pivot_row (row: Vect.t) (c : var) (e : Vect.t) : Vect.t =
+ let g = Vect.get c row in
+ if g =/ Int 0
+ then row
+ else Vect.mul_add g e (Int 1) (Vect.set c (Int 0) row)
+
+let pivot_with (m : tableau) (v: var) (p : Vect.t) =
+ IMap.map (fun (r:Vect.t) -> pivot_row r v p) m
+
+let pivot (m : tableau) (r : var) (c : var) =
+ let row = safe_find "pivot" r m in
+ let piv = solve_column c r row in
+ IMap.add c piv (pivot_with (IMap.remove r m) c piv)
+
+
+let adapt_unbounded vr x rst tbl =
+ if Vect.get_cst (IMap.find vr tbl) >=/ Int 0
+ then tbl
+ else pivot tbl vr x
+
+module BaseSet = Set.Make(struct type t = iset let compare = IMap.compare (fun x y -> 0) end)
+
+let get_base tbl = IMap.mapi (fun k _ -> ()) tbl
+
+let simplex opt vr rst tbl =
+ let b = ref BaseSet.empty in
+
+let rec simplex opt vr rst tbl =
+
+ if debug then begin
+ let base = get_base tbl in
+ if BaseSet.mem base !b
+ then Printf.fprintf stdout "Cycling detected\n"
+ else b := BaseSet.add base !b
+ end;
+
+ if debug && not (is_feasible rst tbl)
+ then
+ begin
+ let m = unfeasible rst tbl in
+ Printf.fprintf stdout "Simplex error\n";
+ Printf.fprintf stdout "The current tableau is not feasible\n";
+ Printf.fprintf stdout "Restricted >= %a\n" Restricted.pp rst ;
+ output_tableau stdout tbl;
+ Printf.fprintf stdout "Error for variables %a\n" output_vars m
+ end;
+
+ if not opt && (Vect.get_cst (IMap.find vr tbl) >=/ Int 0)
+ then Opt(tbl,Feas)
+ else
+ match find_pivot vr rst tbl with
+ | Done r ->
+ begin match r with
+ | Max _ -> Opt(tbl, r)
+ | Ubnd x ->
+ let t' = adapt_unbounded vr x rst tbl in
+ Opt(t',r)
+ | Feas -> raise (Invalid_argument "find_pivot")
+ end
+ | Pivot(i,j,s) ->
+ if debug then begin
+ Printf.fprintf stdout "Find pivot for x%i(%s)\n" vr (string_of_num s);
+ Printf.fprintf stdout "Leaving variable x%i\n" i;
+ Printf.fprintf stdout "Entering variable x%i\n" j;
+ end;
+ let m' = pivot tbl i j in
+ simplex opt vr rst m' in
+
+simplex opt vr rst tbl
+
+
+
+type certificate =
+ | Unsat of Vect.t
+ | Sat of tableau * var option
+
+(** [normalise_row t v]
+ @return a row obtained by pivoting the basic variables of the vector v
+ *)
+
+let normalise_row (t : tableau) (v: Vect.t) =
+ Vect.fold (fun acc vr ai -> try
+ let e = IMap.find vr t in
+ Vect.add (Vect.mul ai e) acc
+ with Not_found -> Vect.add (Vect.set vr ai Vect.null) acc)
+ Vect.null v
+
+let normalise_row (t : tableau) (v: Vect.t) =
+ let v' = normalise_row t v in
+ if debug then Printf.fprintf stdout "Normalised Optimising %a\n" LinPoly.pp v';
+ v'
+
+let add_row (nw :var) (t : tableau) (v : Vect.t) : tableau =
+ IMap.add nw (normalise_row t v) t
+
+(** [push_real] performs reasoning over the rationals *)
+let push_real (opt : bool) (nw : var) (v : Vect.t) (rst: Restricted.t) (t : tableau) : certificate =
+ if debug
+ then begin Printf.fprintf stdout "Current Tableau\n%a\n" output_tableau t;
+ Printf.fprintf stdout "Optimising %a=%a\n" LinPoly.pp_var nw LinPoly.pp v
+ end;
+ match simplex opt nw rst (add_row nw t v) with
+ | Opt(t',r) -> (* Look at the optimal *)
+ match r with
+ | Ubnd x->
+ if debug then Printf.printf "The objective is unbounded (variable %a)\n" LinPoly.pp_var x;
+ Sat (t',Some x) (* This is sat and we can extract a value *)
+ | Feas -> Sat (t',None)
+ | Max n ->
+ if debug then begin
+ Printf.printf "The objective is maximised %s\n" (string_of_num n);
+ Printf.printf "%a = %a\n" LinPoly.pp_var nw pp_row (IMap.find nw t')
+ end;
+
+ if n >=/ Int 0
+ then Sat (t',None)
+ else
+ let v' = safe_find "push_real" nw t' in
+ Unsat (Vect.set nw (Int 1) (Vect.set 0 (Int 0) (Vect.mul (Int (-1)) v')))
+
+
+(** One complication is that equalities needs some pre-processing.contents
+ *)
+open Mutils
+open Polynomial
+
+let fresh_var l =
+ 1 +
+ try
+ (ISet.max_elt (List.fold_left (fun acc c -> ISet.union acc (Vect.variables c.coeffs)) ISet.empty l))
+ with Not_found -> 0
+
+
+(*type varmap = (int * bool) IMap.t*)
+
+
+let make_certificate vm l =
+ Vect.normalise (Vect.fold (fun acc x n ->
+ let (x',b) = IMap.find x vm in
+ Vect.set x' (if b then n else Num.minus_num n) acc) Vect.null l)
+
+
+
+
+
+let eliminate_equalities (vr0:var) (l:Polynomial.cstr list) =
+ let rec elim idx vr vm l acc =
+ match l with
+ | [] -> (vr,vm,acc)
+ | c::l -> match c.op with
+ | Ge -> let v = Vect.set 0 (minus_num c.cst) c.coeffs in
+ elim (idx+1) (vr+1) (IMap.add vr (idx,true) vm) l ((vr,v)::acc)
+ | Eq -> let v1 = Vect.set 0 (minus_num c.cst) c.coeffs in
+ let v2 = Vect.mul (Int (-1)) v1 in
+ let vm = IMap.add vr (idx,true) (IMap.add (vr+1) (idx,false) vm) in
+ elim (idx+1) (vr+2) vm l ((vr,v1)::(vr+1,v2)::acc)
+ | Gt -> raise Strict in
+ elim 0 vr0 IMap.empty l []
+
+let find_solution rst tbl =
+ IMap.fold (fun vr v res -> if Restricted.is_restricted vr rst
+ then res
+ else Vect.set vr (Vect.get_cst v) res) tbl Vect.null
+
+let choose_conflict (sol:Vect.t) (l: (var * Vect.t) list) =
+ let esol = Vect.set 0 (Int 1) sol in
+ let is_conflict (x,v) =
+ if Vect.dotproduct esol v >=/ Int 0
+ then None else Some(x,v) in
+ let (c,r) = extract is_conflict l in
+ match c with
+ | Some (c,_) -> Some (c,r)
+ | None -> match l with
+ | [] -> None
+ | e::l -> Some(e,l)
+
+(*let remove_redundant rst t =
+ IMap.fold (fun k v m -> if Restricted.is_restricted k rst && Vect.for_all (fun x _ -> x == 0 || Restricted.is_restricted x rst) v
+ then begin
+ if debug then
+ Printf.printf "%a is redundant\n" LinPoly.pp_var k;
+ IMap.remove k m
+ end
+ else m) t t
+ *)
+
+
+let rec solve opt l (rst:Restricted.t) (t:tableau) =
+ let sol = find_solution rst t in
+ match choose_conflict sol l with
+ | None -> Inl (rst,t,None)
+ | Some((vr,v),l) ->
+ match push_real opt vr v (Restricted.set_exc vr rst) t with
+ | Sat (t',x) ->
+ (* let t' = remove_redundant rst t' in*)
+ begin
+ match l with
+ | [] -> Inl(rst,t', x)
+ | _ -> solve opt l rst t'
+ end
+ | Unsat c -> Inr c
+
+let find_unsat_certificate (l : Polynomial.cstr list ) =
+ let vr = fresh_var l in
+ let (_,vm,l') = eliminate_equalities vr l in
+
+ match solve false l' (Restricted.make vr) IMap.empty with
+ | Inr c -> Some (make_certificate vm c)
+ | Inl _ -> None
+
+
+
+let find_point (l : Polynomial.cstr list) =
+ let vr = fresh_var l in
+ let (_,vm,l') = eliminate_equalities vr l in
+
+ match solve false l' (Restricted.make vr) IMap.empty with
+ | Inl (rst,t,_) -> Some (find_solution rst t)
+ | _ -> None
+
+
+
+let optimise obj l =
+ let vr0 = fresh_var l in
+ let (_,vm,l') = eliminate_equalities (vr0+1) l in
+
+ let bound pos res =
+ match res with
+ | Opt(_,Max n) -> Some (if pos then n else minus_num n)
+ | Opt(_,Ubnd _) -> None
+ | Opt(_,Feas) -> None
+ in
+
+ match solve false l' (Restricted.make vr0) IMap.empty with
+ | Inl (rst,t,_) ->
+ Some (bound false
+ (simplex true vr0 rst (add_row vr0 t (Vect.uminus obj))),
+ bound true
+ (simplex true vr0 rst (add_row vr0 t obj)))
+ | _ -> None
+
+
+
+open Polynomial
+
+let env_of_list l =
+ List.fold_left (fun (i,m) l -> (i+1, IMap.add i l m)) (0,IMap.empty) l
+
+
+open ProofFormat
+
+let make_farkas_certificate (env: WithProof.t IMap.t) vm v =
+ Vect.fold (fun acc x n ->
+ add_proof acc
+ begin
+ try
+ let (x',b) = IMap.find x vm in
+ (mul_cst_proof
+ (if b then n else (Num.minus_num n))
+ (snd (IMap.find x' env)))
+ with Not_found -> (* This is an introduced hypothesis *)
+ (mul_cst_proof n (snd (IMap.find x env)))
+ end) Zero v
+
+let make_farkas_proof (env: WithProof.t IMap.t) vm v =
+ Vect.fold (fun wp x n ->
+ WithProof.addition wp begin
+ try
+ let (x', b) = IMap.find x vm in
+ let n = if b then n else Num.minus_num n in
+ WithProof.mult (Vect.cst n) (IMap.find x' env)
+ with Not_found ->
+ WithProof.mult (Vect.cst n) (IMap.find x env)
+ end) WithProof.zero v
+
+(*
+let incr_cut rmin x =
+ match rmin with
+ | None -> true
+ | Some r -> Int.compare x r = 1
+ *)
+
+let cut env rmin sol vm (rst:Restricted.t) (x,v) =
+(* if not (incr_cut rmin x)
+ then None
+ else *)
+ let (n,r) = Vect.decomp_cst v in
+
+ let nf = Num.floor_num n in
+ if nf =/ n
+ then None (* The solution is integral *)
+ else
+ (* This is potentially a cut *)
+ let cut = Vect.normalise
+ (Vect.fold (fun acc x n ->
+ if Restricted.is_restricted x rst then
+ Vect.set x (n -/ (Num.floor_num n)) acc
+ else acc
+ ) Vect.null r) in
+ if debug then Printf.fprintf stdout "Cut vector for %a : %a\n" LinPoly.pp_var x LinPoly.pp cut ;
+ let cut = make_farkas_proof env vm cut in
+
+ match WithProof.cutting_plane cut with
+ | None -> None
+ | Some (v,prf) ->
+ if debug then begin
+ Printf.printf "This is a cutting plane:\n" ;
+ Printf.printf "%a -> %a\n" WithProof.output cut WithProof.output (v,prf);
+ end;
+ if Pervasives.(=) (snd v) Eq
+ then (* Unsat *) Some (x,(v,prf))
+ else if eval_op Ge (Vect.dotproduct (fst v) (Vect.set 0 (Int 1) sol)) (Int 0)
+ then begin
+ (* Can this happen? *)
+ if debug then Printf.printf "The cut is feasible - drop it\n";
+ None
+ end
+ else Some(x,(v,prf))
+
+let find_cut env u sol vm rst tbl =
+ (* find first *)
+ IMap.fold (fun x v acc ->
+ match acc with
+ | None -> cut env u sol vm rst (x,v)
+ | Some c -> acc) tbl None
+
+(*
+let find_cut env u sol vm rst tbl =
+ IMap.fold (fun x v acc ->
+ match acc with
+ | Some c -> Some c
+ | None -> cut env u sol vm rst (x,v)
+ ) tbl None
+ *)
+
+let integer_solver lp =
+ let (l,_) = List.split lp in
+ let vr0 = fresh_var l in
+ let (vr,vm,l') = eliminate_equalities vr0 l in
+
+ let _,env = env_of_list (List.map WithProof.of_cstr lp) in
+
+ let insert_row vr v rst tbl =
+ match push_real true vr v rst tbl with
+ | Sat (t',x) -> Inl (Restricted.restrict vr rst,t',x)
+ | Unsat c -> Inr c in
+
+ let rec isolve env cr vr res =
+ match res with
+ | Inr c -> Some (Step (vr, make_farkas_certificate env vm (Vect.normalise c),Done))
+ | Inl (rst,tbl,x) ->
+ if debug then begin
+ Printf.fprintf stdout "Looking for a cut\n";
+ Printf.fprintf stdout "Restricted %a ...\n" Restricted.pp rst;
+ Printf.fprintf stdout "Current Tableau\n%a\n" output_tableau tbl;
+ end;
+ let sol = find_solution rst tbl in
+
+ match find_cut env cr (*x*) sol vm rst tbl with
+ | None -> None
+ | Some(cr,((v,op),cut)) ->
+ if Pervasives.(=) op Eq
+ then (* This is a contradiction *)
+ Some(Step(vr,CutPrf cut, Done))
+ else
+ let res = insert_row vr v (Restricted.set_exc vr rst) tbl in
+ let prf = isolve (IMap.add vr ((v,op),Def vr) env) (Some cr) (vr+1) res in
+ match prf with
+ | None -> None
+ | Some p -> Some (Step(vr,CutPrf cut,p)) in
+
+ let res = solve true l' (Restricted.make vr0) IMap.empty in
+ isolve env None vr res
+
+let integer_solver lp =
+ match integer_solver lp with
+ | None -> None
+ | Some prf -> if debug
+ then Printf.fprintf stdout "Proof %a\n" ProofFormat.output_proof prf ;
+ Some prf
diff --git a/plugins/micromega/simplex.mli b/plugins/micromega/simplex.mli
new file mode 100644
index 0000000000..9f87e745eb
--- /dev/null
+++ b/plugins/micromega/simplex.mli
@@ -0,0 +1,18 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+open Polynomial
+
+val optimise : Vect.t -> cstr list -> (Num.num option * Num.num option) option
+
+val find_point : cstr list -> Vect.t option
+
+val find_unsat_certificate : cstr list -> Vect.t option
+
+val integer_solver : (cstr * ProofFormat.prf_rule) list -> ProofFormat.proof option
diff --git a/plugins/micromega/vect.ml b/plugins/micromega/vect.ml
new file mode 100644
index 0000000000..b188ab4278
--- /dev/null
+++ b/plugins/micromega/vect.ml
@@ -0,0 +1,295 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+
+open Num
+open Mutils
+
+(** [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. *)
+
+let rec equal v1 v2 =
+ match v1 , v2 with
+ | [] , [] -> true
+ | [] , _ -> false
+ | _::_ , [] -> false
+ | (i1,n1)::v1 , (i2,n2)::v2 ->
+ (Int.equal 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 )
+
+
+let null = []
+
+let is_null v =
+ match v with
+ | [] | [0,Int 0] -> true
+ | _ -> false
+
+let pp_var_num pp_var o (v,n) =
+ if Int.equal v 0
+ then if eq_num (Int 0) n then ()
+ else Printf.fprintf o "%s" (string_of_num n)
+ else
+ match n with
+ | Int 1 -> pp_var o v
+ | Int -1 -> Printf.fprintf o "-%a" pp_var v
+ | Int 0 -> ()
+ | _ -> Printf.fprintf o "%s*%a" (string_of_num n) pp_var v
+
+
+let rec pp_gen pp_var o v =
+ match v with
+ | [] -> output_string o "0"
+ | [e] -> pp_var_num pp_var o e
+ | e::l -> Printf.fprintf o "%a + %a" (pp_var_num pp_var) e (pp_gen pp_var) l
+
+
+let pp_var o v = Printf.fprintf o "x%i" v
+
+let pp o v = pp_gen pp_var o v
+
+
+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
+
+ xfrom_list 0 l
+
+let zero_num = Int 0
+
+
+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 Int.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 Int.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 cst n = if n =/ Int 0 then [] else [0,n]
+
+
+let mul z t =
+ match z with
+ | Int 0 -> []
+ | Int 1 -> t
+ | _ -> List.map (fun (i,n) -> (i, mult_num z n)) t
+
+let div z t =
+ if z <>/ Int 1
+ then List.map (fun (x,nx) -> (x,nx // z)) t
+ else t
+
+
+let uminus t = List.map (fun (i,n) -> i, minus_num n) t
+
+
+let rec add (ve1:t) (ve2:t) =
+ match ve1 , ve2 with
+ | [] , v | v , [] -> v
+ | (v1,c1)::l1 , (v2,c2)::l2 ->
+ let cmp = Pervasives.compare v1 v2 in
+ if cmp == 0 then
+ let s = add_num c1 c2 in
+ if eq_num (Int 0) s
+ then add l1 l2
+ else (v1,s)::(add l1 l2)
+ else if cmp < 0 then (v1,c1) :: (add l1 ve2)
+ else (v2,c2) :: (add l2 ve1)
+
+
+let rec xmul_add (n1:num) (ve1:t) (n2:num) (ve2:t) =
+ match ve1 , ve2 with
+ | [] , _ -> mul n2 ve2
+ | _ , [] -> mul n1 ve1
+ | (v1,c1)::l1 , (v2,c2)::l2 ->
+ let cmp = Pervasives.compare v1 v2 in
+ if cmp == 0 then
+ let s = ( n1 */ c1) +/ (n2 */ c2) in
+ if eq_num (Int 0) s
+ then xmul_add n1 l1 n2 l2
+ else (v1,s)::(xmul_add n1 l1 n2 l2)
+ else if cmp < 0 then (v1,n1 */ c1) :: (xmul_add n1 l1 n2 ve2)
+ else (v2,n2 */c2) :: (xmul_add n1 ve1 n2 l2)
+
+let mul_add n1 ve1 n2 ve2 =
+ if n1 =/ Int 1 && n2 =/ Int 1
+ then add ve1 ve2
+ else xmul_add n1 ve1 n2 ve2
+
+
+let compare : t -> t -> int = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical
+ [
+ (fun () -> Int.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 Int.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 -> Int 0
+ | Some(vl,_) -> vl
+
+let is_constant v =
+ match v with
+ | [] | [0,_] -> true
+ | _ -> false
+
+
+
+let get_cst vect =
+ match vect with
+ | (0,v)::_ -> v
+ | _ -> Int 0
+
+let choose v =
+ match v with
+ | [] -> None
+ | (vr,vl)::rst -> Some (vr,vl,rst)
+
+
+let rec fresh v =
+ match v with
+ | [] -> 1
+ | [v,_] -> v + 1
+ | _::v -> fresh v
+
+
+let variables v =
+ List.fold_left (fun acc (x,_) -> ISet.add x acc) ISet.empty v
+
+let decomp_cst v =
+ match v with
+ | (0,vl)::v -> vl,v
+ | _ -> Int 0,v
+
+let fold f acc v =
+ List.fold_left (fun acc (v,i) -> f acc v i) acc v
+
+let fold_error f acc v =
+ let rec fold acc v =
+ match v with
+ | [] -> Some acc
+ | (x,i)::v' -> match f acc x i with
+ | None -> None
+ | Some acc' -> fold acc' v' in
+ fold acc v
+
+
+
+let rec find p v =
+ match v with
+ | [] -> None
+ | (v,n)::v' -> match p v n with
+ | None -> find p v'
+ | Some r -> Some r
+
+
+let for_all p l =
+ List.for_all (fun (v,n) -> p v n) l
+
+
+let decr_var i v = List.map (fun (v,n) -> (v-i,n)) v
+let incr_var i v = List.map (fun (v,n) -> (v+i,n)) v
+
+open Big_int
+
+let gcd v =
+ let res = fold (fun c _ n ->
+ assert (Int.equal (compare_big_int (denominator n) unit_big_int) 0);
+ gcd_big_int c (numerator n)) zero_big_int v in
+ if Int.equal (compare_big_int res zero_big_int) 0
+ then unit_big_int else res
+
+let normalise v =
+ let ppcm = fold (fun c _ n -> ppcm c (denominator n)) unit_big_int v in
+ let gcd =
+ let gcd = fold (fun c _ n -> gcd_big_int c (numerator n)) zero_big_int v in
+ if Int.equal (compare_big_int gcd zero_big_int) 0 then unit_big_int else gcd in
+ List.map (fun (x,v) -> (x, v */ (Big_int ppcm) // (Big_int gcd))) v
+
+let rec exists2 p vect1 vect2 =
+ match vect1 , vect2 with
+ | _ , [] | [], _ -> None
+ | (v1,n1)::vect1' , (v2, n2) :: vect2' ->
+ if Int.equal v1 v2
+ then
+ if p n1 n2
+ then Some (v1,n1,n2)
+ else
+ exists2 p vect1' vect2'
+ else
+ if v1 < v2
+ then exists2 p vect1' vect2
+ else exists2 p vect1 vect2'
+
+let dotproduct v1 v2 =
+ let rec dot acc v1 v2 =
+ match v1, v2 with
+ | [] , _ | _ , [] -> acc
+ | (x1,n1)::v1', (x2,n2)::v2' ->
+ if x1 == x2
+ then dot (acc +/ n1 */ n2) v1' v2'
+ else if x1 < x2
+ then dot acc v1' v2
+ else dot acc v1 v2' in
+ dot (Int 0) v1 v2
diff --git a/plugins/micromega/vect.mli b/plugins/micromega/vect.mli
new file mode 100644
index 0000000000..da6b1e8e9b
--- /dev/null
+++ b/plugins/micromega/vect.mli
@@ -0,0 +1,156 @@
+(************************************************************************)
+(* * The Coq Proof Assistant / The Coq Development Team *)
+(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
+(* <O___,, * (see CREDITS file for the list of authors) *)
+(* \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) *)
+(************************************************************************)
+
+open Num
+open Mutils
+
+type var = int (** Variables are simply (positive) integers. *)
+
+type t (** The type of vectors or equivalently linear expressions.
+ The current implementation is using association lists.
+ A list [(0,c),(x1,ai),...,(xn,an)] represents the linear expression
+ c + a1.xn + ... an.xn where ai are rational constants and xi are variables.
+
+ Note that the variable 0 has a special meaning and represent a constant.
+ Moreover, the representation is spare and variables with a zero coefficient
+ are not represented.
+ *)
+
+(** {1 Generic functions} *)
+
+(** [hash] [equal] and [compare] so that Vect.t can be used as
+ keys for Set Map and Hashtbl *)
+
+val hash : t -> int
+val equal : t -> t -> bool
+val compare : t -> t -> int
+
+(** {1 Basic accessors and utility functions} *)
+
+(** [pp_gen pp_var o v] prints the representation of the vector [v] over the channel [o] *)
+val pp_gen : (out_channel -> var -> unit) -> out_channel -> t -> unit
+
+(** [pp o v] prints the representation of the vector [v] over the channel [o] *)
+val pp : out_channel -> t -> unit
+
+(** [variables v] returns the set of variables with non-zero coefficients *)
+val variables : t -> ISet.t
+
+(** [get_cst v] returns c i.e. the coefficient of the variable zero *)
+val get_cst : t -> num
+
+(** [decomp_cst v] returns the pair (c,a1.x1+...+an.xn) *)
+val decomp_cst : t -> num * t
+
+(** [cst c] returns the vector v=c+0.x1+...+0.xn *)
+val cst : num -> t
+
+(** [is_constant v] holds if [v] is a constant vector i.e. v=c+0.x1+...+0.xn
+ *)
+val is_constant : t -> bool
+
+(** [null] is the empty vector i.e. 0+0.x1+...+0.xn *)
+val null : t
+
+(** [is_null v] returns whether [v] is the [null] vector i.e [equal v null] *)
+val is_null : t -> bool
+
+(** [get xi v] returns the coefficient ai of the variable [xi].
+ [get] is also defined for the variable 0 *)
+val get : var -> t -> num
+
+(** [set xi ai' v] returns the vector c+a1.x1+...ai'.xi+...+an.xn
+ i.e. the coefficient of the variable xi is set to ai' *)
+val set : var -> num -> t -> t
+
+(** [update xi f v] returns c+a1.x1+...+f(ai).xi+...+an.xn *)
+val update : var -> (num -> num) -> t -> t
+
+(** [fresh v] return the fresh variable with inded 1+ max (variables v) *)
+val fresh : t -> int
+
+(** [choose v] decomposes a vector [v] depending on whether it is [null] or not.
+ @return None if v is [null]
+ @return Some(x,n,r) where v = r + n.x x is the smallest variable with non-zero coefficient n <> 0.
+ *)
+val choose : t -> (var * num * t) option
+
+(** [from_list l] returns the vector c+a1.x1...an.xn from the list of coefficient [l=c;a1;...;an] *)
+val from_list : num list -> t
+
+(** [to_list v] returns the list of all coefficient of the vector v i.e. [c;a1;...;an]
+ The list representation is (obviously) not sparsed
+ and therefore certain ai may be 0 *)
+val to_list : t -> num list
+
+(** [decr_var i v] decrements the variables of the vector [v] by the amount [i].
+ Beware, it is only defined if all the variables of v are greater than i
+ *)
+val decr_var : int -> t -> t
+
+(** [incr_var i v] increments the variables of the vector [v] by the amount [i].
+ *)
+val incr_var : int -> t -> t
+
+(** [gcd v] returns gcd(num(c),num(a1),...,num(an)) where num extracts
+ the numerator of a rational value. *)
+val gcd : t -> Big_int.big_int
+
+(** [normalise v] returns a vector with only integer coefficients *)
+val normalise : t -> t
+
+
+(** {1 Linear arithmetics} *)
+
+(** [add v1 v2] is vector addition.
+ @param v1 is of the form c +a1.x1 +...+an.xn
+ @param v2 is of the form c'+a1'.x1 +...+an'.xn
+ @return c1+c1'+ (a1+a1').x1 + ... + (an+an').xn
+ *)
+val add : t -> t -> t
+
+(** [mul a v] is vector multiplication of vector [v] by a scalar [a].
+ @return a.v = a.c+a.a1.x1+...+a.an.xn *)
+val mul : num -> t -> t
+
+(** [mul_add c1 v1 c2 v2] returns the linear combination c1.v1+c2.v2 *)
+val mul_add : num -> t -> num -> t -> t
+
+(** [div c1 v1] returns the mutiplication by the inverse of c1 i.e (1/c1).v1 *)
+val div : num -> t -> t
+
+(** [uminus v] @return -v the opposite vector of v i.e. (-1).v *)
+val uminus : t -> t
+
+(** {1 Iterators} *)
+
+(** [fold f acc v] returns f (f (f acc 0 c ) x1 a1 ) ... xn an *)
+val fold : ('acc -> var -> num -> 'acc) -> 'acc -> t -> 'acc
+
+(** [fold_error f acc v] is the same as
+ [fold (fun acc x i -> match acc with None -> None | Some acc' -> f acc' x i) (Some acc) v]
+ but with early exit...
+ *)
+val fold_error : ('acc -> var -> num -> 'acc option) -> 'acc -> t -> 'acc option
+
+(** [find f v] returns the first [f xi ai] such that [f xi ai <> None].
+ If no such xi ai exists, it returns None *)
+val find : (var -> num -> 'c option) -> t -> 'c option
+
+(** [for_all p v] returns /\_{i>=0} (f xi ai) *)
+val for_all : (var -> num -> bool) -> t -> bool
+
+(** [exists2 p v v'] returns Some(xi,ai,ai')
+ if p(xi,ai,ai') holds and ai,ai' <> 0.
+ It returns None if no such pair of coefficient exists. *)
+val exists2 : (num -> num -> bool) -> t -> t -> (var * num * num) option
+
+(** [dotproduct v1 v2] is the dot product of v1 and v2. *)
+val dotproduct : t -> t -> num
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-31 15:09:54 +0000