From 5b7eafd0f00a16d78f99a27f5c7d5a0de77dc7e6 Mon Sep 17 00:00:00 2001 From: Stephane Glondu Date: Wed, 21 Jul 2010 09:46:51 +0200 Subject: Imported Upstream snapshot 8.3~beta0+13298 --- plugins/micromega/csdpcert.ml | 214 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 plugins/micromega/csdpcert.ml (limited to 'plugins/micromega/csdpcert.ml') diff --git a/plugins/micromega/csdpcert.ml b/plugins/micromega/csdpcert.ml new file mode 100644 index 00000000..d4e6d920 --- /dev/null +++ b/plugins/micromega/csdpcert.ml @@ -0,0 +1,214 @@ +(************************************************************************) +(* v * The Coq Proof Assistant / The Coq Development Team *) +(* 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 + let p2 = expr_to_term p2 in + let res = Mul(p1,p2) in res + + | PEadd(p1,p2) -> Add(expr_to_term p1, expr_to_term p2) + | PEsub(p1,p2) -> Sub(expr_to_term p1, expr_to_term p2) + | PEpow(p,n) -> Pow(expr_to_term p , C2Ml.n n) + | PEopp p -> Opp (expr_to_term p) + + +end +open M + +open List +open Mutils + + + + +let rec canonical_sum_to_string = function s -> failwith "not implemented" + +let print_canonical_sum m = Format.print_string (canonical_sum_to_string m) + +let print_list_term o l = + output_string o "print_list_term\n"; + List.iter (fun (e,k) -> Printf.fprintf o "q: %s %s ;" + (string_of_poly (poly_of_term (expr_to_term e))) + (match k with + Mc.Equal -> "= " + | Mc.Strict -> "> " + | Mc.NonStrict -> ">= " + | _ -> failwith "not_implemented")) (List.map (fun (e, o) -> Mc.denorm e , o) l) ; + output_string o "\n" + + +let partition_expr l = + let rec f i = function + | [] -> ([],[],[]) + | (e,k)::l -> + let (eq,ge,neq) = f (i+1) l in + match k with + | Mc.Equal -> ((e,i)::eq,ge,neq) + | Mc.NonStrict -> (eq,(e,Axiom_le i)::ge,neq) + | Mc.Strict -> (* e > 0 == e >= 0 /\ e <> 0 *) + (eq, (e,Axiom_lt i)::ge,(e,Axiom_lt i)::neq) + | Mc.NonEqual -> (eq,ge,(e,Axiom_eq i)::neq) + (* Not quite sure -- Coq interface has changed *) + in f 0 l + + +let rec sets_of_list l = + match l with + | [] -> [[]] + | e::l -> let s = sets_of_list l in + s@(List.map (fun s0 -> e::s0) s) + +(* The exploration is probably not complete - for simple cases, it works... *) +let real_nonlinear_prover d l = + let l = List.map (fun (e,op) -> (Mc.denorm e,op)) l in + try + let (eq,ge,neq) = partition_expr l in + + let rec elim_const = function + [] -> [] + | (x,y)::l -> let p = poly_of_term (expr_to_term x) in + if poly_isconst p + then elim_const l + else (p,y)::(elim_const l) in + + let eq = elim_const eq in + let peq = List.map fst eq in + + let pge = List.map + (fun (e,psatz) -> poly_of_term (expr_to_term e),psatz) ge in + + let monoids = List.map (fun m -> (List.fold_right (fun (p,kd) y -> + let p = poly_of_term (expr_to_term p) in + match kd with + | Axiom_lt i -> poly_mul p y + | Axiom_eq i -> poly_mul (poly_pow p 2) y + | _ -> failwith "monoids") m (poly_const (Int 1)) , map snd m)) + (sets_of_list neq) in + + let (cert_ideal, cert_cone,monoid) = deepen_until d (fun d -> + list_try_find (fun m -> let (ci,cc) = + real_positivnullstellensatz_general false d peq pge (poly_neg (fst m) ) in + (ci,cc,snd m)) monoids) 0 in + + let proofs_ideal = map2 (fun q i -> Eqmul(term_of_poly q,Axiom_eq i)) + cert_ideal (List.map snd eq) in + + let proofs_cone = map term_of_sos cert_cone in + + let proof_ne = + let (neq , lt) = List.partition + (function Axiom_eq _ -> true | _ -> false ) monoid in + let sq = match + (List.map (function Axiom_eq i -> i | _ -> failwith "error") neq) + with + | [] -> Rational_lt (Int 1) + | l -> Monoid l in + List.fold_right (fun x y -> Product(x,y)) lt sq in + + let proof = list_fold_right_elements + (fun s t -> Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) in + S (Some proof) + with + | Sos_lib.TooDeep -> S None + | x -> F (Printexc.to_string x) + +(* This is somewhat buggy, over Z, strict inequality vanish... *) +let pure_sos l = + let l = List.map (fun (e,o) -> Mc.denorm e, o) l in + + (* If there is no strict inequality, + I should nonetheless be able to try something - over Z > is equivalent to -1 >= *) + try + let l = List.combine l (interval 0 (length l -1)) in + let (lt,i) = try (List.find (fun (x,_) -> snd x = Mc.Strict) l) + with Not_found -> List.hd l in + let plt = poly_neg (poly_of_term (expr_to_term (fst lt))) in + let (n,polys) = sumofsquares plt in (* n * (ci * pi^2) *) + let pos = Product (Rational_lt n, + List.fold_right (fun (c,p) rst -> Sum (Product (Rational_lt c, Square + (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 *) + S (Some proof) + with +(* | Sos.CsdpNotFound -> F "Sos.CsdpNotFound" *) + | x -> (* May be that could be refined *) S None + + + +let run_prover prover pb = + match prover with + | "real_nonlinear_prover", Some d -> real_nonlinear_prover d pb + | "pure_sos", None -> pure_sos pb + | prover, _ -> (Printf.printf "unknown prover: %s\n" prover; exit 1) + + +let output_csdp_certificate o = function + | S None -> output_string o "S None" + | S (Some p) -> Printf.fprintf o "S (Some %a)" output_psatz p + | F s -> Printf.fprintf o "F %s" s + + +let main () = + try + let (prover,poly) = (input_value stdin : provername * micromega_polys) in + let cert = run_prover prover poly in +(* Printf.fprintf chan "%a -> %a" print_list_term poly output_csdp_certificate cert ; + close_out chan ; *) + + output_value stdout (cert:csdp_certificate); + flush stdout ; + Marshal.to_channel chan (cert:csdp_certificate) [] ; + flush chan ; + exit 0 + with x -> (Printf.fprintf chan "error %s" (Printexc.to_string x) ; exit 1) + +;; + +let _ = main () in () + +(* Local Variables: *) +(* coding: utf-8 *) +(* End: *) -- cgit v1.2.3