diff options
Diffstat (limited to 'contrib/micromega')
| -rw-r--r-- | contrib/micromega/Psatz.v (renamed from contrib/micromega/Micromegatac.v) | 42 | ||||
| -rw-r--r-- | contrib/micromega/QMicromega.v | 72 | ||||
| -rw-r--r-- | contrib/micromega/RMicromega.v | 36 | ||||
| -rw-r--r-- | contrib/micromega/ZMicromega.v | 11 | ||||
| -rw-r--r-- | contrib/micromega/certificate.ml | 342 | ||||
| -rw-r--r-- | contrib/micromega/coq_micromega.ml | 334 | ||||
| -rw-r--r-- | contrib/micromega/csdpcert.ml | 158 | ||||
| -rw-r--r-- | contrib/micromega/g_micromega.ml4 | 39 |
8 files changed, 492 insertions, 542 deletions
diff --git a/contrib/micromega/Micromegatac.v b/contrib/micromega/Psatz.v index 13c7eace..b2dd9910 100644 --- a/contrib/micromega/Micromegatac.v +++ b/contrib/micromega/Psatz.v @@ -23,57 +23,53 @@ Require Export RingMicromega. Require Import VarMap. Require Tauto. -Ltac micromegac dom d := +Ltac xpsatz dom d := let tac := lazymatch dom with | Z => - micromegap d ; + (sos_Z || psatz_Z d) ; intros __wit __varmap __ff ; change (Tauto.eval_f (Zeval_formula (@find Z Z0 __varmap)) __ff) ; apply (ZTautoChecker_sound __ff __wit); vm_compute ; reflexivity | R => - rmicromegap d ; + (sos_R || psatz_R d) ; intros __wit __varmap __ff ; change (Tauto.eval_f (Reval_formula (@find R 0%R __varmap)) __ff) ; apply (RTautoChecker_sound __ff __wit); vm_compute ; reflexivity + | Q => + (sos_Q || psatz_Q d) ; + intros __wit __varmap __ff ; + change (Tauto.eval_f (Qeval_formula (@find Q 0%Q __varmap)) __ff) ; + apply (QTautoChecker_sound __ff __wit); vm_compute ; reflexivity | _ => fail "Unsupported domain" end in tac. -Tactic Notation "micromega" constr(dom) int_or_var(n) := micromegac dom n. -Tactic Notation "micromega" constr(dom) := micromegac dom ltac:-1. +Tactic Notation "psatz" constr(dom) int_or_var(n) := xpsatz dom n. +Tactic Notation "psatz" constr(dom) := xpsatz dom ltac:-1. -Ltac zfarkas := omicronp ; - intros __wit __varmap __ff ; - change (Tauto.eval_f (Zeval_formula (@find Z Z0 __varmap)) __ff) ; - apply (ZTautoChecker_sound __ff __wit); vm_compute ; reflexivity. - -Ltac omicron dom := +Ltac psatzl dom := let tac := lazymatch dom with | Z => - zomicronp ; + psatzl_Z ; intros __wit __varmap __ff ; change (Tauto.eval_f (Zeval_formula (@find Z Z0 __varmap)) __ff) ; apply (ZTautoChecker_sound __ff __wit); vm_compute ; reflexivity | Q => - qomicronp ; + psatzl_Q ; intros __wit __varmap __ff ; change (Tauto.eval_f (Qeval_formula (@find Q 0%Q __varmap)) __ff) ; apply (QTautoChecker_sound __ff __wit); vm_compute ; reflexivity | R => - romicronp ; + psatzl_R ; intros __wit __varmap __ff ; change (Tauto.eval_f (Reval_formula (@find R 0%R __varmap)) __ff) ; apply (RTautoChecker_sound __ff __wit); vm_compute ; reflexivity | _ => fail "Unsupported domain" end in tac. -Ltac sos dom := - let tac := lazymatch dom with - | Z => - sosp ; - intros __wit __varmap __ff ; - change (Tauto.eval_f (Zeval_formula (@find Z Z0 __varmap)) __ff) ; - apply (ZTautoChecker_sound __ff __wit); vm_compute ; reflexivity - | _ => fail "Unsupported domain" - end in tac. +Ltac lia := + xlia ; + intros __wit __varmap __ff ; + change (Tauto.eval_f (Zeval_formula (@find Z Z0 __varmap)) __ff) ; + apply (ZTautoChecker_sound __ff __wit); vm_compute ; reflexivity. diff --git a/contrib/micromega/QMicromega.v b/contrib/micromega/QMicromega.v index 9e95f6c4..c054f218 100644 --- a/contrib/micromega/QMicromega.v +++ b/contrib/micromega/QMicromega.v @@ -16,38 +16,11 @@ Require Import OrderedRing. Require Import RingMicromega. Require Import Refl. Require Import QArith. -Require Import Qring. - -(* Qsrt has been removed from the library ? *) -Definition Qsrt : ring_theory 0 1 Qplus Qmult Qminus Qopp Qeq. -Proof. - constructor. - exact Qplus_0_l. - exact Qplus_comm. - exact Qplus_assoc. - exact Qmult_1_l. - exact Qmult_comm. - exact Qmult_assoc. - exact Qmult_plus_distr_l. - reflexivity. - exact Qplus_opp_r. -Qed. - - -Add Ring Qring : Qsrt. - -Lemma Qmult_neutral : forall x , 0 * x == 0. -Proof. - intros. - compute. - reflexivity. -Qed. - -(* Is there any qarith database ? *) +Require Import Qfield. Lemma Qsor : SOR 0 1 Qplus Qmult Qminus Qopp Qeq Qle Qlt. Proof. - constructor; intros ; subst ; try (intuition (subst; auto with qarith)). + constructor; intros ; subst ; try (intuition (subst; auto with qarith)). apply Q_Setoid. rewrite H ; rewrite H0 ; reflexivity. rewrite H ; rewrite H0 ; reflexivity. @@ -67,45 +40,12 @@ Proof. destruct(Q_dec n m) as [[H1 |H1] | H1 ] ; tauto. apply (Qplus_le_compat p p n m (Qle_refl p) H). generalize (Qmult_lt_compat_r 0 n m H0 H). - rewrite Qmult_neutral. + rewrite Qmult_0_l. auto. compute in H. discriminate. Qed. -Definition Qeq_bool (p q : Q) : bool := Zeq_bool (Qnum p * ' Qden q)%Z (Qnum q * ' Qden p)%Z. - -Definition Qle_bool (x y : Q) : bool := Zle_bool (Qnum x * ' Qden y)%Z (Qnum y * ' Qden x)%Z. - -Require ZMicromega. - -Lemma Qeq_bool_ok : forall x y, Qeq_bool x y = true -> x == y. -Proof. - intros. - unfold Qeq_bool in H. - unfold Qeq. - apply (Zeqb_ok _ _ H). -Qed. - - -Lemma Qeq_bool_neq : forall x y, Qeq_bool x y = false -> ~ x == y. -Proof. - unfold Qeq_bool,Qeq. - red ; intros ; subst. - rewrite H0 in H. - apply (ZMicromega.Zeq_bool_neq _ _ H). - reflexivity. -Qed. - -Lemma Qle_bool_imp_le : forall x y : Q, Qle_bool x y = true -> x <= y. -Proof. - unfold Qle_bool, Qle. - intros. - apply Zle_bool_imp_le ; auto. -Qed. - - - Lemma QSORaddon : SORaddon 0 1 Qplus Qmult Qminus Qopp Qeq Qle (* ring elements *) @@ -115,7 +55,7 @@ Lemma QSORaddon : Proof. constructor. constructor ; intros ; try reflexivity. - apply Qeq_bool_ok ; auto. + apply Qeq_bool_eq; auto. constructor. reflexivity. intros x y. @@ -173,9 +113,9 @@ match o with | OpEq => Qeq | OpNEq => fun x y => ~ x == y | OpLe => Qle -| OpGe => Qge +| OpGe => fun x y => Qle y x | OpLt => Qlt -| OpGt => Qgt +| OpGt => fun x y => Qlt y x end. Definition Qeval_formula (e:PolEnv Q) (ff : Formula Q) := diff --git a/contrib/micromega/RMicromega.v b/contrib/micromega/RMicromega.v index ef28db32..7c6969c2 100644 --- a/contrib/micromega/RMicromega.v +++ b/contrib/micromega/RMicromega.v @@ -76,12 +76,12 @@ Proof. apply mult_IZR. apply Ropp_Ropp_IZR. apply IZR_eq. - apply Zeqb_ok ; auto. + apply Zeq_bool_eq ; auto. apply R_power_theory. intros x y. intro. apply IZR_neq. - apply ZMicromega.Zeq_bool_neq ; auto. + apply Zeq_bool_neq ; auto. intros. apply IZR_le. apply Zle_bool_imp_le. auto. Qed. @@ -97,9 +97,34 @@ Definition INZ (n:N) : R := Definition Reval_expr := eval_pexpr Rplus Rmult Rminus Ropp IZR Nnat.nat_of_N pow. -Definition Reval_formula := +Definition Reval_op2 (o:Op2) : R -> R -> Prop := + match o with + | OpEq => @eq R + | OpNEq => fun x y => ~ x = y + | OpLe => Rle + | OpGe => Rge + | OpLt => Rlt + | OpGt => Rgt + end. + + +Definition Reval_formula (e: PolEnv R) (ff : Formula Z) := + let (lhs,o,rhs) := ff in Reval_op2 o (Reval_expr e lhs) (Reval_expr e rhs). + +Definition Reval_formula' := eval_formula Rplus Rmult Rminus Ropp (@eq R) Rle Rlt IZR Nnat.nat_of_N pow. +Lemma Reval_formula_compat : forall env f, Reval_formula env f <-> Reval_formula' env f. +Proof. + intros. + unfold Reval_formula. + destruct f. + unfold Reval_formula'. + unfold Reval_expr. + split ; destruct Fop ; simpl ; auto. + apply Rge_le. + apply Rle_ge. +Qed. Definition Reval_nformula := eval_nformula 0 Rplus Rmult Rminus Ropp (@eq R) Rle Rlt IZR Nnat.nat_of_N pow. @@ -139,8 +164,9 @@ Proof. unfold RTautoChecker. apply (tauto_checker_sound Reval_formula Reval_nformula). apply Reval_nformula_dec. - intros. unfold Reval_formula. now apply (cnf_normalise_correct Rsor). - intros. unfold Reval_formula. now apply (cnf_negate_correct Rsor). + intros. rewrite Reval_formula_compat. + unfold Reval_formula'. now apply (cnf_normalise_correct Rsor). + intros. rewrite Reval_formula_compat. unfold Reval_formula. now apply (cnf_negate_correct Rsor). intros t w0. apply RWeakChecker_sound. Qed. diff --git a/contrib/micromega/ZMicromega.v b/contrib/micromega/ZMicromega.v index 94c83f73..0855925a 100644 --- a/contrib/micromega/ZMicromega.v +++ b/contrib/micromega/ZMicromega.v @@ -39,15 +39,6 @@ Proof. apply Zmult_lt_0_compat ; auto. Qed. -Lemma Zeq_bool_neq : forall x y, Zeq_bool x y = false -> x <> y. -Proof. - red ; intros. - subst. - unfold Zeq_bool in H. - rewrite Zcompare_refl in H. - discriminate. -Qed. - Lemma ZSORaddon : SORaddon 0 1 Zplus Zmult Zminus Zopp (@eq Z) Zle (* ring elements *) 0%Z 1%Z Zplus Zmult Zminus Zopp (* coefficients *) @@ -56,7 +47,7 @@ Lemma ZSORaddon : Proof. constructor. constructor ; intros ; try reflexivity. - apply Zeqb_ok ; auto. + apply Zeq_bool_eq ; auto. constructor. reflexivity. intros x y. diff --git a/contrib/micromega/certificate.ml b/contrib/micromega/certificate.ml index 88e882e6..f4efcd08 100644 --- a/contrib/micromega/certificate.ml +++ b/contrib/micromega/certificate.ml @@ -108,7 +108,7 @@ struct if compare_num v (Int 0) <> 0 then if Monomial.compare Monomial.const k = 0 - then Printf.fprintf o "%s " (string_of_num v) + then Printf.fprintf o "%s " (string_of_num v) else Printf.fprintf o "%s*%a " (string_of_num v) Monomial.pp k) p (* Get the coefficient of monomial mn *) @@ -304,8 +304,8 @@ let factorise_linear_cone c = (match pending with | None -> rebuild_cone l (Some e) | Some p -> (match factorise p e with - | None -> Mc.S_Add(p, rebuild_cone l (Some e)) - | Some f -> rebuild_cone l (Some f) ) + | None -> Mc.S_Add(p, rebuild_cone l (Some e)) + | Some f -> rebuild_cone l (Some f) ) ) in (rebuild_cone (List.sort Pervasives.compare (cone_list c [])) None) @@ -387,10 +387,10 @@ let build_linear_system l = (* 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)); + (List.map (fun (x,y) -> + match y with Mc.Strict -> + Big_int unit_big_int + | _ -> Big_int zero_big_int) l)); op = Ge ; cst = Big_int unit_big_int } in (* Add the positivity constraint *) {coeffs = Vect.from_list ([Big_int unit_big_int]) ; @@ -414,22 +414,22 @@ let make_certificate n_spec cert li = in let rec scalar_product cert l = match cert with - | [] -> Mc.S_Z - | c::cert -> match l with - | [] -> failwith "make_certificate(1)" - | i::l -> - let r = scalar_product cert l in - match compare_big_int c zero_big_int with - | -1 -> Mc.S_Add ( - Mc.S_Ideal (Mc.PEc ( bint_to_cst c), Mc.S_In (Ml2C.nat i)), - r) - | 0 -> r - | _ -> Mc.S_Add ( - Mc.S_Mult (Mc.S_Pos (bint_to_cst c), Mc.S_In (Ml2C.nat i)), - r) in + | [] -> Mc.S_Z + | c::cert -> match l with + | [] -> failwith "make_certificate(1)" + | i::l -> + let r = scalar_product cert l in + match compare_big_int c zero_big_int with + | -1 -> Mc.S_Add ( + Mc.S_Ideal (Mc.PEc ( bint_to_cst c), Mc.S_In (Ml2C.nat i)), + r) + | 0 -> r + | _ -> Mc.S_Add ( + Mc.S_Mult (Mc.S_Pos (bint_to_cst c), Mc.S_In (Ml2C.nat i)), + r) in Some ((factorise_linear_cone - (simplify_cone n_spec (Mc.S_Add (cst, scalar_product cert' li))))) + (simplify_cone n_spec (Mc.S_Add (cst, scalar_product cert' li))))) exception Found of Monomial.t @@ -444,7 +444,7 @@ let raw_certificate l = with x -> if debug then (Printf.printf "raw certificate %s" (Printexc.to_string x); - flush stdout) ; + flush stdout) ; None @@ -462,9 +462,9 @@ let linear_prover n_spec l = (fun (x,_) -> if snd' x = Mc.NonEqual then true else false) li in let l' = List.map (fun (c,i) -> let (Mc.Pair(x,y)) = c in - match y with - Mc.NonEqual -> failwith "cannot happen" - | y -> ((dev_form n_spec x, y),i)) l' in + match y with + Mc.NonEqual -> failwith "cannot happen" + | y -> ((dev_form n_spec x, y),i)) l' in simple_linear_prover n_spec l' @@ -513,106 +513,228 @@ let rec remove_assoc p x l = let eq x y = Vect.compare x y = 0 -(* Beurk... this code is a shame *) +let remove e l = List.fold_left (fun l x -> if eq x e then l else x::l) [] l -let rec zlinear_prover sys = xzlinear_prover [] sys -and xzlinear_prover enum l : (Mc.proofTerm option) = - match linear_prover z_spec l with - | Some prf -> Some (Mc.RatProof prf) - | None -> +(* The prover is (probably) incomplete -- + only searching for naive cutting planes *) + +let candidates sys = + let ll = List.fold_right ( + fun (Mc.Pair(e,k)) r -> + match k with + | Mc.NonStrict -> (dev_form z_spec e , Ge)::r + | Mc.Equal -> (dev_form z_spec e , Eq)::r + (* we already know the bound -- don't compute it again *) + | _ -> failwith "Cannot happen candidates") sys [] in + + let (sys,var_mn) = make_linear_system ll in + let vars = mapi (fun _ i -> Vect.set i (Int 1) Vect.null) var_mn in + (List.fold_left (fun l cstr -> + let gcd = Big_int (Vect.gcd cstr.coeffs) in + if gcd =/ (Int 1) && cstr.op = Eq + then l + else (Vect.mul (Int 1 // gcd) cstr.coeffs)::l) [] sys) @ vars + + +let rec xzlinear_prover planes sys = + match linear_prover z_spec sys with + | Some prf -> Some (Mc.RatProof prf) + | None -> (* find the candidate with the smallest range *) + (* Grrr - linear_prover is also calling 'make_linear_system' *) let ll = List.fold_right (fun (Mc.Pair(e,k)) r -> match k with Mc.NonEqual -> r | k -> (dev_form z_spec e , - match k with - | Mc.Strict | Mc.NonStrict -> Ge - (* Loss of precision -- weakness of fourier*) - | Mc.Equal -> Eq - | Mc.NonEqual -> failwith "Cannot happen") :: r) l [] in - - let (sys,var) = make_linear_system ll in - let res = - match Fourier.find_Q_interval sys with - | Some(i,x,j) -> if i =/ j - then Some(i,Vect.set x (Int 1) Vect.null,i) else None - | None -> None in - let res = match res with - | None -> - begin - let candidates = List.fold_right - (fun cstr acc -> - let gcd = Big_int (Vect.gcd cstr.coeffs) in - let vect = Vect.mul (Int 1 // gcd) cstr.coeffs in - if mem eq vect enum then acc - else ((vect,Fourier.optimise vect sys)::acc)) sys [] in - let candidates = List.fold_left (fun l (x,i) -> - match i with - None -> (x,Empty)::l - | Some i -> (x,i)::l) [] (candidates) in - match List.fold_left (fun (x1,i1) (x2,i2) -> - if smaller_itv i1 i2 - then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates - with - | (i,Empty) -> None - | (x,Itv(Some i, Some j)) -> Some(i,x,j) - | (x,Point n) -> Some(n,x,n) - | x -> match Fourier.find_Q_interval sys with - | None -> None - | Some(i,x,j) -> - if i =/ j - then Some(i,Vect.set x (Int 1) Vect.null,i) - else None - end - | _ -> res in - - match res with + match k with + Mc.NonStrict -> Ge + | Mc.Equal -> Eq + | Mc.Strict | Mc.NonEqual -> failwith "Cannot happen") :: r) sys [] in + let (ll,var) = make_linear_system ll in + let candidates = List.fold_left (fun acc vect -> + match Fourier.optimise vect ll with + | None -> acc + | Some i -> +(* Printf.printf "%s in %s\n" (Vect.string vect) (string_of_intrvl i) ; *) + flush stdout ; + (vect,i) ::acc) [] planes in + + let smallest_interval = + match List.fold_left (fun (x1,i1) (x2,i2) -> + if smaller_itv i1 i2 + then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates + with + | (x,Itv(Some i, Some j)) -> Some(i,x,j) + | (x,Point n) -> Some(n,x,n) + | x -> None (* This might be a cutting plane *) + in + match smallest_interval with | Some (lb,e,ub) -> - let (lbn,lbd) = - (Ml2C.bigint (sub_big_int (numerator lb) unit_big_int), - Ml2C.bigint (denominator lb)) in - let (ubn,ubd) = - (Ml2C.bigint (add_big_int unit_big_int (numerator ub)) , - Ml2C.bigint (denominator ub)) in - let expr = list_to_polynomial var (Vect.to_list e) in - (match - (*x <= ub -> x > ub *) - linear_prover z_spec - (Mc.Pair(pplus (pmult (pconst ubd) expr) (popp (pconst ubn)), - Mc.NonStrict) :: l), - (* lb <= x -> lb > x *) - linear_prover z_spec - (Mc.Pair( pplus (popp (pmult (pconst lbd) expr)) (pconst lbn) , - Mc.NonStrict)::l) - with - | Some cub , Some clb -> - (match zlinear_enum (e::enum) expr - (ceiling_num lb) (floor_num ub) l - with - | None -> None - | Some prf -> - Some (Mc.EnumProof(Ml2C.q lb,expr,Ml2C.q ub,clb,cub,prf))) - | _ -> None - ) + let (lbn,lbd) = + (Ml2C.bigint (sub_big_int (numerator lb) unit_big_int), + Ml2C.bigint (denominator lb)) in + let (ubn,ubd) = + (Ml2C.bigint (add_big_int unit_big_int (numerator ub)) , + Ml2C.bigint (denominator ub)) in + let expr = list_to_polynomial var (Vect.to_list e) in + (match + (*x <= ub -> x > ub *) + linear_prover z_spec + (Mc.Pair(pplus (pmult (pconst ubd) expr) (popp (pconst ubn)), + Mc.NonStrict) :: sys), + (* lb <= x -> lb > x *) + linear_prover z_spec + (Mc.Pair( pplus (popp (pmult (pconst lbd) expr)) (pconst lbn) , + Mc.NonStrict)::sys) + with + | Some cub , Some clb -> + (match zlinear_enum (remove e planes) expr + (ceiling_num lb) (floor_num ub) sys + with + | None -> None + | Some prf -> + Some (Mc.EnumProof(Ml2C.q lb,expr,Ml2C.q ub,clb,cub,prf))) + | _ -> None + ) | _ -> None -and xzlinear_enum enum expr clb cub l = +and zlinear_enum planes expr clb cub l = if clb >/ cub then Some Mc.Nil else let pexpr = pplus (popp (pconst (Ml2C.bigint (numerator clb)))) expr in let sys' = (Mc.Pair(pexpr, Mc.Equal))::l in - match xzlinear_prover enum sys' with + (*let enum = *) + match xzlinear_prover planes sys' with | None -> if debug then print_string "zlp?"; None | Some prf -> if debug then print_string "zlp!"; - match zlinear_enum enum expr (clb +/ (Int 1)) cub l with + match zlinear_enum planes expr (clb +/ (Int 1)) cub l with | None -> None | Some prfl -> Some (Mc.Cons(prf,prfl)) -and zlinear_enum enum expr clb cub l = - let res = xzlinear_enum enum expr clb cub l in - if debug then Printf.printf "zlinear_enum %s %s -> %s\n" - (string_of_num clb) - (string_of_num cub) - (match res with - | None -> "None" - | Some r -> "Some") ; res +let zlinear_prover sys = + let candidates = candidates sys in + (* Printf.printf "candidates %d" (List.length candidates) ; *) + xzlinear_prover candidates sys + +open Sos + +let rec scale_term t = + match t with + | Zero -> unit_big_int , Zero + | Const n -> (denominator n) , Const (Big_int (numerator n)) + | Var n -> unit_big_int , Var n + | Inv _ -> failwith "scale_term : not implemented" + | Opp t -> let s, t = scale_term t in s, Opp t + | Add(t1,t2) -> let s1,y1 = scale_term t1 and s2,y2 = scale_term t2 in + let g = gcd_big_int s1 s2 in + let s1' = div_big_int s1 g in + let s2' = div_big_int s2 g in + let e = mult_big_int g (mult_big_int s1' s2') in + if (compare_big_int e unit_big_int) = 0 + then (unit_big_int, Add (y1,y2)) + else e, Add (Mul(Const (Big_int s2'), y1), + Mul (Const (Big_int s1'), y2)) + | Sub _ -> failwith "scale term: not implemented" + | Mul(y,z) -> let s1,y1 = scale_term y and s2,y2 = scale_term z in + mult_big_int s1 s2 , Mul (y1, y2) + | Pow(t,n) -> let s,t = scale_term t in + power_big_int_positive_int s n , Pow(t,n) + | _ -> failwith "scale_term : not implemented" + +let scale_term t = + let (s,t') = scale_term t in + s,t' + + +let get_index_of_ith_match f i l = + let rec get j res l = + match l with + | [] -> failwith "bad index" + | e::l -> if f e + then + (if j = i then res else get (j+1) (res+1) l ) + else get j (res+1) l in + get 0 0 l + + +let rec scale_certificate pos = match pos with + | Axiom_eq i -> unit_big_int , Axiom_eq i + | Axiom_le i -> unit_big_int , Axiom_le i + | Axiom_lt i -> unit_big_int , Axiom_lt i + | Monoid l -> unit_big_int , Monoid l + | Rational_eq n -> (denominator n) , Rational_eq (Big_int (numerator n)) + | Rational_le n -> (denominator n) , Rational_le (Big_int (numerator n)) + | Rational_lt n -> (denominator n) , Rational_lt (Big_int (numerator n)) + | Square t -> let s,t' = scale_term t in + mult_big_int s s , Square t' + | Eqmul (t, y) -> let s1,y1 = scale_term t and s2,y2 = scale_certificate y in + mult_big_int s1 s2 , Eqmul (y1,y2) + | Sum (y, z) -> let s1,y1 = scale_certificate y + and s2,y2 = scale_certificate z in + let g = gcd_big_int s1 s2 in + let s1' = div_big_int s1 g in + let s2' = div_big_int s2 g in + mult_big_int g (mult_big_int s1' s2'), + Sum (Product(Rational_le (Big_int s2'), y1), + Product (Rational_le (Big_int s1'), y2)) + | Product (y, z) -> + let s1,y1 = scale_certificate y and s2,y2 = scale_certificate z in + mult_big_int s1 s2 , Product (y1,y2) + + +open Micromega + let rec term_to_q_expr = function + | Const n -> PEc (Ml2C.q n) + | Zero -> PEc ( Ml2C.q (Int 0)) + | Var s -> PEX (Ml2C.index + (int_of_string (String.sub s 1 (String.length s - 1)))) + | Mul(p1,p2) -> PEmul(term_to_q_expr p1, term_to_q_expr p2) + | Add(p1,p2) -> PEadd(term_to_q_expr p1, term_to_q_expr p2) + | Opp p -> PEopp (term_to_q_expr p) + | Pow(t,n) -> PEpow (term_to_q_expr t,Ml2C.n n) + | Sub(t1,t2) -> PEsub (term_to_q_expr t1, term_to_q_expr t2) + | _ -> failwith "term_to_q_expr: not implemented" + +let q_cert_of_pos pos = + let rec _cert_of_pos = function + Axiom_eq i -> Mc.S_In (Ml2C.nat i) + | Axiom_le i -> Mc.S_In (Ml2C.nat i) + | Axiom_lt i -> Mc.S_In (Ml2C.nat i) + | Monoid l -> Mc.S_Monoid (Ml2C.list Ml2C.nat l) + | Rational_eq n | Rational_le n | Rational_lt n -> + if compare_num n (Int 0) = 0 then Mc.S_Z else + Mc.S_Pos (Ml2C.q n) + | Square t -> Mc.S_Square (term_to_q_expr t) + | Eqmul (t, y) -> Mc.S_Ideal(term_to_q_expr t, _cert_of_pos y) + | Sum (y, z) -> Mc.S_Add (_cert_of_pos y, _cert_of_pos z) + | Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in + simplify_cone q_spec (_cert_of_pos pos) + + + let rec term_to_z_expr = function + | Const n -> PEc (Ml2C.bigint (big_int_of_num n)) + | Zero -> PEc ( Z0) + | Var s -> PEX (Ml2C.index + (int_of_string (String.sub s 1 (String.length s - 1)))) + | Mul(p1,p2) -> PEmul(term_to_z_expr p1, term_to_z_expr p2) + | Add(p1,p2) -> PEadd(term_to_z_expr p1, term_to_z_expr p2) + | Opp p -> PEopp (term_to_z_expr p) + | Pow(t,n) -> PEpow (term_to_z_expr t,Ml2C.n n) + | Sub(t1,t2) -> PEsub (term_to_z_expr t1, term_to_z_expr t2) + | _ -> failwith "term_to_z_expr: not implemented" + +let z_cert_of_pos pos = + let s,pos = (scale_certificate pos) in + let rec _cert_of_pos = function + Axiom_eq i -> Mc.S_In (Ml2C.nat i) + | Axiom_le i -> Mc.S_In (Ml2C.nat i) + | Axiom_lt i -> Mc.S_In (Ml2C.nat i) + | Monoid l -> Mc.S_Monoid (Ml2C.list Ml2C.nat l) + | Rational_eq n | Rational_le n | Rational_lt n -> + if compare_num n (Int 0) = 0 then Mc.S_Z else + Mc.S_Pos (Ml2C.bigint (big_int_of_num n)) + | Square t -> Mc.S_Square (term_to_z_expr t) + | Eqmul (t, y) -> Mc.S_Ideal(term_to_z_expr t, _cert_of_pos y) + | Sum (y, z) -> Mc.S_Add (_cert_of_pos y, _cert_of_pos z) + | Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in + simplify_cone z_spec (_cert_of_pos pos) diff --git a/contrib/micromega/coq_micromega.ml b/contrib/micromega/coq_micromega.ml index 29e2a183..02ff61a1 100644 --- a/contrib/micromega/coq_micromega.ml +++ b/contrib/micromega/coq_micromega.ml @@ -106,6 +106,7 @@ struct ["Coq" ; "micromega" ; "EnvRing"]; ["Coq";"QArith"; "QArith_base"]; ["Coq";"Reals" ; "Rdefinitions"]; + ["Coq";"Reals" ; "Rpow_def"]; ["LRing_normalise"]] let constant = gen_constant_in_modules "ZMicromega" coq_modules @@ -163,6 +164,9 @@ struct let coq_Qmake = lazy (constant "Qmake") + let coq_R0 = lazy (constant "R0") + let coq_R1 = lazy (constant "R1") + let coq_proofTerm = lazy (constant "ProofTerm") let coq_ratProof = lazy (constant "RatProof") @@ -179,10 +183,36 @@ struct let coq_Zminus = lazy (constant "Zminus") let coq_Zopp = lazy (constant "Zopp") let coq_Zmult = lazy (constant "Zmult") + let coq_Zpower = lazy (constant "Zpower") let coq_N_of_Z = lazy (gen_constant_in_modules "ZArithRing" [["Coq";"setoid_ring";"ZArithRing"]] "N_of_Z") + let coq_Qgt = lazy (constant "Qgt") + let coq_Qge = lazy (constant "Qge") + let coq_Qle = lazy (constant "Qle") + let coq_Qlt = lazy (constant "Qlt") + let coq_Qeq = lazy (constant "Qeq") + + + let coq_Qplus = lazy (constant "Qplus") + let coq_Qminus = lazy (constant "Qminus") + let coq_Qopp = lazy (constant "Qopp") + let coq_Qmult = lazy (constant "Qmult") + let coq_Qpower = lazy (constant "Qpower") + + + let coq_Rgt = lazy (constant "Rgt") + let coq_Rge = lazy (constant "Rge") + let coq_Rle = lazy (constant "Rle") + let coq_Rlt = lazy (constant "Rlt") + + let coq_Rplus = lazy (constant "Rplus") + let coq_Rminus = lazy (constant "Rminus") + let coq_Ropp = lazy (constant "Ropp") + let coq_Rmult = lazy (constant "Rmult") + let coq_Rpower = lazy (constant "pow") + let coq_PEX = lazy (constant "PEX" ) let coq_PEc = lazy (constant"PEc") @@ -225,6 +255,7 @@ struct (gen_constant_in_modules "RingMicromega" [["Coq" ; "micromega" ; "RingMicromega"] ; ["RingMicromega"] ] "Formula") + type parse_error = | Ukn | BadStr of string @@ -347,16 +378,11 @@ let dump_q q = let parse_q term = match Term.kind_of_term term with - | Term.App(c, args) -> - ( - match Term.kind_of_term c with - Term.Construct((n,j),i) -> - if Names.string_of_kn n = "Coq.QArith.QArith_base#<>#Q" - then {Mc.qnum = parse_z args.(0) ; Mc.qden = parse_positive args.(1) } + | Term.App(c, args) -> if c = Lazy.force coq_Qmake then + {Mc.qnum = parse_z args.(0) ; Mc.qden = parse_positive args.(1) } else raise ParseError - | _ -> raise ParseError - ) - | _ -> raise ParseError + | _ -> raise ParseError + let rec parse_list parse_elt term = let (i,c) = get_left_construct term in @@ -466,19 +492,6 @@ let parse_q term = pp_cone o e - - - let rec parse_op term = - let (i,c) = get_left_construct term in - match i with - | 1 -> Mc.OpEq - | 2 -> Mc.OpLe - | 3 -> Mc.OpGe - | 4 -> Mc.OpGt - | 5 -> Mc.OpLt - | i -> raise ParseError - - let rec dump_op = function | Mc.OpEq-> Lazy.force coq_OpEq | Mc.OpNEq-> Lazy.force coq_OpNEq @@ -510,68 +523,52 @@ let parse_q term = dump_op o ; dump_expr typ dump_constant e2|]) + let assoc_const x l = + try + snd (List.find (fun (x',y) -> x = Lazy.force x') l) + with + Not_found -> raise ParseError + + let zop_table = [ + coq_Zgt, Mc.OpGt ; + coq_Zge, Mc.OpGe ; + coq_Zlt, Mc.OpLt ; + coq_Zle, Mc.OpLe ] + + let rop_table = [ + coq_Rgt, Mc.OpGt ; + coq_Rge, Mc.OpGe ; + coq_Rlt, Mc.OpLt ; + coq_Rle, Mc.OpLe ] + + let qop_table = [ + coq_Qlt, Mc.OpLt ; + coq_Qle, Mc.OpLe ; + coq_Qeq, Mc.OpEq + ] let parse_zop (op,args) = match kind_of_term op with - | Const x -> - (match Names.string_of_con x with - | "Coq.ZArith.BinInt#<>#Zgt" -> (Mc.OpGt, args.(0), args.(1)) - | "Coq.ZArith.BinInt#<>#Zge" -> (Mc.OpGe, args.(0), args.(1)) - | "Coq.ZArith.BinInt#<>#Zlt" -> (Mc.OpLt, args.(0), args.(1)) - | "Coq.ZArith.BinInt#<>#Zle" -> (Mc.OpLe, args.(0), args.(1)) - (*| "Coq.Init.Logic#<>#not" -> Mc.OpNEq (* for backward compat *)*) - | s -> raise ParseError - ) + | Const x -> (assoc_const op zop_table, args.(0) , args.(1)) | Ind(n,0) -> - (match Names.string_of_kn n with - | "Coq.Init.Logic#<>#eq" -> - if args.(0) <> Lazy.force coq_Z - then raise ParseError - else (Mc.OpEq, args.(1), args.(2)) - | _ -> raise ParseError) + if op = Lazy.force coq_Eq && args.(0) = Lazy.force coq_Z + then (Mc.OpEq, args.(1), args.(2)) + else raise ParseError | _ -> failwith "parse_zop" let parse_rop (op,args) = - try match kind_of_term op with - | Const x -> - (match Names.string_of_con x with - | "Coq.Reals.Rdefinitions#<>#Rgt" -> (Mc.OpGt, args.(0), args.(1)) - | "Coq.Reals.Rdefinitions#<>#Rge" -> (Mc.OpGe, args.(0), args.(1)) - | "Coq.Reals.Rdefinitions#<>#Rlt" -> (Mc.OpLt, args.(0), args.(1)) - | "Coq.Reals.Rdefinitions#<>#Rle" -> (Mc.OpLe, args.(0), args.(1)) - (*| "Coq.Init.Logic#<>#not"-> Mc.OpNEq (* for backward compat *)*) - | s -> raise ParseError - ) + | Const x -> (assoc_const op rop_table, args.(0) , args.(1)) | Ind(n,0) -> - (match Names.string_of_kn n with - | "Coq.Init.Logic#<>#eq" -> - (* if args.(0) <> Lazy.force coq_R - then raise ParseError - else*) (Mc.OpEq, args.(1), args.(2)) - | _ -> raise ParseError) - | _ -> failwith "parse_rop" - with x -> - (Pp.pp (Pp.str "parse_rop failure ") ; - Pp.pp (Printer.prterm op) ; Pp.pp_flush ()) - ; raise x - + if op = Lazy.force coq_Eq && args.(0) = Lazy.force coq_R + then (Mc.OpEq, args.(1), args.(2)) + else raise ParseError + | _ -> failwith "parse_zop" let parse_qop (op,args) = - ( - (match kind_of_term op with - | Const x -> - (match Names.string_of_con x with - | "Coq.QArith.QArith_base#<>#Qgt" -> Mc.OpGt - | "Coq.QArith.QArith_base#<>#Qge" -> Mc.OpGe - | "Coq.QArith.QArith_base#<>#Qlt" -> Mc.OpLt - | "Coq.QArith.QArith_base#<>#Qle" -> Mc.OpLe - | "Coq.QArith.QArith_base#<>#Qeq" -> Mc.OpEq - | s -> raise ParseError - ) - | _ -> failwith "parse_zop") , args.(0) , args.(1)) + (assoc_const op qop_table, args.(0) , args.(1)) module Env = @@ -612,6 +609,14 @@ let parse_q term = | Ukn of string + let assoc_ops x l = + try + snd (List.find (fun (x',y) -> x = Lazy.force x') l) + with + Not_found -> Ukn "Oups" + + + let parse_expr parse_constant parse_exp ops_spec env term = if debug then (Pp.pp (Pp.str "parse_expr: "); @@ -634,7 +639,7 @@ let parse_q term = ( match kind_of_term t with | Const c -> - ( match ops_spec (Names.string_of_con c) with + ( match assoc_ops t ops_spec with | Binop f -> combine env f (args.(0),args.(1)) | Opp -> let (expr,env) = parse_expr env args.(0) in (Mc.PEopp expr, env) @@ -653,29 +658,29 @@ let parse_q term = parse_expr env term -let zop_spec = function - | "Coq.ZArith.BinInt#<>#Zplus" -> Binop (fun x y -> Mc.PEadd(x,y)) - | "Coq.ZArith.BinInt#<>#Zminus" -> Binop (fun x y -> Mc.PEsub(x,y)) - | "Coq.ZArith.BinInt#<>#Zmult" -> Binop (fun x y -> Mc.PEmul (x,y)) - | "Coq.ZArith.BinInt#<>#Zopp" -> Opp - | "Coq.ZArith.Zpow_def#<>#Zpower" -> Power - | s -> Ukn s + let zop_spec = + [ + coq_Zplus , Binop (fun x y -> Mc.PEadd(x,y)) ; + coq_Zminus , Binop (fun x y -> Mc.PEsub(x,y)) ; + coq_Zmult , Binop (fun x y -> Mc.PEmul (x,y)) ; + coq_Zopp , Opp ; + coq_Zpower , Power] -let qop_spec = function - | "Coq.QArith.QArith_base#<>#Qplus" -> Binop (fun x y -> Mc.PEadd(x,y)) - | "Coq.QArith.QArith_base#<>#Qminus" -> Binop (fun x y -> Mc.PEsub(x,y)) - | "Coq.QArith.QArith_base#<>#Qmult" -> Binop (fun x y -> Mc.PEmul (x,y)) - | "Coq.QArith.QArith_base#<>#Qopp" -> Opp - | "Coq.QArith.QArith_base#<>#Qpower" -> Power - | s -> Ukn s +let qop_spec = + [ + coq_Qplus , Binop (fun x y -> Mc.PEadd(x,y)) ; + coq_Qminus , Binop (fun x y -> Mc.PEsub(x,y)) ; + coq_Qmult , Binop (fun x y -> Mc.PEmul (x,y)) ; + coq_Qopp , Opp ; + coq_Qpower , Power] -let rop_spec = function - | "Coq.Reals.Rdefinitions#<>#Rplus" -> Binop (fun x y -> Mc.PEadd(x,y)) - | "Coq.Reals.Rdefinitions#<>#Rminus" -> Binop (fun x y -> Mc.PEsub(x,y)) - | "Coq.Reals.Rdefinitions#<>#Rmult" -> Binop (fun x y -> Mc.PEmul (x,y)) - | "Coq.Reals.Rdefinitions#<>#Ropp" -> Opp - | "Coq.Reals.Rpow_def#<>#pow" -> Power - | s -> Ukn s +let rop_spec = + [ + coq_Rplus , Binop (fun x y -> Mc.PEadd(x,y)) ; + coq_Rminus , Binop (fun x y -> Mc.PEsub(x,y)) ; + coq_Rmult , Binop (fun x y -> Mc.PEmul (x,y)) ; + coq_Ropp , Opp ; + coq_Rpower , Power] @@ -691,12 +696,12 @@ let rconstant term = Pp.pp (Pp.str "rconstant: "); Pp.pp (Printer.prterm term); Pp.pp_flush ()); match Term.kind_of_term term with - | Const x -> - (match Names.string_of_con x with - | "Coq.Reals.Rdefinitions#<>#R0" -> Mc.Z0 - | "Coq.Reals.Rdefinitions#<>#R1" -> Mc.Zpos Mc.XH - | _ -> raise ParseError - ) + | Const x -> + if term = Lazy.force coq_R0 + then Mc.Z0 + else if term = Lazy.force coq_R1 + then Mc.Zpos Mc.XH + else raise ParseError | _ -> raise ParseError @@ -731,23 +736,6 @@ let parse_rexpr = (* generic parsing of arithmetic expressions *) - let rec parse_conj parse_arith env term = - match kind_of_term term with - | App(l,rst) -> - (match kind_of_term l with - | Ind (n,_) -> - ( match Names.string_of_kn n with - | "Coq.Init.Logic#<>#and" -> - let (e1,env) = parse_arith env rst.(0) in - let (e2,env) = parse_conj parse_arith env rst.(1) in - (Mc.Cons(e1,e2),env) - | _ -> (* This might be an equality *) - let (e,env) = parse_arith env term in - (Mc.Cons(e,Mc.Nil),env)) - | _ -> (* This is an arithmetic expression *) - let (e,env) = parse_arith env term in - (Mc.Cons(e,Mc.Nil),env)) - | _ -> failwith "parse_conj(2)" @@ -850,46 +838,6 @@ let parse_rexpr = xdump f - (* Backward compat *) - - let rec parse_concl parse_arith env term = - match kind_of_term term with - | Prod(_,expr,rst) -> (* a -> b *) - let (lhs,rhs,env) = parse_concl parse_arith env rst in - let (e,env) = parse_arith env expr in - (Mc.Cons(e,lhs),rhs,env) - | App(_,_) -> - let (conj, env) = parse_conj parse_arith env term in - (Mc.Nil,conj,env) - | Ind(n,_) -> - (match (Names.string_of_kn n) with - | "Coq.Init.Logic#<>#False" -> (Mc.Nil,Mc.Nil,env) - | s -> - print_string s ; flush stdout; - failwith "parse_concl") - | _ -> failwith "parse_concl" - - - let rec parse_hyps parse_arith env goal_hyps hyps = - match hyps with - | [] -> ([],goal_hyps,env) - | (i,t)::l -> - let (li,lt,env) = parse_hyps parse_arith env goal_hyps l in - try - let (c,env) = parse_arith env t in - (i::li, Mc.Cons(c,lt), env) - with x -> - (*(if debug then Printf.printf "parse_arith : %s\n" x);*) - (li,lt,env) - - - let parse_goal parse_arith env hyps term = - try - let (lhs,rhs,env) = parse_concl parse_arith env term in - let (li,lt,env) = parse_hyps parse_arith env lhs hyps in - (li,lt,rhs,env) - with Failure x -> raise ParseError - (* backward compat *) (* ! reverse the list of bindings *) @@ -1235,8 +1183,8 @@ let lift_ratproof prover l = | Some c -> Some (Mc.RatProof c) -type csdpcert = Certificate.Mc.z Certificate.Mc.coneMember option -type micromega_polys = (Micromega.z Mc.pExpr, Mc.op1) Micromega.prod list +type csdpcert = Sos.positivstellensatz option +type micromega_polys = (Micromega.q Mc.pExpr, Mc.op1) Micromega.prod list type provername = string * int option let call_csdpcert provername poly = @@ -1255,36 +1203,84 @@ let call_csdpcert provername poly = close_in ch_from; Sys.remove tmp_from; cert -let omicron gl = +let rec z_to_q_expr e = + match e with + | Mc.PEc z -> Mc.PEc {Mc.qnum = z ; Mc.qden = Mc.XH} + | Mc.PEX x -> Mc.PEX x + | Mc.PEadd(e1,e2) -> Mc.PEadd(z_to_q_expr e1, z_to_q_expr e2) + | Mc.PEsub(e1,e2) -> Mc.PEsub(z_to_q_expr e1, z_to_q_expr e2) + | Mc.PEmul(e1,e2) -> Mc.PEmul(z_to_q_expr e1, z_to_q_expr e2) + | Mc.PEopp(e) -> Mc.PEopp(z_to_q_expr e) + | Mc.PEpow(e,n) -> Mc.PEpow(z_to_q_expr e,n) + + +let call_csdpcert_q provername poly = + match call_csdpcert provername poly with + | None -> None + | Some cert -> + let cert = Certificate.q_cert_of_pos cert in + match Mc.qWeakChecker (CamlToCoq.list (fun x -> x) poly) cert with + | Mc.True -> Some cert + | Mc.False -> (print_string "buggy certificate" ; flush stdout) ;None + + +let call_csdpcert_z provername poly = + let l = List.map (fun (Mc.Pair(e,o)) -> (Mc.Pair(z_to_q_expr e,o))) poly in + match call_csdpcert provername l with + | None -> None + | Some cert -> + let cert = Certificate.z_cert_of_pos cert in + match Mc.zWeakChecker (CamlToCoq.list (fun x -> x) poly) cert with + | Mc.True -> Some cert + | Mc.False -> (print_string "buggy certificate" ; flush stdout) ;None + + + + +let psatzl_Z gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec [lift_ratproof (Certificate.linear_prover Certificate.z_spec), "fourier refutation" ] gl -let qomicron gl = +let psatzl_Q gl = micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec [ Certificate.linear_prover Certificate.q_spec, "fourier refutation" ] gl -let romicron gl = +let psatz_Q i gl = + micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec + [ call_csdpcert_q ("real_nonlinear_prover", Some i), "fourier refutation" ] gl + +let psatzl_R gl = micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec [ Certificate.linear_prover Certificate.z_spec, "fourier refutation" ] gl -let rmicromega i gl = - micromega_gen parse_rarith Mc.negate Mc.normalise rz_domain_spec - [ call_csdpcert ("real_nonlinear_prover", Some i), "fourier refutation" ] gl +let psatz_R i gl = + micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec + [ call_csdpcert_z ("real_nonlinear_prover", Some i), "fourier refutation" ] gl -let micromega i gl = +let psatz_Z i gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec - [lift_ratproof (call_csdpcert ("real_nonlinear_prover",Some i)), + [lift_ratproof (call_csdpcert_z ("real_nonlinear_prover",Some i)), "fourier refutation" ] gl -let sos gl = +let sos_Z gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec - [lift_ratproof (call_csdpcert ("pure_sos", None)), "pure sos refutation"] gl + [lift_ratproof (call_csdpcert_z ("pure_sos", None)), "pure sos refutation"] gl + +let sos_Q gl = + micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec + [call_csdpcert_q ("pure_sos", None), "pure sos refutation"] gl + +let sos_R gl = + micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec + [call_csdpcert_z ("pure_sos", None), "pure sos refutation"] gl + + -let zomicron gl = +let xlia gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec [Certificate.zlinear_prover, "zprover"] gl diff --git a/contrib/micromega/csdpcert.ml b/contrib/micromega/csdpcert.ml index cfaf6ae1..e451a38f 100644 --- a/contrib/micromega/csdpcert.ml +++ b/contrib/micromega/csdpcert.ml @@ -27,7 +27,7 @@ struct open Mc let rec expr_to_term = function - | PEc z -> Const (Big_int (C2Ml.z_big_int z)) + | PEc z -> Const (C2Ml.q_to_num z) | PEX v -> Var ("x"^(string_of_int (C2Ml.index v))) | PEmul(p1,p2) -> let p1 = expr_to_term p1 in @@ -40,25 +40,15 @@ struct | PEopp p -> Opp (expr_to_term p) - let rec term_to_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_expr p1, term_to_expr p2) - | Add(p1,p2) -> PEadd(term_to_expr p1, term_to_expr p2) - | Opp p -> PEopp (term_to_expr p) - | Pow(t,n) -> PEpow (term_to_expr t,Ml2C.n n) - | Sub(t1,t2) -> PEsub (term_to_expr t1, term_to_expr t2) - | _ -> failwith "term_to_expr: not implemented" - - let term_to_expr e = + + +(* let term_to_expr e = let e' = term_to_expr e in if debug then Printf.printf "term_to_expr : %s - %s\n" (string_of_poly (poly_of_term e)) (string_of_poly (poly_of_term (expr_to_term e'))); - e' + e' *) end open M @@ -66,110 +56,8 @@ open M open List open Mutils -let rec scale_term t = - match t with - | Zero -> unit_big_int , Zero - | Const n -> (denominator n) , Const (Big_int (numerator n)) - | Var n -> unit_big_int , Var n - | Inv _ -> failwith "scale_term : not implemented" - | Opp t -> let s, t = scale_term t in s, Opp t - | Add(t1,t2) -> let s1,y1 = scale_term t1 and s2,y2 = scale_term t2 in - let g = gcd_big_int s1 s2 in - let s1' = div_big_int s1 g in - let s2' = div_big_int s2 g in - let e = mult_big_int g (mult_big_int s1' s2') in - if (compare_big_int e unit_big_int) = 0 - then (unit_big_int, Add (y1,y2)) - else e, Add (Mul(Const (Big_int s2'), y1), - Mul (Const (Big_int s1'), y2)) - | Sub _ -> failwith "scale term: not implemented" - | Mul(y,z) -> let s1,y1 = scale_term y and s2,y2 = scale_term z in - mult_big_int s1 s2 , Mul (y1, y2) - | Pow(t,n) -> let s,t = scale_term t in - power_big_int_positive_int s n , Pow(t,n) - | _ -> failwith "scale_term : not implemented" - -let scale_term t = - let (s,t') = scale_term t in - s,t' - - -let 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) - - -let is_eq = function Mc.Equal -> true | _ -> false -let is_le = function Mc.NonStrict -> true | _ -> false -let is_lt = function Mc.Strict -> true | _ -> false - -let get_index_of_ith_match f i l = - let rec get j res l = - match l with - | [] -> failwith "bad index" - | e::l -> if f e - then - (if j = i then res else get (j+1) (res+1) l ) - else get j (res+1) l in - get 0 0 l - - -let cert_of_pos eq le lt ll l pos = - let s,pos = (scale_certificate pos) in - let rec _cert_of_pos = function - Axiom_eq i -> let idx = get_index_of_ith_match is_eq i l in - Mc.S_In (Ml2C.nat idx) - | Axiom_le i -> let idx = get_index_of_ith_match is_le i l in - Mc.S_In (Ml2C.nat idx) - | Axiom_lt i -> let idx = get_index_of_ith_match is_lt i l in - Mc.S_In (Ml2C.nat idx) - | Monoid l -> Mc.S_Monoid (Ml2C.list Ml2C.nat l) - | Rational_eq n | Rational_le n | Rational_lt n -> - if compare_num n (Int 0) = 0 then Mc.S_Z else - Mc.S_Pos (Ml2C.bigint (big_int_of_num n)) - | Square t -> Mc.S_Square (term_to_expr t) - | Eqmul (t, y) -> Mc.S_Ideal(term_to_expr t, _cert_of_pos y) - | Sum (y, z) -> Mc.S_Add (_cert_of_pos y, _cert_of_pos z) - | Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in - s, Certificate.simplify_cone Certificate.z_spec (_cert_of_pos pos) - - -let term_of_cert l pos = - let l = List.map fst' l in - let rec _cert_of_pos = function - | Mc.S_In i -> expr_to_term (List.nth l (C2Ml.nat i)) - | Mc.S_Pos p -> Const (C2Ml.num p) - | Mc.S_Z -> Const (Int 0) - | Mc.S_Square t -> Mul(expr_to_term t, expr_to_term t) - | Mc.S_Monoid m -> List.fold_right - (fun x m -> Mul (expr_to_term (List.nth l (C2Ml.nat x)),m)) - (C2Ml.list (fun x -> x) m) (Const (Int 1)) - | Mc.S_Ideal (t, y) -> Mul(expr_to_term t, _cert_of_pos y) - | Mc.S_Add (y, z) -> Add (_cert_of_pos y, _cert_of_pos z) - | Mc.S_Mult (y, z) -> Mul (_cert_of_pos y, _cert_of_pos z) in - (_cert_of_pos pos) let rec canonical_sum_to_string = function s -> failwith "not implemented" @@ -208,22 +96,6 @@ let rec sets_of_list l = | e::l -> let s = sets_of_list l in s@(List.map (fun s0 -> e::s0) s) -let cert_of_pos pos = - let s,pos = (scale_certificate pos) in - let rec _cert_of_pos = function - Axiom_eq i -> Mc.S_In (Ml2C.nat i) - | Axiom_le i -> Mc.S_In (Ml2C.nat i) - | Axiom_lt i -> Mc.S_In (Ml2C.nat i) - | Monoid l -> Mc.S_Monoid (Ml2C.list Ml2C.nat l) - | Rational_eq n | Rational_le n | Rational_lt n -> - if compare_num n (Int 0) = 0 then Mc.S_Z else - Mc.S_Pos (Ml2C.bigint (big_int_of_num n)) - | Square t -> Mc.S_Square (term_to_expr t) - | Eqmul (t, y) -> Mc.S_Ideal(term_to_expr t, _cert_of_pos y) - | Sum (y, z) -> Mc.S_Add (_cert_of_pos y, _cert_of_pos z) - | Product (y, z) -> Mc.S_Mult (_cert_of_pos y, _cert_of_pos z) in - s, Certificate.simplify_cone Certificate.z_spec (_cert_of_pos pos) - (* The exploration is probably not complete - for simple cases, it works... *) let real_nonlinear_prover d l = try @@ -272,15 +144,7 @@ let real_nonlinear_prover d l = let proof = list_fold_right_elements (fun s t -> Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) in - - let s,proof' = scale_certificate proof in - let cert = snd (cert_of_pos proof') in - if debug - then Printf.printf "cert poly : %s\n" - (string_of_poly (poly_of_term (term_of_cert l cert))); - match Mc.zWeakChecker (Ml2C.list (fun x -> x) l) cert with - | Mc.True -> Some cert - | Mc.False -> (print_string "buggy certificate" ; flush stdout) ;None + Some proof with | Sos.TooDeep -> None @@ -300,16 +164,16 @@ let pure_sos l = (term_of_poly p)), rst)) polys (Rational_lt (Int 0))) in let proof = Sum(Axiom_lt i, pos) in - let s,proof' = scale_certificate proof in - let cert = snd (cert_of_pos proof') in - Some cert +(* let s,proof' = scale_certificate proof in + let cert = snd (cert_of_pos proof') in *) + Some proof with | Not_found -> (* This is no strict inequality *) None | x -> None -type micromega_polys = (Micromega.z Mc.pExpr, Mc.op1) Micromega.prod list -type csdp_certificate = Certificate.Mc.z Certificate.Mc.coneMember option +type micromega_polys = (Micromega.q Mc.pExpr, Mc.op1) Micromega.prod list +type csdp_certificate = Sos.positivstellensatz option type provername = string * int option let main () = diff --git a/contrib/micromega/g_micromega.ml4 b/contrib/micromega/g_micromega.ml4 index 259b5d4b..50024e78 100644 --- a/contrib/micromega/g_micromega.ml4 +++ b/contrib/micromega/g_micromega.ml4 @@ -14,7 +14,7 @@ (*i camlp4deps: "parsing/grammar.cma" i*) -(* $Id: g_micromega.ml4 10947 2008-05-19 19:10:40Z herbelin $ *) +(* $Id: g_micromega.ml4 11306 2008-08-05 16:51:08Z notin $ *) open Quote open Ring @@ -26,34 +26,49 @@ let out_arg = function | ArgVar _ -> anomaly "Unevaluated or_var variable" | ArgArg x -> x -TACTIC EXTEND Micromega -| [ "micromegap" int_or_var(i) ] -> [ Coq_micromega.micromega (out_arg i) ] -| [ "micromegap" ] -> [ Coq_micromega.micromega (-1) ] +TACTIC EXTEND PsatzZ +| [ "psatz_Z" int_or_var(i) ] -> [ Coq_micromega.psatz_Z (out_arg i) ] +| [ "psatz_Z" ] -> [ Coq_micromega.psatz_Z (-1) ] END -TACTIC EXTEND Sos -[ "sosp" ] -> [ Coq_micromega.sos] +TACTIC EXTEND Sos_Z +| [ "sos_Z" ] -> [ Coq_micromega.sos_Z] + END + +TACTIC EXTEND Sos_Q +| [ "sos_Q" ] -> [ Coq_micromega.sos_Q] + END + +TACTIC EXTEND Sos_R +| [ "sos_R" ] -> [ Coq_micromega.sos_R] END TACTIC EXTEND Omicron -[ "omicronp" ] -> [ Coq_micromega.omicron] +[ "psatzl_Z" ] -> [ Coq_micromega.psatzl_Z] END TACTIC EXTEND QOmicron -[ "qomicronp" ] -> [ Coq_micromega.qomicron] +[ "psatzl_Q" ] -> [ Coq_micromega.psatzl_Q] END TACTIC EXTEND ZOmicron -[ "zomicronp" ] -> [ Coq_micromega.zomicron] +[ "xlia" ] -> [ Coq_micromega.xlia] END TACTIC EXTEND ROmicron -[ "romicronp" ] -> [ Coq_micromega.romicron] +[ "psatzl_R" ] -> [ Coq_micromega.psatzl_R] END TACTIC EXTEND RMicromega -| [ "rmicromegap" int_or_var(i) ] -> [ Coq_micromega.rmicromega (out_arg i) ] -| [ "rmicromegap" ] -> [ Coq_micromega.rmicromega (-1) ] +| [ "psatz_R" int_or_var(i) ] -> [ Coq_micromega.psatz_R (out_arg i) ] +| [ "psatz_R" ] -> [ Coq_micromega.psatz_R (-1) ] +END + + +TACTIC EXTEND QMicromega +| [ "psatz_Q" int_or_var(i) ] -> [ Coq_micromega.psatz_Q (out_arg i) ] +| [ "psatz_Q" ] -> [ Coq_micromega.psatz_Q (-1) ] END + |