From 19f49cb0088868cfb6c91d04924c44ffbf8d9fc7 Mon Sep 17 00:00:00 2001 From: Maxime Dénès Date: Thu, 3 May 2018 15:46:26 +0200 Subject: Micromega clean-up We add .mli files, removed dead code and use standard combinators instead of redefined ad-hoc ones in a few places. A lot of cleaning still has to be done on this code: documenting the interfaces, resolving the many abstraction leaks. I suspect there is still a lot of code duplication. --- plugins/micromega/certificate.ml | 194 +---------- plugins/micromega/certificate.mli | 22 ++ plugins/micromega/coq_micromega.ml | 284 +-------------- plugins/micromega/coq_micromega.mli | 22 ++ plugins/micromega/csdpcert.ml | 36 +- plugins/micromega/csdpcert.mli | 9 + plugins/micromega/g_micromega.mli | 9 + plugins/micromega/mfourier.ml | 85 +---- plugins/micromega/mfourier.mli | 49 +++ plugins/micromega/mutils.ml | 109 ------ plugins/micromega/mutils.mli | 70 ++++ plugins/micromega/persistent_cache.mli | 47 +++ plugins/micromega/polynomial.ml | 68 +--- plugins/micromega/polynomial.mli | 118 +++++++ plugins/micromega/sos.ml | 616 +++------------------------------ plugins/micromega/sos_lib.ml | 105 +----- plugins/micromega/sos_lib.mli | 79 +++++ 17 files changed, 516 insertions(+), 1406 deletions(-) create mode 100644 plugins/micromega/certificate.mli create mode 100644 plugins/micromega/coq_micromega.mli create mode 100644 plugins/micromega/csdpcert.mli create mode 100644 plugins/micromega/g_micromega.mli create mode 100644 plugins/micromega/mfourier.mli create mode 100644 plugins/micromega/mutils.mli create mode 100644 plugins/micromega/persistent_cache.mli create mode 100644 plugins/micromega/polynomial.mli create mode 100644 plugins/micromega/sos_lib.mli (limited to 'plugins/micromega') diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml index 9f39191f8..3a9709b6c 100644 --- a/plugins/micromega/certificate.ml +++ b/plugins/micromega/certificate.ml @@ -17,10 +17,9 @@ (* We take as input a list of polynomials [p1...pn] and return an unfeasibility certificate polynomial. *) -type var = int - - +let debug = false +open Util open Big_int open Num open Polynomial @@ -59,9 +58,6 @@ let q_spec = { eqb = Mc.qeq_bool } -let r_spec = z_spec - - let dev_form n_spec p = let rec dev_form p = match p with @@ -84,38 +80,6 @@ let dev_form n_spec p = pow n in dev_form p - -let monomial_to_polynomial mn = - Monomial.fold - (fun v i acc -> - let v = Ml2C.positive v in - let mn = if Int.equal i 1 then Mc.PEX v else Mc.PEpow (Mc.PEX v ,Ml2C.n i) in - if Pervasives.(=) acc (Mc.PEc (Mc.Zpos Mc.XH)) (** FIXME *) - then mn - else Mc.PEmul(mn,acc)) - mn - (Mc.PEc (Mc.Zpos Mc.XH)) - - - -let list_to_polynomial vars l = - assert (List.for_all (fun x -> ceiling_num x =/ x) l); - let var x = monomial_to_polynomial (List.nth vars x) in - - let rec xtopoly p i = function - | [] -> p - | c::l -> if c =/ (Int 0) then xtopoly p (i+1) l - else let c = Mc.PEc (Ml2C.bigint (numerator c)) in - let mn = - if Pervasives.(=) c (Mc.PEc (Mc.Zpos Mc.XH)) - then var i - else Mc.PEmul (c,var i) in - let p' = if Pervasives.(=) p (Mc.PEc Mc.Z0) then mn else - Mc.PEadd (mn, p) in - xtopoly p' (i+1) l in - - xtopoly (Mc.PEc Mc.Z0) 0 l - let rec fixpoint f x = let y' = f x in if Pervasives.(=) y' x then y' @@ -135,15 +99,6 @@ let rec_simpl_cone n_spec e = let simplify_cone n_spec c = fixpoint (rec_simpl_cone n_spec) c - -type cone_prod = - Const of cone -| Ideal of cone *cone -| Mult of cone * cone -| Other of cone -and cone = Mc.zWitness - - let factorise_linear_cone c = @@ -224,14 +179,6 @@ let positivity l = in xpositivity 0 l - -let string_of_op = function - | Mc.Strict -> "> 0" - | Mc.NonStrict -> ">= 0" - | Mc.Equal -> "= 0" - | Mc.NonEqual -> "<> 0" - - module MonSet = Set.Make(Monomial) (* If the certificate includes at least one strict inequality, @@ -261,9 +208,6 @@ let build_linear_system l = op = Ge ; cst = Big_int zero_big_int}::(strict::(positivity l)@s0) - -let big_int_to_z = Ml2C.bigint - (* 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) = @@ -296,8 +240,6 @@ let make_certificate n_spec (cert,li) = (simplify_cone n_spec (scalar_product cert' li))) -exception Found of Monomial.t - exception Strict module MonMap = Map.Make(Monomial) @@ -367,7 +309,7 @@ let simple_linear_prover l = let linear_prover n_spec l = let build_system n_spec l = - let li = List.combine l (interval 0 (List.length l -1)) in + 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 @@ -397,7 +339,7 @@ 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 = mapi (fun c i -> (c,Mc.PsatzIn (Ml2C.nat i))) sys in + let sys = List.mapi (fun i c -> (c,Mc.PsatzIn (Ml2C.nat i))) sys in (* Add all the product of hypotheses *) @@ -452,39 +394,6 @@ let nlinear_prover prfdepth (sys: (Mc.q Mc.pExpr * Mc.op1) list) = | Mc.PsatzZ -> Mc.PsatzZ in Some (map_psatz cert) - - -let make_linear_system l = - let l' = List.map fst l in - let monomials = List.fold_left (fun acc p -> Poly.addition p acc) - (Poly.constant (Int 0)) l' in - let monomials = Poly.fold - (fun mn _ l -> if Pervasives.(=) mn Monomial.const then l else mn::l) monomials [] in - (List.map (fun (c,op) -> - {coeffs = Vect.from_list (List.map (fun mn -> (Poly.get mn c)) monomials) ; - op = op ; - cst = minus_num ( (Poly.get Monomial.const c))}) l - ,monomials) - - -let pplus x y = Mc.PEadd(x,y) -let pmult x y = Mc.PEmul(x,y) -let pconst x = Mc.PEc x -let popp x = Mc.PEopp x - -(* keep track of enumerated vectors *) -let rec mem p x l = - match l with [] -> false | e::l -> if p x e then true else mem p x l - -let rec remove_assoc p x l = - match l with [] -> [] | e::l -> if p x (fst e) then - remove_assoc p x l else e::(remove_assoc p x l) - -let eq x y = Int.equal (Vect.compare x y) 0 - -let remove e l = List.fold_left (fun l x -> if eq x e then l else x::l) [] l - - (* The prover is (probably) incomplete -- only searching for naive cutting planes *) @@ -494,38 +403,6 @@ let develop_constraint z_spec (e,k) = | Mc.Equal -> (dev_form z_spec e , Eq) | _ -> assert false - -let op_of_op_compat = function - | Ge -> Mc.NonStrict - | Eq -> Mc.Equal - - -let integer_vector coeffs = - let vars , coeffs = List.split coeffs in - List.combine vars (List.map (fun x -> Big_int x) (rats_to_ints coeffs)) - -let integer_cstr {coeffs = coeffs ; op = op ; cst = cst } = - let vars , coeffs = List.split coeffs in - match rats_to_ints (cst::coeffs) with - | cst :: coeffs -> - { - coeffs = List.combine vars (List.map (fun x -> Big_int x) coeffs) ; - op = op ; cst = Big_int cst} - | _ -> assert false - - -let pexpr_of_cstr_compat var cstr = - let {coeffs = coeffs ; op = op ; cst = cst } = integer_cstr cstr in - try - let expr = list_to_polynomial var (Vect.to_list coeffs) in - let d = Ml2C.bigint (denominator cst) in - let n = Ml2C.bigint (numerator cst) in - (pplus (pmult (pconst d) expr) (popp (pconst n)), op_of_op_compat op) - with Failure _ -> failwith "pexpr_of_cstr_compat" - - - - open Sos_types let rec scale_term t = @@ -555,18 +432,6 @@ 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 Int.equal 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 @@ -681,8 +546,6 @@ open Polynomial module Env = struct - type t = int list - let id_of_hyp hyp l = let rec xid_of_hyp i l = match l with @@ -749,9 +612,6 @@ let xlinear_prover sys = | Inl _ -> None -let output_num o n = output_string o (string_of_num n) -let output_bigint o n = output_string o (string_of_big_int n) - 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 = @@ -894,23 +754,6 @@ let rec ext_gcd a b = let (s,t) = ext_gcd b r in (t, sub_big_int s (mult_big_int q t)) - -let pp_ext_gcd a b = - let a' = big_int_of_int a in - let b' = big_int_of_int b in - - let (x,y) = ext_gcd a' b' in - Printf.fprintf stdout "%s * %s + %s * %s = %s\n" - (string_of_big_int x) (string_of_big_int a') - (string_of_big_int y) (string_of_big_int b') - (string_of_big_int (add_big_int (mult_big_int x a') (mult_big_int y b'))) - -exception Result of (int * (proof * cstr_compat)) - -let split_equations psys = - List.partition (fun (c,p) -> c.op == Eq) - - let extract_coprime (c1,p1) (c2,p2) = let rec exist2 vect1 vect2 = match vect1 , vect2 with @@ -1058,29 +901,6 @@ let reduce_var_change psys = Some (apply_and_normalise pivot_eq sys) - - - -let reduce_pivot psys = - let is_equation (cstr,prf) = - if cstr.op == Eq - then - try - Some (fst (List.hd cstr.coeffs)) - with Not_found -> None - else None in - let (oeq,sys) = extract is_equation psys in - match oeq with - | None -> None (* Nothing to do *) - | Some(v,pc) -> - if debug then - Printf.printf "Bad news : loss of completeness %a=%s" Vect.pp_vect (fst pc).coeffs (string_of_num (fst pc).cst); - Some(pivot_sys v pc sys) - - - - - let iterate_until_stable f x = let rec iter x = match f x with @@ -1225,7 +1045,7 @@ let xlia (can_enum:bool) reduction_equations sys = | None -> None | Some prf -> (*Printf.printf "direct proof %a\n" output_proof prf ; *) - let env = mapi (fun _ i -> i) sys in + let env = List.mapi (fun i _ -> i) sys in let prf = compile_proof env prf in (*try if Mc.zChecker sys' prf then Some prf else @@ -1244,7 +1064,7 @@ let lia (can_enum:bool) (prfdepth:int) sys = 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 = mapi (fun c i -> (c,Hyp i)) sys in + let sys = List.mapi (fun i c -> (c,Hyp i)) sys in xlia can_enum reduction_equations sys @@ -1252,7 +1072,7 @@ 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 = mapi (fun c i -> (c,Hyp i)) 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 diff --git a/plugins/micromega/certificate.mli b/plugins/micromega/certificate.mli new file mode 100644 index 000000000..13d50d1ee --- /dev/null +++ b/plugins/micromega/certificate.mli @@ -0,0 +1,22 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* Mc.q Mc.psatz +val z_cert_of_pos : Sos_types.positivstellensatz -> Mc.z Mc.psatz +val lia : bool -> int -> (Mc.z Mc.pExpr * Mc.op1) list -> Mc.zArithProof option +val nlia : bool -> int -> (Mc.z Mc.pExpr * Mc.op1) list -> Mc.zArithProof option +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 c7abd58b0..68620dbfc 100644 --- a/plugins/micromega/coq_micromega.ml +++ b/plugins/micromega/coq_micromega.ml @@ -20,9 +20,9 @@ open Pp open Names -open Constr open Goptions open Mutils +open Constr (** * Debug flag @@ -30,19 +30,6 @@ open Mutils let debug = false -(** - * Time function - *) - -let time str f x = - let t0 = (Unix.times()).Unix.tms_utime in - let res = f x in - let t1 = (Unix.times()).Unix.tms_utime in - (*if debug then*) (Printf.printf "time %s %f\n" str (t1 -. t0) ; - flush stdout); - res - - (* Limit the proof search *) let max_depth = max_int @@ -305,8 +292,7 @@ let rec add_term t0 = function *) module ISet = Set.Make(Int) -module IMap = Map.Make(Int) - + (** * Given a set of integers s=\{i0,...,iN\} and a list m, return the list of * elements of m that are at position i0,...,iN. @@ -395,16 +381,10 @@ struct let coq_O = lazy (init_constant "O") let coq_S = lazy (init_constant "S") - let coq_nat = lazy (init_constant "nat") let coq_N0 = lazy (bin_constant "N0") let coq_Npos = lazy (bin_constant "Npos") - let coq_pair = lazy (init_constant "pair") - let coq_None = lazy (init_constant "None") - let coq_option = lazy (init_constant "option") - - let coq_positive = lazy (bin_constant "positive") let coq_xH = lazy (bin_constant "xH") let coq_xO = lazy (bin_constant "xO") let coq_xI = lazy (bin_constant "xI") @@ -417,8 +397,6 @@ struct let coq_Q = lazy (constant "Q") let coq_R = lazy (constant "R") - let coq_Build_Witness = lazy (constant "Build_Witness") - let coq_Qmake = lazy (constant "Qmake") let coq_Rcst = lazy (constant "Rcst") @@ -455,8 +433,6 @@ struct let coq_Zmult = lazy (z_constant "Z.mul") let coq_Zpower = lazy (z_constant "Z.pow") - 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") @@ -476,7 +452,6 @@ struct let coq_Rminus = lazy (r_constant "Rminus") let coq_Ropp = lazy (r_constant "Ropp") let coq_Rmult = lazy (r_constant "Rmult") - let coq_Rdiv = lazy (r_constant "Rdiv") let coq_Rinv = lazy (r_constant "Rinv") let coq_Rpower = lazy (r_constant "pow") let coq_IZR = lazy (r_constant "IZR") @@ -509,12 +484,6 @@ struct let coq_PsatzAdd = lazy (constant "PsatzAdd") let coq_PsatzC = lazy (constant "PsatzC") let coq_PsatzZ = lazy (constant "PsatzZ") - let coq_coneMember = lazy (constant "coneMember") - - let coq_make_impl = lazy - (gen_constant_in_modules "Zmicromega" [["Refl"]] "make_impl") - let coq_make_conj = lazy - (gen_constant_in_modules "Zmicromega" [["Refl"]] "make_conj") let coq_TT = lazy (gen_constant_in_modules "ZMicromega" @@ -552,13 +521,6 @@ struct let coq_QWitness = lazy (gen_constant_in_modules "QMicromega" [["Coq"; "micromega"; "QMicromega"]] "QWitness") - let coq_ZWitness = lazy - (gen_constant_in_modules "QMicromega" - [["Coq"; "micromega"; "ZMicromega"]] "ZWitness") - - let coq_N_of_Z = lazy - (gen_constant_in_modules "ZArithRing" - [["Coq";"setoid_ring";"ZArithRing"]] "N_of_Z") let coq_Build = lazy (gen_constant_in_modules "RingMicromega" @@ -577,24 +539,6 @@ struct * pp_* functions pretty-print Coq terms. *) - (* Error datastructures *) - - type parse_error = - | Ukn - | BadStr of string - | BadNum of int - | BadTerm of constr - | Msg of string - | Goal of (constr list ) * constr * parse_error - - let string_of_error = function - | Ukn -> "ukn" - | BadStr s -> s - | BadNum i -> string_of_int i - | BadTerm _ -> "BadTerm" - | Msg s -> s - | Goal _ -> "Goal" - exception ParseError (* A simple but useful getter function *) @@ -648,19 +592,6 @@ struct | Mc.N0 -> Lazy.force coq_N0 | Mc.Npos p -> EConstr.mkApp(Lazy.force coq_Npos,[| dump_positive p|]) - let rec dump_index x = - match x with - | Mc.XH -> Lazy.force coq_xH - | Mc.XO p -> EConstr.mkApp(Lazy.force coq_xO,[| dump_index p |]) - | Mc.XI p -> EConstr.mkApp(Lazy.force coq_xI,[| dump_index p |]) - - let pp_index o x = Printf.fprintf o "%i" (CoqToCaml.index x) - - let pp_n o x = output_string o (string_of_int (CoqToCaml.n x)) - - let dump_pair t1 t2 dump_t1 dump_t2 (x,y) = - EConstr.mkApp(Lazy.force coq_pair,[| t1 ; t2 ; dump_t1 x ; dump_t2 y|]) - let parse_z sigma term = let (i,c) = get_left_construct sigma term in match i with @@ -677,11 +608,6 @@ struct let pp_z o x = Printf.fprintf o "%s" (Big_int.string_of_big_int (CoqToCaml.z_big_int x)) - let dump_num bd1 = - EConstr.mkApp(Lazy.force coq_Qmake, - [|dump_z (CamlToCoq.bigint (numerator bd1)) ; - dump_positive (CamlToCoq.positive_big_int (denominator bd1)) |]) - let dump_q q = EConstr.mkApp(Lazy.force coq_Qmake, [| dump_z q.Micromega.qnum ; dump_positive q.Micromega.qden|]) @@ -719,29 +645,6 @@ struct | Mc.CInv t -> EConstr.mkApp(Lazy.force coq_CInv, [| dump_Rcst t |]) | Mc.COpp t -> EConstr.mkApp(Lazy.force coq_COpp, [| dump_Rcst t |]) - let rec parse_Rcst sigma term = - let (i,c) = get_left_construct sigma term in - match i with - | 1 -> Mc.C0 - | 2 -> Mc.C1 - | 3 -> Mc.CQ (parse_q sigma c.(0)) - | 4 -> Mc.CPlus(parse_Rcst sigma c.(0), parse_Rcst sigma c.(1)) - | 5 -> Mc.CMinus(parse_Rcst sigma c.(0), parse_Rcst sigma c.(1)) - | 6 -> Mc.CMult(parse_Rcst sigma c.(0), parse_Rcst sigma c.(1)) - | 7 -> Mc.CInv(parse_Rcst sigma c.(0)) - | 8 -> Mc.COpp(parse_Rcst sigma c.(0)) - | _ -> raise ParseError - - - - - let rec parse_list sigma parse_elt term = - let (i,c) = get_left_construct sigma term in - match i with - | 1 -> [] - | 2 -> parse_elt sigma c.(1) :: parse_list sigma parse_elt c.(2) - | i -> raise ParseError - let rec dump_list typ dump_elt l = match l with | [] -> EConstr.mkApp(Lazy.force coq_nil,[| typ |]) @@ -756,22 +659,8 @@ struct | e::l -> Printf.fprintf o "%a ,%a" elt e _pp l in Printf.fprintf o "%s%a%s" op _pp l cl - let pp_var = pp_positive - let dump_var = dump_positive - let pp_expr pp_z o e = - let rec pp_expr o e = - match e with - | Mc.PEX n -> Printf.fprintf o "V %a" pp_var n - | Mc.PEc z -> pp_z o z - | Mc.PEadd(e1,e2) -> Printf.fprintf o "(%a)+(%a)" pp_expr e1 pp_expr e2 - | Mc.PEmul(e1,e2) -> Printf.fprintf o "%a*(%a)" pp_expr e1 pp_expr e2 - | Mc.PEopp e -> Printf.fprintf o "-(%a)" pp_expr e - | Mc.PEsub(e1,e2) -> Printf.fprintf o "(%a)-(%a)" pp_expr e1 pp_expr e2 - | Mc.PEpow(e,n) -> Printf.fprintf o "(%a)^(%a)" pp_expr e pp_n n in - pp_expr o e - let dump_expr typ dump_z e = let rec dump_expr e = match e with @@ -854,18 +743,6 @@ struct | Mc.OpGt-> Lazy.force coq_OpGt | Mc.OpLt-> Lazy.force coq_OpLt - let pp_op o e= - match e with - | Mc.OpEq-> Printf.fprintf o "=" - | Mc.OpNEq-> Printf.fprintf o "<>" - | Mc.OpLe -> Printf.fprintf o "=<" - | Mc.OpGe -> Printf.fprintf o ">=" - | Mc.OpGt-> Printf.fprintf o ">" - | Mc.OpLt-> Printf.fprintf o "<" - - let pp_cstr pp_z o {Mc.flhs = l ; Mc.fop = op ; Mc.frhs = r } = - Printf.fprintf o"(%a %a %a)" (pp_expr pp_z) l pp_op op (pp_expr pp_z) r - let dump_cstr typ dump_constant {Mc.flhs = e1 ; Mc.fop = o ; Mc.frhs = e2} = EConstr.mkApp(Lazy.force coq_Build, [| typ; dump_expr typ dump_constant e1 ; @@ -924,11 +801,6 @@ struct let parse_qop gl (op,args) = (assoc_const gl.sigma op qop_table, args.(0) , args.(1)) - let is_constant sigma t = (* This is an approx *) - match EConstr.kind sigma t with - | Construct(i,_) -> true - | _ -> false - type 'a op = | Binop of ('a Mc.pExpr -> 'a Mc.pExpr -> 'a Mc.pExpr) | Opp @@ -947,8 +819,6 @@ struct module Env = struct - type t = EConstr.constr list - let compute_rank_add env sigma v = let rec _add env n v = match env with @@ -1168,17 +1038,6 @@ struct (* generic parsing of arithmetic expressions *) - let rec f2f = function - | TT -> Mc.TT - | FF -> Mc.FF - | X _ -> Mc.X - | A (x,_,_) -> Mc.A x - | C (a,b) -> Mc.Cj(f2f a,f2f b) - | D (a,b) -> Mc.D(f2f a,f2f b) - | N (a) -> Mc.N(f2f a) - | I(a,_,b) -> Mc.I(f2f a,f2f b) - - let mkC f1 f2 = C(f1,f2) let mkD f1 f2 = D(f1,f2) let mkIff f1 f2 = C(I(f1,None,f2),I(f2,None,f1)) @@ -1323,31 +1182,6 @@ let dump_qexpr = lazy dump_op = List.map (fun (x,y) -> (y,Lazy.force x)) qop_table } - let dump_positive_as_R p = - let mult = Lazy.force coq_Rmult in - let add = Lazy.force coq_Rplus in - - let one = Lazy.force coq_R1 in - let mk_add x y = EConstr.mkApp(add,[|x;y|]) in - let mk_mult x y = EConstr.mkApp(mult,[|x;y|]) in - - let two = mk_add one one in - - let rec dump_positive p = - match p with - | Mc.XH -> one - | Mc.XO p -> mk_mult two (dump_positive p) - | Mc.XI p -> mk_add one (mk_mult two (dump_positive p)) in - - dump_positive p - -let dump_n_as_R n = - let z = CoqToCaml.n n in - if z = 0 - then Lazy.force coq_R0 - else dump_positive_as_R (CamlToCoq.positive z) - - let rec dump_Rcst_as_R cst = match cst with | Mc.C0 -> Lazy.force coq_R0 @@ -1481,54 +1315,6 @@ end (** open M -let rec sig_of_cone = function - | Mc.PsatzIn n -> [CoqToCaml.nat n] - | Mc.PsatzMulE(w1,w2) -> (sig_of_cone w1)@(sig_of_cone w2) - | Mc.PsatzMulC(w1,w2) -> (sig_of_cone w2) - | Mc.PsatzAdd(w1,w2) -> (sig_of_cone w1)@(sig_of_cone w2) - | _ -> [] - -let same_proof sg cl1 cl2 = - let rec xsame_proof sg = - match sg with - | [] -> true - | n::sg -> - (try Int.equal (List.nth cl1 n) (List.nth cl2 n) with Invalid_argument _ -> false) - && (xsame_proof sg ) in - xsame_proof sg - -let tags_of_clause tgs wit clause = - let rec xtags tgs = function - | Mc.PsatzIn n -> Names.Id.Set.union tgs - (snd (List.nth clause (CoqToCaml.nat n) )) - | Mc.PsatzMulC(e,w) -> xtags tgs w - | Mc.PsatzMulE (w1,w2) | Mc.PsatzAdd(w1,w2) -> xtags (xtags tgs w1) w2 - | _ -> tgs in - xtags tgs wit - -(*let tags_of_cnf wits cnf = - List.fold_left2 (fun acc w cl -> tags_of_clause acc w cl) - Names.Id.Set.empty wits cnf *) - -let find_witness prover polys1 = try_any prover polys1 - -let rec witness prover l1 l2 = - match l2 with - | [] -> Some [] - | e :: l2 -> - match find_witness prover (e::l1) with - | None -> None - | Some w -> - (match witness prover l1 l2 with - | None -> None - | Some l -> Some (w::l) - ) - -let rec apply_ids t ids = - match ids with - | [] -> t - | i::ids -> apply_ids (mkApp(t,[| mkVar i |])) ids - let coq_Node = lazy (gen_constant_in_modules "VarMap" [["Coq" ; "micromega" ; "VarMap"];["VarMap"]] "Node") @@ -1559,15 +1345,6 @@ let vm_of_list env = List.fold_left (fun vm (c,i) -> Mc.vm_add d (CamlToCoq.positive i) c vm) Mc.Empty env - -let rec pp_varmap o vm = - match vm with - | Mc.Empty -> output_string o "[]" - | Mc.Leaf z -> Printf.fprintf o "[%a]" pp_z z - | Mc.Node(l,z,r) -> Printf.fprintf o "[%a, %a, %a]" pp_varmap l pp_z z pp_varmap r - - - let rec dump_proof_term = function | Micromega.DoneProof -> Lazy.force coq_doneProof | Micromega.RatProof(cone,rst) -> @@ -1662,45 +1439,11 @@ let qq_domain_spec = lazy { dump_proof = dump_psatz coq_Q dump_q } -let rcst_domain_spec = lazy { - typ = Lazy.force coq_R; - coeff = Lazy.force coq_Rcst; - dump_coeff = dump_Rcst; - proof_typ = Lazy.force coq_QWitness ; - dump_proof = dump_psatz coq_Q dump_q -} - (** Naive topological sort of constr according to the subterm-ordering *) (* An element is minimal x is minimal w.r.t y if x <= y or (x and y are incomparable) *) -let is_min le x y = - if le x y then true - else if le y x then false else true - -let is_minimal le l c = List.for_all (is_min le c) l - -let find_rem p l = - let rec xfind_rem acc l = - match l with - | [] -> (None, acc) - | x :: l -> if p x then (Some x, acc @ l) - else xfind_rem (x::acc) l in - xfind_rem [] l - -let find_minimal le l = find_rem (is_minimal le l) l - -let rec mk_topo_order le l = - match find_minimal le l with - | (None , _) -> [] - | (Some v,l') -> v :: (mk_topo_order le l') - - -let topo_sort_constr l = - mk_topo_order (fun c t -> Termops.dependent Evd.empty (** FIXME *) (EConstr.of_constr c) (EConstr.of_constr t)) l - - (** * Instanciate the current Coq goal with a Micromega formula, a varmap, and a * witness. @@ -1778,13 +1521,6 @@ let witness_list prover l = let witness_list_tags = witness_list -(* *Deprecated* let is_singleton = function [] -> true | [e] -> true | _ -> false *) - -let pp_ml_list pp_elt o l = - output_string o "[" ; - List.iter (fun x -> Printf.fprintf o "%a ;" pp_elt x) l ; - output_string o "]" - (** * Prune the proof object, according to the 'diff' between two cnf formulas. *) @@ -1792,7 +1528,7 @@ let pp_ml_list pp_elt o l = let compact_proofs (cnf_ff: 'cst cnf) res (cnf_ff': 'cst cnf) = let compact_proof (old_cl:'cst clause) (prf,prover) (new_cl:'cst clause) = - let new_cl = Mutils.mapi (fun (f,_) i -> (f,i)) new_cl in + let new_cl = List.mapi (fun i (f,_) -> (f,i)) new_cl in let remap i = let formula = try fst (List.nth old_cl i) with Failure _ -> failwith "bad old index" in List.assoc formula new_cl in @@ -2158,7 +1894,11 @@ let lift_ratproof prover l = | Some c -> Some (Mc.RatProof( c,Mc.DoneProof)) type micromega_polys = (Micromega.q Mc.pol * Mc.op1) list + +[@@@ocaml.warning "-37"] type csdp_certificate = S of Sos_types.positivstellensatz option | F of string +(* Used to read the result of the execution of csdpcert *) + type provername = string * int option (** @@ -2406,16 +2146,6 @@ let nlinear_Z = { pp_f = fun o x -> pp_pol pp_z o (fst x) } - - -let tauto_lia ff = - let prover = linear_Z in - let cnf_ff,_ = cnf Mc.negate Mc.normalise Mc.zunsat Mc.zdeduce ff in - match witness_list_tags [prover] cnf_ff with - | None -> None - | Some l -> Some (List.map fst l) - - (** * Functions instantiating micromega_gen with the appropriate theories and * solvers diff --git a/plugins/micromega/coq_micromega.mli b/plugins/micromega/coq_micromega.mli new file mode 100644 index 000000000..b91feb398 --- /dev/null +++ b/plugins/micromega/coq_micromega.mli @@ -0,0 +1,22 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* unit Proofview.tactic -> unit Proofview.tactic +val psatz_Q : int -> unit Proofview.tactic -> unit Proofview.tactic +val psatz_R : int -> unit Proofview.tactic -> unit Proofview.tactic +val xlia : unit Proofview.tactic -> unit Proofview.tactic +val xnlia : unit Proofview.tactic -> unit Proofview.tactic +val nra : unit Proofview.tactic -> unit Proofview.tactic +val nqa : unit Proofview.tactic -> unit Proofview.tactic +val sos_Z : unit Proofview.tactic -> unit Proofview.tactic +val sos_Q : unit Proofview.tactic -> unit Proofview.tactic +val sos_R : unit Proofview.tactic -> unit Proofview.tactic +val lra_Q : unit Proofview.tactic -> unit Proofview.tactic +val lra_R : unit Proofview.tactic -> unit Proofview.tactic diff --git a/plugins/micromega/csdpcert.ml b/plugins/micromega/csdpcert.ml index a1245b7cc..9c1b4810d 100644 --- a/plugins/micromega/csdpcert.ml +++ b/plugins/micromega/csdpcert.ml @@ -20,7 +20,6 @@ open Sos_types open Sos_lib module Mc = Micromega -module Ml2C = Mutils.CamlToCoq module C2Ml = Mutils.CoqToCaml type micromega_polys = (Micromega.q Mc.pol * Mc.op1) list @@ -28,7 +27,6 @@ type csdp_certificate = S of Sos_types.positivstellensatz option | F of string type provername = string * int option -let debug = false let flags = [Open_append;Open_binary;Open_creat] let chan = open_out_gen flags 0o666 "trace" @@ -55,27 +53,6 @@ struct end open M -open Mutils - - - - -let 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 | [] -> ([],[],[]) @@ -125,7 +102,7 @@ let real_nonlinear_prover d l = (sets_of_list neq) in let (cert_ideal, cert_cone,monoid) = deepen_until d (fun d -> - list_try_find (fun m -> let (ci,cc) = + tryfind (fun m -> let (ci,cc) = real_positivnullstellensatz_general false d peq pge (poly_neg (fst m) ) in (ci,cc,snd m)) monoids) 0 in @@ -144,7 +121,7 @@ let real_nonlinear_prover d l = | l -> Monoid l in List.fold_right (fun x y -> Product(x,y)) lt sq in - let proof = list_fold_right_elements + let proof = end_itlist (fun s t -> Sum(s,t)) (proof_ne :: proofs_ideal @ proofs_cone) in S (Some proof) with @@ -158,7 +135,7 @@ let pure_sos l = (* 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 (List.length l -1)) in + let l = List.combine l (CList.interval 0 (List.length l -1)) in let (lt,i) = try (List.find (fun (x,_) -> Pervasives.(=) (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 @@ -183,13 +160,6 @@ let run_prover prover 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 diff --git a/plugins/micromega/csdpcert.mli b/plugins/micromega/csdpcert.mli new file mode 100644 index 000000000..7c3ee6004 --- /dev/null +++ b/plugins/micromega/csdpcert.mli @@ -0,0 +1,9 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* if a <=/ b then Some itv else None | _ -> Some itv - (** [opp_itv itv] computes the opposite interval *) - let opp_itv itv = - let (l,r) = itv in - (map_option minus_num r, map_option minus_num l) - - - - (** [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 = @@ -92,10 +77,6 @@ type vector = Vect.t module ISet = Set.Make(Int) - -module PSet = ISet - - module System = Hashtbl.Make(Vect) type proof = @@ -131,14 +112,6 @@ and cstr_info = { (** To be thrown when a system has no solution *) exception SystemContradiction of proof -let hyps prf = - let rec hyps prf acc = - match prf with - | Assum i -> ISet.add i acc - | Elim(_,prf1,prf2) - | And(prf1,prf2) -> hyps prf1 (hyps prf2 acc) in - hyps prf ISet.empty - (** Pretty printing *) let rec pp_proof o prf = @@ -147,26 +120,6 @@ let hyps prf = | Elim(v, prf1,prf2) -> Printf.fprintf o "E(%i,%a,%a)" v pp_proof prf1 pp_proof prf2 | And(prf1,prf2) -> Printf.fprintf o "A(%a,%a)" pp_proof prf1 pp_proof prf2 -let pp_bound o = function - | None -> output_string o "oo" - | Some a -> output_string o (string_of_num a) - -let pp_itv o (l,r) = Printf.fprintf o "(%a,%a)" pp_bound l pp_bound r - - -let pp_iset o s = - output_string o "{" ; - ISet.fold (fun i _ -> Printf.fprintf o "%i " i) s (); - output_string o "}" - -let pp_pset o s = - output_string o "{" ; - PSet.fold (fun i _ -> Printf.fprintf o "%i " i) s (); - output_string o "}" - - -let pp_info o i = pp_itv o i.bound - let pp_cstr o (vect,bnd) = let (l,r) = bnd in (match l with @@ -183,11 +136,6 @@ let pp_system o sys= System.iter (fun vect ibnd -> pp_cstr o (vect,(!ibnd).bound)) sys - - -let pp_split_cstr o (vl,v,c,_) = - Printf.fprintf o "(val x = %s ,%a,%s)" (string_of_num vl) pp_vect v (string_of_num c) - (** [merge_cstr_info] takes: - the intersection of bounds and - the union of proofs @@ -243,8 +191,8 @@ let normalise_cstr vect cinfo = (if n <>/ Int 1 then List.map (fun (x,nx) -> (x,nx // n)) vect else vect), let divn x = x // n in if Int.equal (sign_num n) 1 - then{cinfo with bound = (map_option divn l , map_option divn r) } - else {cinfo with pos = cinfo.neg ; neg = cinfo.pos ; bound = (map_option divn r , map_option divn l)}) + 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 *) @@ -281,7 +229,7 @@ let load_system l = let sys = System.create 1000 in - let li = Mutils.mapi (fun e i -> (e,i)) l in + let li = List.mapi (fun i e -> (e,i)) l in let vars = List.fold_left (fun vrs (cstr,i) -> match norm_cstr cstr i with @@ -335,9 +283,6 @@ let add (v1,c1) (v2,c2) = (* Printf.printf "add(%a,%s,%a,%s) -> %a\n" pp_vect v1 (string_of_num c1) pp_vect v2 (string_of_num c2) pp_vect (fst res) ;*) res -type tlr = (num * vector * cstr_info) list -type tm = (vector * cstr_info ) list - (** To perform Fourier elimination, constraints are categorised depending on the sign of the variable to eliminate. *) (** [split x vect info (l,m,r)] @@ -381,8 +326,8 @@ let project vr sys = let {neg = n1 ; pos = p1 ; bound = bound1 ; prf = prf1} = info1 and {neg = n2 ; pos = p2 ; bound = bound2 ; prf = prf2} = info2 in - let bnd1 = from_option (fst bound1) - and bnd2 = from_option (fst bound2) in + let bnd1 = Option.get (fst bound1) + and bnd2 = Option.get (fst bound2) in let bound = (bnd1 // v1) +/ (bnd2 // minus_num v2) in let vres,(n,p) = add (vect1,v1) (vect2,minus_num v2) in (vres,{neg = n ; pos = p ; bound = (Some bound, None); prf = Elim(vr,info1.prf,info2.prf)}) in @@ -419,13 +364,13 @@ let project_using_eq vr c vect bound prf (vect',info') = let bndres = let f x = cst +/ x // c2 in let (l,r) = info'.bound in - (map_option f l , map_option f r) in + (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 = from_option (get vr vect) in + let c = Option.get (get vr vect) in let elim_var = project_using_eq vr c vect cst prf in @@ -444,9 +389,7 @@ let elim_var_using_eq vr vect cst prf sys = (** [size sys] computes the number of entries in the system of constraints *) let size sys = System.fold (fun v iref s -> s + (!iref).neg + (!iref).pos) sys 0 -module IMap = Map.Make(Int) - -let pp_map o map = IMap.fold (fun k elt () -> Printf.fprintf o "%i -> %s\n" k (string_of_num elt)) map () +module IMap = CMap.Make(Int) (** [eval_vect map vect] evaluates vector [vect] using the values of [map]. If [map] binds all the variables of [vect], we get @@ -475,8 +418,8 @@ let restrict_bound n sum (itv:interval) = | 0 -> if in_bound itv sum then (None,None) (* redundant *) else failwith "SystemContradiction" - | 1 -> map_option f l , map_option f r - | _ -> map_option f r , map_option f l + | 1 -> Option.map f l , Option.map f r + | _ -> Option.map f r , Option.map f l (** [bound_of_variable map v sys] computes the interval of [v] in @@ -613,12 +556,6 @@ struct |(Some a, Some b) -> a =/ b | _ -> false - let eq_bound bnd c = - match bnd with - |(Some a, Some b) -> a =/ b && c =/ b - | _ -> false - - let rec unroll_until v l = match l with | [] -> (false,[]) diff --git a/plugins/micromega/mfourier.mli b/plugins/micromega/mfourier.mli new file mode 100644 index 000000000..f1d8edeab --- /dev/null +++ b/plugins/micromega/mfourier.mli @@ -0,0 +1,49 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* 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 -> + Itv.interval option + +end + +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 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/mutils.ml b/plugins/micromega/mutils.ml index 82367c0b2..9d03560b7 100644 --- a/plugins/micromega/mutils.ml +++ b/plugins/micromega/mutils.ml @@ -19,8 +19,6 @@ (* *) (************************************************************************) -let debug = false - let rec pp_list f o l = match l with | [] -> () @@ -36,15 +34,6 @@ let finally f rst = with any -> raise reraise ); raise reraise -let map_option f x = - match x with - | None -> None - | Some v -> Some (f v) - -let from_option = function - | None -> failwith "from_option" - | Some v -> v - let rec try_any l x = match l with | [] -> None @@ -52,13 +41,6 @@ let rec try_any l x = | None -> try_any l x | x -> x -let iteri f l = - let rec xiter i l = - match l with - | [] -> () - | e::l -> f i e ; xiter (i+1) l in - xiter 0 l - let all_sym_pairs f l = let pair_with acc e l = List.fold_left (fun acc x -> (f e x) ::acc) acc l in @@ -77,14 +59,6 @@ let all_pairs f l = | e::lx -> xpairs (pair_with acc e l) lx in xpairs [] l - - -let rec map3 f l1 l2 l3 = - match l1 , l2 ,l3 with - | [] , [] , [] -> [] - | e1::l1 , e2::l2 , e3::l3 -> (f e1 e2 e3)::(map3 f l1 l2 l3) - | _ -> invalid_arg "map3" - let rec is_sublist f l1 l2 = match l1 ,l2 with | [] ,_ -> true @@ -93,26 +67,6 @@ let rec is_sublist f l1 l2 = if f e e' then is_sublist f l1' l2' else is_sublist f l1 l2' -let list_try_find f = - let rec try_find_f = function - | [] -> failwith "try_find" - | h::t -> try f h with Failure _ -> try_find_f t - in - try_find_f - -let list_fold_right_elements f l = - let rec aux = function - | [] -> invalid_arg "list_fold_right_elements" - | [x] -> x - | x::l -> f x (aux l) in - aux l - -let interval n m = - let rec interval_n (l,m) = - if n > m then l else interval_n (m::l,pred m) - in - interval_n ([],m) - let extract pred l = List.fold_left (fun (fd,sys) e -> match fd with @@ -163,51 +117,7 @@ let rats_to_ints l = List.map (fun x -> (div_big_int (mult_big_int (numerator x) c) (denominator x))) l -(* Nasty reordering of lists - useful to trim certificate down *) -let mapi f l = - let rec xmapi i l = - match l with - | [] -> [] - | e::l -> (f e i)::(xmapi (i+1) l) in - xmapi 0 l - -let concatMapi f l = List.rev (mapi (fun e i -> (i,f e)) l) - (* assoc_pos j [a0...an] = [j,a0....an,j+n],j+n+1 *) -let assoc_pos j l = (mapi (fun e i -> e,i+j) l, j + (List.length l)) - -let assoc_pos_assoc l = - let rec xpos i l = - match l with - | [] -> [] - | (x,l) ::rst -> let (l',j) = assoc_pos i l in - (x,l')::(xpos j rst) in - xpos 0 l - -let filter_pos f l = - (* Could sort ... take care of duplicates... *) - let rec xfilter l = - match l with - | [] -> [] - | (x,e)::l -> - if List.exists (fun ee -> List.mem ee f) (List.map snd e) - then (x,e)::(xfilter l) - else xfilter l in - xfilter l - -let select_pos lpos l = - let rec xselect i lpos l = - match lpos with - | [] -> [] - | j::rpos -> - match l with - | [] -> failwith "select_pos" - | e::l -> - if Int.equal i j - then e:: (xselect (i+1) rpos l) - else xselect (i+1) lpos l in - xselect 0 lpos l - (** * MODULE: Coq to Caml data-structure mappings *) @@ -238,12 +148,6 @@ struct | XI i -> 1+(2*(index i)) | XO i -> 2*(index i) - let z x = - match x with - | Z0 -> 0 - | Zpos p -> (positive p) - | Zneg p -> - (positive p) - open Big_int let rec positive_big_int p = @@ -258,8 +162,6 @@ struct | Zpos p -> (positive_big_int p) | Zneg p -> minus_big_int (positive_big_int p) - let num x = Num.Big_int (z_big_int x) - let q_to_num {qnum = x ; qden = y} = Big_int (z_big_int x) // (Big_int (z_big_int (Zpos y))) @@ -352,17 +254,6 @@ struct let c = cmp e1 e2 in if Int.equal c 0 then compare_list cmp l1 l2 else c -(** - * hash_list takes a hash function and a list, and computes an integer which - * is the hash value of the list. - *) - let hash_list hash l = - let rec _hash_list l h = - match l with - | [] -> h lxor (Hashtbl.hash []) - | e::l -> _hash_list l ((hash e) lxor h) - in _hash_list l 0 - end (** diff --git a/plugins/micromega/mutils.mli b/plugins/micromega/mutils.mli new file mode 100644 index 000000000..7b7a090de --- /dev/null +++ b/plugins/micromega/mutils.mli @@ -0,0 +1,70 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* Big_int.big_int +val denominator : Num.num -> Big_int.big_int + +module Cmp : sig + + val compare_list : ('a -> 'b -> int) -> 'a list -> 'b list -> int + val compare_lexical : (unit -> int) list -> int + +end + +module Tag : sig + + type t + + val pp : out_channel -> t -> unit + val next : t -> t + val from : int -> t + +end + +module TagSet : CSig.SetS with type elt = Tag.t + +val pp_list : (out_channel -> 'a -> 'b) -> out_channel -> 'a list -> unit + +module CamlToCoq : sig + + val positive : int -> Micromega.positive + val bigint : Big_int.big_int -> Micromega.z + val n : int -> Micromega.n + val nat : int -> Micromega.nat + val q : Num.num -> Micromega.q + val index : int -> Micromega.positive + val z : int -> Micromega.z + val positive_big_int : Big_int.big_int -> Micromega.positive + +end + +module CoqToCaml : sig + + val z_big_int : Micromega.z -> Big_int.big_int + val q_to_num : Micromega.q -> Num.num + val positive : Micromega.positive -> int + val n : Micromega.n -> int + val nat : Micromega.nat -> int + val index : Micromega.positive -> int + +end + +val rats_to_ints : Num.num list -> Big_int.big_int list + +val all_pairs : ('a -> 'a -> 'b) -> 'a list -> 'b list +val all_sym_pairs : ('a -> 'a -> 'b) -> 'a list -> 'b list +val try_any : (('a -> 'b option) * 'c) list -> 'a -> 'b option +val is_sublist : ('a -> 'b -> bool) -> 'a list -> 'b list -> bool + +val gcd_list : Num.num list -> Big_int.big_int + +val extract : ('a -> 'b option) -> 'a list -> ('b * 'a) option * 'a list + +val command : string -> string array -> 'a -> 'b diff --git a/plugins/micromega/persistent_cache.mli b/plugins/micromega/persistent_cache.mli new file mode 100644 index 000000000..240fa490f --- /dev/null +++ b/plugins/micromega/persistent_cache.mli @@ -0,0 +1,47 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* string -> 'a t + (** [create i f] creates an empty persistent table + with initial size i associated with file [f] *) + + + val open_in : string -> 'a t + (** [open_in f] rebuilds a table from the records stored in file [f]. + As marshaling is not type-safe, it migth segault. + *) + + val find : 'a t -> key -> 'a + (** find has the specification of Hashtable.find *) + + val add : 'a t -> key -> 'a -> unit + (** [add tbl key elem] adds the binding [key] [elem] to the table [tbl]. + (and writes the binding to the file associated with [tbl].) + If [key] is already bound, raises KeyAlreadyBound *) + + val close : 'a t -> unit + (** [close tbl] is closing the table. + Once closed, a table cannot be used. + i.e, find,add will raise UnboundTable *) + + val memo : string -> (key -> 'a) -> (key -> 'a) + (** [memo cache f] returns a memo function for [f] using file [cache] as persistent table. + Note that the cache will only be loaded when the function is used for the first time *) + + end + +module PHashtable(Key:HashedType) : PHashtable with type key = Key.t diff --git a/plugins/micromega/polynomial.ml b/plugins/micromega/polynomial.ml index db8b73a20..1d18a26f3 100644 --- a/plugins/micromega/polynomial.ml +++ b/plugins/micromega/polynomial.ml @@ -20,9 +20,9 @@ open Utils type var = int +let debug = false let (<+>) = add_num -let (<->) = minus_num let (<*>) = mult_num @@ -33,8 +33,6 @@ sig val is_const : t -> bool val var : var -> t val is_var : t -> bool - val find : var -> t -> int - val mult : var -> t -> t val prod : t -> t -> t val exp : t -> int -> t val div : t -> t -> t * int @@ -99,9 +97,6 @@ struct (* Get the degre of a variable in a monomial *) let find x m = try find x m with Not_found -> 0 - (* Multiply a monomial by a variable *) - let mult x m = add x ( (find x m) + 1) m - (* Product of monomials *) let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2 @@ -145,14 +140,10 @@ sig val variable : var -> t val add : Monomial.t -> num -> t -> t val constant : num -> t - val mult : Monomial.t -> num -> t -> 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 pp : out_channel -> t -> unit - val compare : t -> t -> int - val is_null : t -> bool val is_linear : t -> bool end = struct @@ -162,12 +153,6 @@ struct type t = num P.t - let pp o p = P.iter - (fun k v -> - if Monomial.compare Monomial.const k = 0 - 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 *) let get : Monomial.t -> t -> num = fun mn p -> try find mn p with Not_found -> (Int 0) @@ -220,10 +205,6 @@ struct let fold = P.fold - let is_null p = fold (fun mn vl b -> b && sign_num vl = 0) p true - - let compare = compare compare_num - let is_linear p = P.fold (fun m _ acc -> acc && (Monomial.is_const m || Monomial.is_var m)) p true (* let is_linear p = @@ -277,7 +258,6 @@ module Vect = xfrom_list 0 l let zero_num = Int 0 - let unit_num = Int 1 let to_list m = @@ -311,11 +291,6 @@ module Vect = | 1 -> (k,v) :: (set i n l) | _ -> failwith "compare_num" - let gcd m = - let res = List.fold_left (fun x (i,e) -> Big_int.gcd_big_int x (Utils.numerator e)) Big_int.zero_big_int m in - if Big_int.compare_big_int res Big_int.zero_big_int = 0 - then Big_int.unit_big_int else res - let mul z t = match z with | Int 0 -> [] @@ -345,7 +320,7 @@ module Vect = - let compare : t -> t -> int = Utils.Cmp.compare_list (fun x y -> Utils.Cmp.compare_lexical + 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))]) @@ -395,18 +370,8 @@ let opMult o1 o2 = | Eq , Ge | Ge , Eq -> Ge | Ge , Ge -> Ge -let opAdd o1 o2 = - match o1 , o2 with - | Eq , _ | _ , Eq -> Eq - | Ge , Ge -> Ge - - - - open Big_int -type index = int - type prf_rule = | Hyp of int | Def of int @@ -550,35 +515,6 @@ let mul_proof_ext (p,c) prf = | _ -> MulC((p,c),prf) - -(* - let rec scale_prf_rule = function - | Hyp i -> (unit_big_int, Hyp i) - | Def i -> (unit_big_int, Def i) - | Cst c -> (unit_big_int, Cst i) - | Zero -> (unit_big_int, Zero) - | Square p -> (unit_big_int,Square p) - | Div(c,pr) -> - let (bi,pr') = scale_prf_rule pr in - (mult_big_int c bi , pr') - | MulC(p,pr) -> - let bi,pr' = scale_prf_rule pr in - (bi,MulC p,pr') - | MulPrf(p1,p2) -> - let b1,p1 = scale_prf_rule p1 in - let b2,p2 = scale_prf_rule p2 in - - - | AddPrf(p1,p2) -> - let b1,p1 = scale_prf_rule p1 in - let b2,p2 = scale_prf_rule p2 in - let g = gcd_big_int -*) - - - - - module LinPoly = struct type t = Vect.t * num diff --git a/plugins/micromega/polynomial.mli b/plugins/micromega/polynomial.mli new file mode 100644 index 000000000..4c095202a --- /dev/null +++ b/plugins/micromega/polynomial.mli @@ -0,0 +1,118 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* int -> 'a -> 'a) -> t -> 'a -> 'a + val const : t + val sqrt : t -> t option + val is_var : t -> bool + val div : t -> t -> t * int + + val compare : t -> t -> int + +end + +module Poly : sig + + type t + + val constant : Num.num -> t + val variable : var -> t + val addition : t -> t -> t + val product : t -> t -> t + val uminus : t -> t + val get : Monomial.t -> t -> Num.num + val fold : (Monomial.t -> Num.num -> 'a -> 'a) -> t -> 'a -> 'a + + val is_linear : t -> bool + + val add : Monomial.t -> Num.num -> t -> t + +end + +module Vect : sig + + 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 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 + + val from_list : Num.num list -> t + val to_list : t -> Num.num list + + val add : t -> t -> t + val mul : Num.num -> t -> t + +end + +type cstr_compat = {coeffs : Vect.t ; op : op ; cst : Num.num} +and op = Eq | Ge + +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 + +type proof = + | Done + | Step of int * prf_rule * proof + | Enum of int * prf_rule * Vect.t * prf_rule * proof list + +val proof_max_id : proof -> int + +val normalise_proof : int -> proof -> int * proof + +val output_proof : out_channel -> proof -> unit + +val add_proof : prf_rule -> prf_rule -> prf_rule +val mul_proof : Big_int.big_int -> prf_rule -> prf_rule + +module LinPoly : sig + + type t = Vect.t * Num.num + + module MonT : sig + + val clear : unit -> unit + val retrieve : int -> Monomial.t + + end + + val pivot_eq : Vect.var -> + cstr_compat * prf_rule -> + cstr_compat * prf_rule -> (cstr_compat * prf_rule) option + + val linpol_of_pol : Poly.t -> t + +end + +val output_cstr : out_channel -> cstr_compat -> unit + +val opMult : op -> op -> op diff --git a/plugins/micromega/sos.ml b/plugins/micromega/sos.ml index e1ceabe9e..42a41e176 100644 --- a/plugins/micromega/sos.ml +++ b/plugins/micromega/sos.ml @@ -95,7 +95,7 @@ let dim (v:vector) = fst v;; let vector_const c n = if c =/ Int 0 then vector_0 n - else (n,itlist (fun k -> k |-> c) (1--n) undefined :vector);; + else (n,List.fold_right (fun k -> k |-> c) (1--n) undefined :vector);; let vector_cmul c (v:vector) = let n = dim v in @@ -104,7 +104,7 @@ let vector_cmul c (v:vector) = let vector_of_list l = let n = List.length l in - (n,itlist2 (|->) (1--n) l undefined :vector);; + (n,List.fold_right2 (|->) (1--n) l undefined :vector);; (* ------------------------------------------------------------------------- *) (* Matrices; again rows and columns indexed from 1. *) @@ -242,7 +242,7 @@ let string_of_monomial m = if m = monomial_1 then "1" else let vps = List.fold_right (fun (x,k) a -> string_of_varpow x k :: a) (sort humanorder_varpow (graph m)) [] in - end_itlist (fun s t -> s^"*"^t) vps;; + String.concat "*" vps;; let string_of_cmonomial (c,m) = if m = monomial_1 then string_of_num c @@ -310,7 +310,7 @@ let rec poly_of_term t = match t with let sdpa_of_vector (v:vector) = let n = dim v in let strs = List.map (o (decimalize 20) (element v)) (1--n) in - end_itlist (fun x y -> x ^ " " ^ y) strs ^ "\n";; + String.concat " " strs ^ "\n";; (* ------------------------------------------------------------------------- *) (* String for a matrix numbered k, in SDPA sparse format. *) @@ -321,7 +321,7 @@ let sdpa_of_matrix k (m:matrix) = let ms = foldr (fun (i,j) c a -> if i > j then a else ((i,j),c)::a) (snd m) [] in let mss = sort (increasing fst) ms in - itlist (fun ((i,j),c) a -> + List.fold_right (fun ((i,j),c) a -> pfx ^ string_of_int i ^ " " ^ string_of_int j ^ " " ^ decimalize 20 c ^ "\n" ^ a) mss "";; @@ -340,7 +340,7 @@ let sdpa_of_problem comment obj mats = "1\n" ^ string_of_int n ^ "\n" ^ sdpa_of_vector obj ^ - itlist2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a) + List.fold_right2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a) (1--List.length mats) mats "";; (* ------------------------------------------------------------------------- *) @@ -489,11 +489,11 @@ let scale_then = and maximal_element amat acc = foldl (fun maxa m c -> max_num maxa (abs_num c)) acc amat in fun solver obj mats -> - let cd1 = itlist common_denominator mats (Int 1) + let cd1 = List.fold_right common_denominator mats (Int 1) and cd2 = common_denominator (snd obj) (Int 1) in let mats' = List.map (mapf (fun x -> cd1 */ x)) mats and obj' = vector_cmul cd2 obj in - let max1 = itlist maximal_element mats' (Int 0) + let max1 = List.fold_right maximal_element mats' (Int 0) and max2 = maximal_element (snd obj') (Int 0) in let scal1 = pow2 (20-int_of_float(log(float_of_num max1) /. log 2.0)) and scal2 = pow2 (20-int_of_float(log(float_of_num max2) /. log 2.0)) in @@ -551,7 +551,7 @@ let minimal_convex_hull = | (m::ms) -> if in_convex_hull ms m then ms else ms@[m] in let augment m ms = funpow 3 augment1 (m::ms) in fun mons -> - let mons' = itlist augment (List.tl mons) [List.hd mons] in + let mons' = List.fold_right augment (List.tl mons) [List.hd mons] in funpow (List.length mons') augment1 mons';; (* ------------------------------------------------------------------------- *) @@ -612,11 +612,11 @@ let newton_polytope pol = let vars = poly_variables pol in let mons = List.map (fun m -> List.map (fun x -> monomial_degree x m) vars) (dom pol) and ds = List.map (fun x -> (degree x pol + 1) / 2) vars in - let all = itlist (fun n -> allpairs (fun h t -> h::t) (0--n)) ds [[]] + let all = List.fold_right (fun n -> allpairs (fun h t -> h::t) (0--n)) ds [[]] and mons' = minimal_convex_hull mons in let all' = List.filter (fun m -> in_convex_hull mons' (List.map (fun x -> 2 * x) m)) all in - List.map (fun m -> itlist2 (fun v i a -> if i = 0 then a else (v |-> i) a) + List.map (fun m -> List.fold_right2 (fun v i a -> if i = 0 then a else (v |-> i) a) vars m monomial_1) (List.rev all');; (* ------------------------------------------------------------------------- *) @@ -657,8 +657,8 @@ let deration d = foldl (fun a i c -> gcd_num a (numerator c)) (Int 0) (snd l) in (c // (a */ a)),mapa (fun x -> a */ x) l in let d' = List.map adj d in - let a = itlist ((o) lcm_num ( (o) denominator fst)) d' (Int 1) // - itlist ((o) gcd_num ( (o) numerator fst)) d' (Int 0) in + let a = List.fold_right ((o) lcm_num ( (o) denominator fst)) d' (Int 1) // + List.fold_right ((o) gcd_num ( (o) numerator fst)) d' (Int 0) in (Int 1 // a),List.map (fun (c,l) -> (a */ c,l)) d';; (* ------------------------------------------------------------------------- *) @@ -719,7 +719,7 @@ let sdpa_of_blockdiagonal k m = let ents = foldl (fun a (b,i,j) c -> if i > j then a else ((b,i,j),c)::a) [] m in let entss = sort (increasing fst) ents in - itlist (fun ((b,i,j),c) a -> + List.fold_right (fun ((b,i,j),c) a -> pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ " " ^ decimalize 20 c ^ "\n" ^ a) entss "";; @@ -732,10 +732,10 @@ let sdpa_of_blockproblem comment nblocks blocksizes obj mats = "\"" ^ comment ^ "\"\n" ^ string_of_int m ^ "\n" ^ string_of_int nblocks ^ "\n" ^ - (end_itlist (fun s t -> s^" "^t) (List.map string_of_int blocksizes)) ^ + (String.concat " " (List.map string_of_int blocksizes)) ^ "\n" ^ sdpa_of_vector obj ^ - itlist2 (fun k m a -> sdpa_of_blockdiagonal (k - 1) m ^ a) + List.fold_right2 (fun k m a -> sdpa_of_blockdiagonal (k - 1) m ^ a) (1--List.length mats) mats "";; (* ------------------------------------------------------------------------- *) @@ -791,14 +791,14 @@ let blocks blocksizes bm = (fun a (b,i,j) c -> if b = b0 then ((i,j) |-> c) a else a) undefined bm in (((bs,bs),m):matrix)) - (zip blocksizes (1--List.length blocksizes));; + (List.combine blocksizes (1--List.length blocksizes));; (* ------------------------------------------------------------------------- *) (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) (* ------------------------------------------------------------------------- *) let real_positivnullstellensatz_general linf d eqs leqs pol = - let vars = itlist ((o) union poly_variables) (pol::eqs @ List.map fst leqs) [] in + let vars = List.fold_right ((o) union poly_variables) (pol::eqs @ List.map fst leqs) [] in let monoid = if linf then (poly_const num_1,Rational_lt num_1):: @@ -808,16 +808,16 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = let mk_idmultiplier k p = let e = d - multidegree p in let mons = enumerate_monomials e vars in - let nons = zip mons (1--List.length mons) in + let nons = List.combine mons (1--List.length mons) in mons, - itlist (fun (m,n) -> (m |-> ((-k,-n,n) |=> Int 1))) nons undefined in + List.fold_right (fun (m,n) -> (m |-> ((-k,-n,n) |=> Int 1))) nons undefined in let mk_sqmultiplier k (p,c) = let e = (d - multidegree p) / 2 in let mons = enumerate_monomials e vars in - let nons = zip mons (1--List.length mons) in + let nons = List.combine mons (1--List.length mons) in mons, - itlist (fun (m1,n1) -> - itlist (fun (m2,n2) a -> + List.fold_right (fun (m1,n1) -> + List.fold_right (fun (m2,n2) a -> let m = monomial_mul m1 m2 in if n1 > n2 then a else let c = if n1 = n2 then Int 1 else Int 2 in @@ -825,17 +825,17 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = (m |-> equation_add ((k,n1,n2) |=> c) e) a) nons) nons undefined in - let sqmonlist,sqs = unzip(List.map2 mk_sqmultiplier (1--List.length monoid) monoid) - and idmonlist,ids = unzip(List.map2 mk_idmultiplier (1--List.length eqs) eqs) in + let sqmonlist,sqs = List.split(List.map2 mk_sqmultiplier (1--List.length monoid) monoid) + and idmonlist,ids = List.split(List.map2 mk_idmultiplier (1--List.length eqs) eqs) in let blocksizes = List.map List.length sqmonlist in let bigsum = - itlist2 (fun p q a -> epoly_pmul p q a) eqs ids - (itlist2 (fun (p,c) s a -> epoly_pmul p s a) monoid sqs + List.fold_right2 (fun p q a -> epoly_pmul p q a) eqs ids + (List.fold_right2 (fun (p,c) s a -> epoly_pmul p s a) monoid sqs (epoly_of_poly(poly_neg pol))) in let eqns = foldl (fun a m e -> e::a) [] bigsum in let pvs,assig = eliminate_all_equations (0,0,0) eqns in let qvars = (0,0,0)::pvs in - let allassig = itlist (fun v -> (v |-> (v |=> Int 1))) pvs assig in + let allassig = List.fold_right (fun v -> (v |-> (v |=> Int 1))) pvs assig in let mk_matrix v = foldl (fun m (b,i,j) ass -> if b < 0 then m else let c = tryapplyd ass v (Int 0) in @@ -858,8 +858,8 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = else ()); let vec = nice_vector d raw_vec in let blockmat = iter (1,dim vec) - (fun i a -> bmatrix_add (bmatrix_cmul (element vec i) (el i mats)) a) - (bmatrix_neg (el 0 mats)) in + (fun i a -> bmatrix_add (bmatrix_cmul (element vec i) (List.nth mats i)) a) + (bmatrix_neg (List.nth mats 0)) in let allmats = blocks blocksizes blockmat in vec,List.map diag allmats in let vec,ratdias = @@ -867,7 +867,7 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = else tryfind find_rounding (List.map Num.num_of_int (1--31) @ List.map pow2 (5--66)) in let newassigs = - itlist (fun k -> el (k - 1) pvs |-> element vec k) + List.fold_right (fun k -> List.nth pvs (k - 1) |-> element vec k) (1--dim vec) ((0,0,0) |=> Int(-1)) in let finalassigs = foldl (fun a v e -> (v |-> equation_eval newassigs e) a) newassigs @@ -877,17 +877,17 @@ let real_positivnullstellensatz_general linf d eqs leqs pol = undefined p in let mk_sos mons = let mk_sq (c,m) = - c,itlist (fun k a -> (el (k - 1) mons |--> element m k) a) + c,List.fold_right (fun k a -> (List.nth mons (k - 1) |--> element m k) a) (1--List.length mons) undefined in List.map mk_sq in let sqs = List.map2 mk_sos sqmonlist ratdias and cfs = List.map poly_of_epoly ids in let msq = List.filter (fun (a,b) -> b <> []) (List.map2 (fun a b -> a,b) monoid sqs) in - let eval_sq sqs = itlist + let eval_sq sqs = List.fold_right (fun (c,q) -> poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 in let sanity = - itlist (fun ((p,c),s) -> poly_add (poly_mul p (eval_sq s))) msq - (itlist2 (fun p q -> poly_add (poly_mul p q)) cfs eqs + List.fold_right (fun ((p,c),s) -> poly_add (poly_mul p (eval_sq s))) msq + (List.fold_right2 (fun p q -> poly_add (poly_mul p q)) cfs eqs (poly_neg pol)) in if not(is_undefined sanity) then raise Sanity else cfs,List.map (fun (a,b) -> snd a,b) msq;; @@ -913,8 +913,8 @@ let monomial_order = fun m1 m2 -> if m2 = monomial_1 then true else if m1 = monomial_1 then false else let mon1 = dest_monomial m1 and mon2 = dest_monomial m2 in - let deg1 = itlist ((o) (+) snd) mon1 0 - and deg2 = itlist ((o) (+) snd) mon2 0 in + let deg1 = List.fold_right ((o) (+) snd) mon1 0 + and deg2 = List.fold_right ((o) (+) snd) mon2 0 in if deg1 < deg2 then false else if deg1 > deg2 then true else lexorder mon1 mon2;; @@ -929,7 +929,7 @@ let term_of_varpow = let term_of_monomial = fun m -> if m = monomial_1 then Const num_1 else let m' = dest_monomial m in - let vps = itlist (fun (x,k) a -> term_of_varpow x k :: a) m' [] in + let vps = List.fold_right (fun (x,k) a -> term_of_varpow x k :: a) m' [] in end_itlist (fun s t -> Mul (s,t)) vps;; let term_of_cmonomial = @@ -952,203 +952,13 @@ let term_of_sos (pr,sqs) = if sqs = [] then pr else Product(pr,end_itlist (fun a b -> Sum(a,b)) (List.map term_of_sqterm sqs));; -(* ------------------------------------------------------------------------- *) -(* Interface to HOL. *) -(* ------------------------------------------------------------------------- *) -(* -let REAL_NONLINEAR_PROVER translator (eqs,les,lts) = - let eq0 = map (poly_of_term o lhand o concl) eqs - and le0 = map (poly_of_term o lhand o concl) les - and lt0 = map (poly_of_term o lhand o concl) lts in - let eqp0 = map (fun (t,i) -> t,Axiom_eq i) (zip eq0 (0--(length eq0 - 1))) - and lep0 = map (fun (t,i) -> t,Axiom_le i) (zip le0 (0--(length le0 - 1))) - and ltp0 = map (fun (t,i) -> t,Axiom_lt i) (zip lt0 (0--(length lt0 - 1))) in - let keq,eq = partition (fun (p,_) -> multidegree p = 0) eqp0 - and klep,lep = partition (fun (p,_) -> multidegree p = 0) lep0 - and kltp,ltp = partition (fun (p,_) -> multidegree p = 0) ltp0 in - let trivial_axiom (p,ax) = - match ax with - Axiom_eq n when eval undefined p <>/ num_0 -> el n eqs - | Axiom_le n when eval undefined p el n les - | Axiom_lt n when eval undefined p <=/ num_0 -> el n lts - | _ -> failwith "not a trivial axiom" in - try let th = tryfind trivial_axiom (keq @ klep @ kltp) in - CONV_RULE (LAND_CONV REAL_POLY_CONV THENC REAL_RAT_RED_CONV) th - with Failure _ -> - let pol = itlist poly_mul (map fst ltp) (poly_const num_1) in - let leq = lep @ ltp in - let tryall d = - let e = multidegree pol in - let k = if e = 0 then 0 else d / e in - let eq' = map fst eq in - tryfind (fun i -> d,i,real_positivnullstellensatz_general false d eq' leq - (poly_neg(poly_pow pol i))) - (0--k) in - let d,i,(cert_ideal,cert_cone) = deepen tryall 0 in - let proofs_ideal = - map2 (fun q (p,ax) -> Eqmul(term_of_poly q,ax)) cert_ideal eq - and proofs_cone = map term_of_sos cert_cone - and proof_ne = - if ltp = [] then Rational_lt num_1 else - let p = end_itlist (fun s t -> Product(s,t)) (map snd ltp) in - funpow i (fun q -> Product(p,q)) (Rational_lt num_1) in - let proof = end_itlist (fun s t -> Sum(s,t)) - (proof_ne :: proofs_ideal @ proofs_cone) in - print_string("Translating proof certificate to HOL"); - print_newline(); - translator (eqs,les,lts) proof;; -*) -(* ------------------------------------------------------------------------- *) -(* A wrapper that tries to substitute away variables first. *) -(* ------------------------------------------------------------------------- *) -(* -let REAL_NONLINEAR_SUBST_PROVER = - let zero = `&0:real` - and mul_tm = `( * ):real->real->real` - and shuffle1 = - CONV_RULE(REWR_CONV(REAL_ARITH `a + x = (y:real) <=> x = y - a`)) - and shuffle2 = - CONV_RULE(REWR_CONV(REAL_ARITH `x + a = (y:real) <=> x = y - a`)) in - let rec substitutable_monomial fvs tm = - match tm with - Var(_,Tyapp("real",[])) when not (mem tm fvs) -> Int 1,tm - | Comb(Comb(Const("real_mul",_),c),(Var(_,_) as t)) - when is_ratconst c && not (mem t fvs) - -> rat_of_term c,t - | Comb(Comb(Const("real_add",_),s),t) -> - (try substitutable_monomial (union (frees t) fvs) s - with Failure _ -> substitutable_monomial (union (frees s) fvs) t) - | _ -> failwith "substitutable_monomial" - and isolate_variable v th = - match lhs(concl th) with - x when x = v -> th - | Comb(Comb(Const("real_add",_),(Var(_,Tyapp("real",[])) as x)),t) - when x = v -> shuffle2 th - | Comb(Comb(Const("real_add",_),s),t) -> - isolate_variable v(shuffle1 th) in - let make_substitution th = - let (c,v) = substitutable_monomial [] (lhs(concl th)) in - let th1 = AP_TERM (mk_comb(mul_tm,term_of_rat(Int 1 // c))) th in - let th2 = CONV_RULE(BINOP_CONV REAL_POLY_MUL_CONV) th1 in - CONV_RULE (RAND_CONV REAL_POLY_CONV) (isolate_variable v th2) in - fun translator -> - let rec substfirst(eqs,les,lts) = - try let eth = tryfind make_substitution eqs in - let modify = - CONV_RULE(LAND_CONV(SUBS_CONV[eth] THENC REAL_POLY_CONV)) in - substfirst(filter (fun t -> lhand(concl t) <> zero) (map modify eqs), - map modify les,map modify lts) - with Failure _ -> REAL_NONLINEAR_PROVER translator (eqs,les,lts) in - substfirst;; -*) -(* ------------------------------------------------------------------------- *) -(* Overall function. *) -(* ------------------------------------------------------------------------- *) -(* -let REAL_SOS = - let init = GEN_REWRITE_CONV ONCE_DEPTH_CONV [DECIMAL] - and pure = GEN_REAL_ARITH REAL_NONLINEAR_SUBST_PROVER in - fun tm -> let th = init tm in EQ_MP (SYM th) (pure(rand(concl th)));; -*) -(* ------------------------------------------------------------------------- *) -(* Add hacks for division. *) -(* ------------------------------------------------------------------------- *) -(* -let REAL_SOSFIELD = - let inv_tm = `inv:real->real` in - let prenex_conv = - TOP_DEPTH_CONV BETA_CONV THENC - PURE_REWRITE_CONV[FORALL_SIMP; EXISTS_SIMP; real_div; - REAL_INV_INV; REAL_INV_MUL; GSYM REAL_POW_INV] THENC - NNFC_CONV THENC DEPTH_BINOP_CONV `(/\)` CONDS_CELIM_CONV THENC - PRENEX_CONV - and setup_conv = NNF_CONV THENC WEAK_CNF_CONV THENC CONJ_CANON_CONV - and core_rule t = - try REAL_ARITH t - with Failure _ -> try REAL_RING t - with Failure _ -> REAL_SOS t - and is_inv = - let is_div = is_binop `(/):real->real->real` in - fun tm -> (is_div tm or (is_comb tm && rator tm = inv_tm)) && - not(is_ratconst(rand tm)) in - let BASIC_REAL_FIELD tm = - let is_freeinv t = is_inv t && free_in t tm in - let itms = setify(map rand (find_terms is_freeinv tm)) in - let hyps = map (fun t -> SPEC t REAL_MUL_RINV) itms in - let tm' = itlist (fun th t -> mk_imp(concl th,t)) hyps tm in - let itms' = map (curry mk_comb inv_tm) itms in - let gvs = map (genvar o type_of) itms' in - let tm'' = subst (zip gvs itms') tm' in - let th1 = setup_conv tm'' in - let cjs = conjuncts(rand(concl th1)) in - let ths = map core_rule cjs in - let th2 = EQ_MP (SYM th1) (end_itlist CONJ ths) in - rev_itlist (C MP) hyps (INST (zip itms' gvs) th2) in - fun tm -> - let th0 = prenex_conv tm in - let tm0 = rand(concl th0) in - let avs,bod = strip_forall tm0 in - let th1 = setup_conv bod in - let ths = map BASIC_REAL_FIELD (conjuncts(rand(concl th1))) in - EQ_MP (SYM th0) (GENL avs (EQ_MP (SYM th1) (end_itlist CONJ ths)));; -*) -(* ------------------------------------------------------------------------- *) -(* Integer version. *) -(* ------------------------------------------------------------------------- *) -(* -let INT_SOS = - let atom_CONV = - let pth = prove - (`(~(x <= y) <=> y + &1 <= x:int) /\ - (~(x < y) <=> y <= x) /\ - (~(x = y) <=> x + &1 <= y \/ y + &1 <= x) /\ - (x < y <=> x + &1 <= y)`, - REWRITE_TAC[INT_NOT_LE; INT_NOT_LT; INT_NOT_EQ; INT_LT_DISCRETE]) in - GEN_REWRITE_CONV I [pth] - and bub_CONV = GEN_REWRITE_CONV TOP_SWEEP_CONV - [int_eq; int_le; int_lt; int_ge; int_gt; - int_of_num_th; int_neg_th; int_add_th; int_mul_th; - int_sub_th; int_pow_th; int_abs_th; int_max_th; int_min_th] in - let base_CONV = TRY_CONV atom_CONV THENC bub_CONV in - let NNF_NORM_CONV = GEN_NNF_CONV false - (base_CONV,fun t -> base_CONV t,base_CONV(mk_neg t)) in - let init_CONV = - GEN_REWRITE_CONV DEPTH_CONV [FORALL_SIMP; EXISTS_SIMP] THENC - GEN_REWRITE_CONV DEPTH_CONV [INT_GT; INT_GE] THENC - CONDS_ELIM_CONV THENC NNF_NORM_CONV in - let p_tm = `p:bool` - and not_tm = `(~)` in - let pth = TAUT(mk_eq(mk_neg(mk_neg p_tm),p_tm)) in - fun tm -> - let th0 = INST [tm,p_tm] pth - and th1 = NNF_NORM_CONV(mk_neg tm) in - let th2 = REAL_SOS(mk_neg(rand(concl th1))) in - EQ_MP th0 (EQ_MP (AP_TERM not_tm (SYM th1)) th2);; -*) -(* ------------------------------------------------------------------------- *) -(* Natural number version. *) -(* ------------------------------------------------------------------------- *) -(* -let SOS_RULE tm = - let avs = frees tm in - let tm' = list_mk_forall(avs,tm) in - let th1 = NUM_TO_INT_CONV tm' in - let th2 = INT_SOS (rand(concl th1)) in - SPECL avs (EQ_MP (SYM th1) th2);; -*) -(* ------------------------------------------------------------------------- *) -(* Now pure SOS stuff. *) -(* ------------------------------------------------------------------------- *) - -(*prioritize_real();;*) - (* ------------------------------------------------------------------------- *) (* Some combinatorial helper functions. *) (* ------------------------------------------------------------------------- *) let rec allpermutations l = if l = [] then [[]] else - itlist (fun h acc -> List.map (fun t -> h::t) + List.fold_right (fun h acc -> List.map (fun t -> h::t) (allpermutations (subtract l [h])) @ acc) l [];; let changevariables_monomial zoln (m:monomial) = @@ -1165,14 +975,14 @@ let changevariables zoln pol = let sdpa_of_vector (v:vector) = let n = dim v in let strs = List.map (o (decimalize 20) (element v)) (1--n) in - end_itlist (fun x y -> x ^ " " ^ y) strs ^ "\n";; + String.concat " " strs ^ "\n";; let sdpa_of_matrix k (m:matrix) = let pfx = string_of_int k ^ " 1 " in let ms = foldr (fun (i,j) c a -> if i > j then a else ((i,j),c)::a) (snd m) [] in let mss = sort (increasing fst) ms in - itlist (fun ((i,j),c) a -> + List.fold_right (fun ((i,j),c) a -> pfx ^ string_of_int i ^ " " ^ string_of_int j ^ " " ^ decimalize 20 c ^ "\n" ^ a) mss "";; @@ -1184,7 +994,7 @@ let sdpa_of_problem comment obj mats = "1\n" ^ string_of_int n ^ "\n" ^ sdpa_of_vector obj ^ - itlist2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a) + List.fold_right2 (fun k m a -> sdpa_of_matrix (k - 1) m ^ a) (1--List.length mats) mats "";; let run_csdp dbg obj mats = @@ -1224,9 +1034,9 @@ let sumofsquares_general_symmetry tool pol = let sym_eqs = let invariants = List.filter (fun vars' -> - is_undefined(poly_sub pol (changevariables (zip vars vars') pol))) + is_undefined(poly_sub pol (changevariables (List.combine vars vars') pol))) (allpermutations vars) in - let lpns = zip lpps (1--List.length lpps) in + let lpns = List.combine lpps (1--List.length lpps) in let lppcs = List.filter (fun (m,(n1,n2)) -> n1 <= n2) (allpairs @@ -1234,8 +1044,8 @@ let sumofsquares_general_symmetry tool pol = let clppcs = end_itlist (@) (List.map (fun ((m1,m2),(n1,n2)) -> List.map (fun vars' -> - (changevariables_monomial (zip vars vars') m1, - changevariables_monomial (zip vars vars') m2),(n1,n2)) + (changevariables_monomial (List.combine vars vars') m1, + changevariables_monomial (List.combine vars vars') m2),(n1,n2)) invariants) lppcs) in let clppcs_dom = setify(List.map fst clppcs) in @@ -1247,7 +1057,7 @@ let sumofsquares_general_symmetry tool pol = [] -> raise Sanity | [h] -> acc | h::t -> List.map (fun k -> (k |-> Int(-1)) (h |=> Int 1)) t @ acc in - itlist mk_eq eqvcls [] in + List.fold_right mk_eq eqvcls [] in let eqs = foldl (fun a x y -> y::a) [] (itern 1 lpps (fun m1 n1 -> itern 1 lpps (fun m2 n2 f -> @@ -1259,7 +1069,7 @@ let sumofsquares_general_symmetry tool pol = undefined pol)) @ sym_eqs in let pvs,assig = eliminate_all_equations (0,0) eqs in - let allassig = itlist (fun v -> (v |-> (v |=> Int 1))) pvs assig in + let allassig = List.fold_right (fun v -> (v |-> (v |=> Int 1))) pvs assig in let qvars = (0,0)::pvs in let diagents = end_itlist equation_add (List.map (fun i -> apply allassig (i,i)) (1--n)) in @@ -1281,18 +1091,18 @@ let sumofsquares_general_symmetry tool pol = else ()); let vec = nice_vector d raw_vec in let mat = iter (1,dim vec) - (fun i a -> matrix_add (matrix_cmul (element vec i) (el i mats)) a) - (matrix_neg (el 0 mats)) in + (fun i a -> matrix_add (matrix_cmul (element vec i) (List.nth mats i)) a) + (matrix_neg (List.nth mats 0)) in deration(diag mat) in let rat,dia = if pvs = [] then - let mat = matrix_neg (el 0 mats) in + let mat = matrix_neg (List.nth mats 0) in deration(diag mat) else tryfind find_rounding (List.map Num.num_of_int (1--31) @ List.map pow2 (5--66)) in let poly_of_lin(d,v) = - d,foldl(fun a i c -> (el (i - 1) lpps |-> c) a) undefined (snd v) in + d,foldl(fun a i c -> (List.nth lpps (i - 1) |-> c) a) undefined (snd v) in let lins = List.map poly_of_lin dia in let sqs = List.map (fun (d,l) -> poly_mul (poly_const d) (poly_pow l 2)) lins in let sos = poly_cmul rat (end_itlist poly_add sqs) in @@ -1300,325 +1110,3 @@ let sumofsquares_general_symmetry tool pol = let sumofsquares = sumofsquares_general_symmetry csdp;; -(* ------------------------------------------------------------------------- *) -(* Pure HOL SOS conversion. *) -(* ------------------------------------------------------------------------- *) -(* -let SOS_CONV = - let mk_square = - let pow_tm = `(pow)` and two_tm = `2` in - fun tm -> mk_comb(mk_comb(pow_tm,tm),two_tm) - and mk_prod = mk_binop `( * )` - and mk_sum = mk_binop `(+)` in - fun tm -> - let k,sos = sumofsquares(poly_of_term tm) in - let mk_sqtm(c,p) = - mk_prod (term_of_rat(k */ c)) (mk_square(term_of_poly p)) in - let tm' = end_itlist mk_sum (map mk_sqtm sos) in - let th = REAL_POLY_CONV tm and th' = REAL_POLY_CONV tm' in - TRANS th (SYM th');; -*) -(* ------------------------------------------------------------------------- *) -(* Attempt to prove &0 <= x by direct SOS decomposition. *) -(* ------------------------------------------------------------------------- *) -(* -let PURE_SOS_TAC = - let tac = - MATCH_ACCEPT_TAC(REWRITE_RULE[GSYM REAL_POW_2] REAL_LE_SQUARE) ORELSE - MATCH_ACCEPT_TAC REAL_LE_SQUARE ORELSE - (MATCH_MP_TAC REAL_LE_ADD THEN CONJ_TAC) ORELSE - (MATCH_MP_TAC REAL_LE_MUL THEN CONJ_TAC) ORELSE - CONV_TAC(RAND_CONV REAL_RAT_REDUCE_CONV THENC REAL_RAT_LE_CONV) in - REPEAT GEN_TAC THEN REWRITE_TAC[real_ge] THEN - GEN_REWRITE_TAC I [GSYM REAL_SUB_LE] THEN - CONV_TAC(RAND_CONV SOS_CONV) THEN - REPEAT tac THEN NO_TAC;; - -let PURE_SOS tm = prove(tm,PURE_SOS_TAC);; -*) -(* ------------------------------------------------------------------------- *) -(* Examples. *) -(* ------------------------------------------------------------------------- *) - -(***** - -time REAL_SOS - `a1 >= &0 /\ a2 >= &0 /\ - (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + &2) /\ - (a1 * b1 + a2 * b2 = &0) - ==> a1 * a2 - b1 * b2 >= &0`;; - -time REAL_SOS `&3 * x + &7 * a < &4 /\ &3 < &2 * x ==> a < &0`;; - -time REAL_SOS - `b pow 2 < &4 * a * c ==> ~(a * x pow 2 + b * x + c = &0)`;; - -time REAL_SOS - `(a * x pow 2 + b * x + c = &0) ==> b pow 2 >= &4 * a * c`;; - -time REAL_SOS - `&0 <= x /\ x <= &1 /\ &0 <= y /\ y <= &1 - ==> x pow 2 + y pow 2 < &1 \/ - (x - &1) pow 2 + y pow 2 < &1 \/ - x pow 2 + (y - &1) pow 2 < &1 \/ - (x - &1) pow 2 + (y - &1) pow 2 < &1`;; - -time REAL_SOS - `&0 <= b /\ &0 <= c /\ &0 <= x /\ &0 <= y /\ - (x pow 2 = c) /\ (y pow 2 = a pow 2 * c + b) - ==> a * c <= y * x`;; - -time REAL_SOS - `&0 <= x /\ &0 <= y /\ &0 <= z /\ x + y + z <= &3 - ==> x * y + x * z + y * z >= &3 * x * y * z`;; - -time REAL_SOS - `(x pow 2 + y pow 2 + z pow 2 = &1) ==> (x + y + z) pow 2 <= &3`;; - -time REAL_SOS - `(w pow 2 + x pow 2 + y pow 2 + z pow 2 = &1) - ==> (w + x + y + z) pow 2 <= &4`;; - -time REAL_SOS - `x >= &1 /\ y >= &1 ==> x * y >= x + y - &1`;; - -time REAL_SOS - `x > &1 /\ y > &1 ==> x * y > x + y - &1`;; - -time REAL_SOS - `abs(x) <= &1 - ==> abs(&64 * x pow 7 - &112 * x pow 5 + &56 * x pow 3 - &7 * x) <= &1`;; - -time REAL_SOS - `abs(x - z) <= e /\ abs(y - z) <= e /\ &0 <= u /\ &0 <= v /\ (u + v = &1) - ==> abs((u * x + v * y) - z) <= e`;; - -(* ------------------------------------------------------------------------- *) -(* One component of denominator in dodecahedral example. *) -(* ------------------------------------------------------------------------- *) - -time REAL_SOS - `&2 <= x /\ x <= &125841 / &50000 /\ - &2 <= y /\ y <= &125841 / &50000 /\ - &2 <= z /\ z <= &125841 / &50000 - ==> &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= &0`;; - -(* ------------------------------------------------------------------------- *) -(* Over a larger but simpler interval. *) -(* ------------------------------------------------------------------------- *) - -time REAL_SOS - `&2 <= x /\ x <= &4 /\ &2 <= y /\ y <= &4 /\ &2 <= z /\ z <= &4 - ==> &0 <= &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)`;; - -(* ------------------------------------------------------------------------- *) -(* We can do 12. I think 12 is a sharp bound; see PP's certificate. *) -(* ------------------------------------------------------------------------- *) - -time REAL_SOS - `&2 <= x /\ x <= &4 /\ &2 <= y /\ y <= &4 /\ &2 <= z /\ z <= &4 - ==> &12 <= &2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)`;; - -(* ------------------------------------------------------------------------- *) -(* Gloptipoly example. *) -(* ------------------------------------------------------------------------- *) - -(*** This works but normalization takes minutes - -time REAL_SOS - `(x - y - &2 * x pow 4 = &0) /\ &0 <= x /\ x <= &2 /\ &0 <= y /\ y <= &3 - ==> y pow 2 - &7 * y - &12 * x + &17 >= &0`;; - - ***) - -(* ------------------------------------------------------------------------- *) -(* Inequality from sci.math (see "Leon-Sotelo, por favor"). *) -(* ------------------------------------------------------------------------- *) - -time REAL_SOS - `&0 <= x /\ &0 <= y /\ (x * y = &1) - ==> x + y <= x pow 2 + y pow 2`;; - -time REAL_SOS - `&0 <= x /\ &0 <= y /\ (x * y = &1) - ==> x * y * (x + y) <= x pow 2 + y pow 2`;; - -time REAL_SOS - `&0 <= x /\ &0 <= y ==> x * y * (x + y) pow 2 <= (x pow 2 + y pow 2) pow 2`;; - -(* ------------------------------------------------------------------------- *) -(* Some examples over integers and natural numbers. *) -(* ------------------------------------------------------------------------- *) - -time SOS_RULE `!m n. 2 * m + n = (n + m) + m`;; -time SOS_RULE `!n. ~(n = 0) ==> (0 MOD n = 0)`;; -time SOS_RULE `!m n. m < n ==> (m DIV n = 0)`;; -time SOS_RULE `!n:num. n <= n * n`;; -time SOS_RULE `!m n. n * (m DIV n) <= m`;; -time SOS_RULE `!n. ~(n = 0) ==> (0 DIV n = 0)`;; -time SOS_RULE `!m n p. ~(p = 0) /\ m <= n ==> m DIV p <= n DIV p`;; -time SOS_RULE `!a b n. ~(a = 0) ==> (n <= b DIV a <=> a * n <= b)`;; - -(* ------------------------------------------------------------------------- *) -(* This is particularly gratifying --- cf hideous manual proof in arith.ml *) -(* ------------------------------------------------------------------------- *) - -(*** This doesn't now seem to work as well as it did; what changed? - -time SOS_RULE - `!a b c d. ~(b = 0) /\ b * c < (a + 1) * d ==> c DIV d <= a DIV b`;; - - ***) - -(* ------------------------------------------------------------------------- *) -(* Key lemma for injectivity of Cantor-type pairing functions. *) -(* ------------------------------------------------------------------------- *) - -time SOS_RULE - `!x1 y1 x2 y2. ((x1 + y1) EXP 2 + x1 + 1 = (x2 + y2) EXP 2 + x2 + 1) - ==> (x1 + y1 = x2 + y2)`;; - -time SOS_RULE - `!x1 y1 x2 y2. ((x1 + y1) EXP 2 + x1 + 1 = (x2 + y2) EXP 2 + x2 + 1) /\ - (x1 + y1 = x2 + y2) - ==> (x1 = x2) /\ (y1 = y2)`;; - -time SOS_RULE - `!x1 y1 x2 y2. - (((x1 + y1) EXP 2 + 3 * x1 + y1) DIV 2 = - ((x2 + y2) EXP 2 + 3 * x2 + y2) DIV 2) - ==> (x1 + y1 = x2 + y2)`;; - -time SOS_RULE - `!x1 y1 x2 y2. - (((x1 + y1) EXP 2 + 3 * x1 + y1) DIV 2 = - ((x2 + y2) EXP 2 + 3 * x2 + y2) DIV 2) /\ - (x1 + y1 = x2 + y2) - ==> (x1 = x2) /\ (y1 = y2)`;; - -(* ------------------------------------------------------------------------- *) -(* Reciprocal multiplication (actually just ARITH_RULE does these). *) -(* ------------------------------------------------------------------------- *) - -time SOS_RULE `x <= 127 ==> ((86 * x) DIV 256 = x DIV 3)`;; - -time SOS_RULE `x < 2 EXP 16 ==> ((104858 * x) DIV (2 EXP 20) = x DIV 10)`;; - -(* ------------------------------------------------------------------------- *) -(* This is more impressive since it's really nonlinear. See REMAINDER_DECODE *) -(* ------------------------------------------------------------------------- *) - -time SOS_RULE `0 < m /\ m < n ==> ((m * ((n * x) DIV m + 1)) DIV n = x)`;; - -(* ------------------------------------------------------------------------- *) -(* Some conversion examples. *) -(* ------------------------------------------------------------------------- *) - -time SOS_CONV - `&2 * x pow 4 + &2 * x pow 3 * y - x pow 2 * y pow 2 + &5 * y pow 4`;; - -time SOS_CONV - `x pow 4 - (&2 * y * z + &1) * x pow 2 + - (y pow 2 * z pow 2 + &2 * y * z + &2)`;; - -time SOS_CONV `&4 * x pow 4 + - &4 * x pow 3 * y - &7 * x pow 2 * y pow 2 - &2 * x * y pow 3 + - &10 * y pow 4`;; - -time SOS_CONV `&4 * x pow 4 * y pow 6 + x pow 2 - x * y pow 2 + y pow 2`;; - -time SOS_CONV - `&4096 * (x pow 4 + x pow 2 + z pow 6 - &3 * x pow 2 * z pow 2) + &729`;; - -time SOS_CONV - `&120 * x pow 2 - &63 * x pow 4 + &10 * x pow 6 + - &30 * x * y - &120 * y pow 2 + &120 * y pow 4 + &31`;; - -time SOS_CONV - `&9 * x pow 2 * y pow 4 + &9 * x pow 2 * z pow 4 + &36 * x pow 2 * y pow 3 + - &36 * x pow 2 * y pow 2 - &48 * x * y * z pow 2 + &4 * y pow 4 + - &4 * z pow 4 - &16 * y pow 3 + &16 * y pow 2`;; - -time SOS_CONV - `(x pow 2 + y pow 2 + z pow 2) * - (x pow 4 * y pow 2 + x pow 2 * y pow 4 + - z pow 6 - &3 * x pow 2 * y pow 2 * z pow 2)`;; - -time SOS_CONV - `x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z + &3`;; - -(*** I think this will work, but normalization is slow - -time SOS_CONV - `&100 * (x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z) + &212`;; - - ***) - -time SOS_CONV - `&100 * ((&2 * x - &2) pow 2 + (x pow 3 - &8 * x - &2) pow 2) - &588`;; - -time SOS_CONV - `x pow 2 * (&120 - &63 * x pow 2 + &10 * x pow 4) + &30 * x * y + - &30 * y pow 2 * (&4 * y pow 2 - &4) + &31`;; - -(* ------------------------------------------------------------------------- *) -(* Example of basic rule. *) -(* ------------------------------------------------------------------------- *) - -time PURE_SOS - `!x. x pow 4 + y pow 4 + z pow 4 - &4 * x * y * z + x + y + z + &3 - >= &1 / &7`;; - -time PURE_SOS - `&0 <= &98 * x pow 12 + - -- &980 * x pow 10 + - &3038 * x pow 8 + - -- &2968 * x pow 6 + - &1022 * x pow 4 + - -- &84 * x pow 2 + - &2`;; - -time PURE_SOS - `!x. &0 <= &2 * x pow 14 + - -- &84 * x pow 12 + - &1022 * x pow 10 + - -- &2968 * x pow 8 + - &3038 * x pow 6 + - -- &980 * x pow 4 + - &98 * x pow 2`;; - -(* ------------------------------------------------------------------------- *) -(* From Zeng et al, JSC vol 37 (2004), p83-99. *) -(* All of them work nicely with pure SOS_CONV, except (maybe) the one noted. *) -(* ------------------------------------------------------------------------- *) - -PURE_SOS - `x pow 6 + y pow 6 + z pow 6 - &3 * x pow 2 * y pow 2 * z pow 2 >= &0`;; - -PURE_SOS `x pow 4 + y pow 4 + z pow 4 + &1 - &4*x*y*z >= &0`;; - -PURE_SOS `x pow 4 + &2*x pow 2*z + x pow 2 - &2*x*y*z + &2*y pow 2*z pow 2 + -&2*y*z pow 2 + &2*z pow 2 - &2*x + &2* y*z + &1 >= &0`;; - -(**** This is harder. Interestingly, this fails the pure SOS test, it seems. - Yet only on rounding(!?) Poor Newton polytope optimization or something? - But REAL_SOS does finally converge on the second run at level 12! - -REAL_SOS -`x pow 4*y pow 4 - &2*x pow 5*y pow 3*z pow 2 + x pow 6*y pow 2*z pow 4 + &2*x -pow 2*y pow 3*z - &4* x pow 3*y pow 2*z pow 3 + &2*x pow 4*y*z pow 5 + z pow -2*y pow 2 - &2*z pow 4*y*x + z pow 6*x pow 2 >= &0`;; - - ****) - -PURE_SOS -`x pow 4 + &4*x pow 2*y pow 2 + &2*x*y*z pow 2 + &2*x*y*w pow 2 + y pow 4 + z -pow 4 + w pow 4 + &2*z pow 2*w pow 2 + &2*x pow 2*w + &2*y pow 2*w + &2*x*y + -&3*w pow 2 + &2*z pow 2 + &1 >= &0`;; - -PURE_SOS -`w pow 6 + &2*z pow 2*w pow 3 + x pow 4 + y pow 4 + z pow 4 + &2*x pow 2*w + -&2*x pow 2*z + &3*x pow 2 + w pow 2 + &2*z*w + z pow 2 + &2*z + &2*w + &1 >= -&0`;; - -*****) diff --git a/plugins/micromega/sos_lib.ml b/plugins/micromega/sos_lib.ml index 6b8b820ac..6aebc4ca9 100644 --- a/plugins/micromega/sos_lib.ml +++ b/plugins/micromega/sos_lib.ml @@ -9,8 +9,6 @@ open Num -let debugging = ref false;; - (* ------------------------------------------------------------------------- *) (* Comparisons that are reflexive on NaN and also short-circuiting. *) (* ------------------------------------------------------------------------- *) @@ -21,7 +19,6 @@ let (=?) = fun x y -> cmp x y = 0;; let ( cmp x y < 0;; let (<=?) = fun x y -> cmp x y <= 0;; let (>?) = fun x y -> cmp x y > 0;; -let (>=?) = fun x y -> cmp x y >= 0;; (* ------------------------------------------------------------------------- *) (* Combinators. *) @@ -58,49 +55,30 @@ let lcm_num x y = else abs_num((x */ y) // gcd_num x y);; -(* ------------------------------------------------------------------------- *) -(* List basics. *) -(* ------------------------------------------------------------------------- *) - -let rec el n l = - if n = 0 then List.hd l else el (n - 1) (List.tl l);; - - (* ------------------------------------------------------------------------- *) (* Various versions of list iteration. *) (* ------------------------------------------------------------------------- *) -let rec itlist f l b = - match l with - [] -> b - | (h::t) -> f h (itlist f t b);; - let rec end_itlist f l = match l with [] -> failwith "end_itlist" | [x] -> x | (h::t) -> f h (end_itlist f t);; -let rec itlist2 f l1 l2 b = - match (l1,l2) with - ([],[]) -> b - | (h1::t1,h2::t2) -> f h1 h2 (itlist2 f t1 t2 b) - | _ -> failwith "itlist2";; - (* ------------------------------------------------------------------------- *) (* All pairs arising from applying a function over two lists. *) (* ------------------------------------------------------------------------- *) let rec allpairs f l1 l2 = match l1 with - h1::t1 -> itlist (fun x a -> f h1 x :: a) l2 (allpairs f t1 l2) + h1::t1 -> List.fold_right (fun x a -> f h1 x :: a) l2 (allpairs f t1 l2) | [] -> [];; (* ------------------------------------------------------------------------- *) (* String operations (surely there is a better way...) *) (* ------------------------------------------------------------------------- *) -let implode l = itlist (^) l "";; +let implode l = List.fold_right (^) l "";; let explode s = let rec exap n l = @@ -109,13 +87,6 @@ let explode s = exap (String.length s - 1) [];; -(* ------------------------------------------------------------------------- *) -(* Attempting function or predicate applications. *) -(* ------------------------------------------------------------------------- *) - -let can f x = try (f x; true) with Failure _ -> false;; - - (* ------------------------------------------------------------------------- *) (* Repetition of a function. *) (* ------------------------------------------------------------------------- *) @@ -126,36 +97,20 @@ let rec funpow n f x = (* ------------------------------------------------------------------------- *) -(* Replication and sequences. *) +(* Sequences. *) (* ------------------------------------------------------------------------- *) -let rec replicate x n = - if n < 1 then [] - else x::(replicate x (n - 1));; - let rec (--) = fun m n -> if m > n then [] else m::((m + 1) -- n);; (* ------------------------------------------------------------------------- *) (* Various useful list operations. *) (* ------------------------------------------------------------------------- *) -let rec forall p l = - match l with - [] -> true - | h::t -> p(h) && forall p t;; - let rec tryfind f l = match l with [] -> failwith "tryfind" | (h::t) -> try f h with Failure _ -> tryfind f t;; -let index x = - let rec ind n l = - match l with - [] -> failwith "index" - | (h::t) -> if x =? h then n else ind (n + 1) t in - ind 0;; - (* ------------------------------------------------------------------------- *) (* "Set" operations on lists. *) (* ------------------------------------------------------------------------- *) @@ -168,46 +123,16 @@ let rec mem x lis = let insert x l = if mem x l then l else x::l;; -let union l1 l2 = itlist insert l1 l2;; +let union l1 l2 = List.fold_right insert l1 l2;; let subtract l1 l2 = List.filter (fun x -> not (mem x l2)) l1;; -(* ------------------------------------------------------------------------- *) -(* Merging and bottom-up mergesort. *) -(* ------------------------------------------------------------------------- *) - -let rec merge ord l1 l2 = - match l1 with - [] -> l2 - | h1::t1 -> match l2 with - [] -> l1 - | h2::t2 -> if ord h1 h2 then h1::(merge ord t1 l2) - else h2::(merge ord l1 t2);; - - (* ------------------------------------------------------------------------- *) (* Common measure predicates to use with "sort". *) (* ------------------------------------------------------------------------- *) let increasing f x y = f x ? f y;; - -(* ------------------------------------------------------------------------- *) -(* Zipping, unzipping etc. *) -(* ------------------------------------------------------------------------- *) - -let rec zip l1 l2 = - match (l1,l2) with - ([],[]) -> [] - | (h1::t1,h2::t2) -> (h1,h2)::(zip t1 t2) - | _ -> failwith "zip";; - -let rec unzip = - function [] -> [],[] - | ((a,b)::rest) -> let alist,blist = unzip rest in - (a::alist,b::blist);; - (* ------------------------------------------------------------------------- *) (* Iterating functions over lists. *) (* ------------------------------------------------------------------------- *) @@ -443,8 +368,6 @@ let apply f = applyd f (fun x -> failwith "apply");; let tryapplyd f a d = applyd f (fun x -> d) a;; -let defined f x = try apply f x; true with Failure _ -> false;; - (* ------------------------------------------------------------------------- *) (* Undefinition. *) (* ------------------------------------------------------------------------- *) @@ -490,8 +413,6 @@ let graph f = setify (foldl (fun a x y -> (x,y)::a) [] f);; let dom f = setify(foldl (fun a x y -> x::a) [] f);; -let ran f = setify(foldl (fun a x y -> y::a) [] f);; - (* ------------------------------------------------------------------------- *) (* More parser basics. *) (* ------------------------------------------------------------------------- *) @@ -499,7 +420,7 @@ let ran f = setify(foldl (fun a x y -> y::a) [] f);; exception Noparse;; -let isspace,issep,isbra,issymb,isalpha,isnum,isalnum = +let isspace,isnum = let charcode s = Char.code(String.get s 0) in let spaces = " \t\n\r" and separators = ",;" @@ -508,7 +429,7 @@ let isspace,issep,isbra,issymb,isalpha,isnum,isalnum = and alphas = "'abcdefghijklmnopqrstuvwxyz_ABCDEFGHIJKLMNOPQRSTUVWXYZ" and nums = "0123456789" in let allchars = spaces^separators^brackets^symbs^alphas^nums in - let csetsize = itlist ((o) max charcode) (explode allchars) 256 in + let csetsize = List.fold_right ((o) max charcode) (explode allchars) 256 in let ctable = Array.make csetsize 0 in do_list (fun c -> Array.set ctable (charcode c) 1) (explode spaces); do_list (fun c -> Array.set ctable (charcode c) 2) (explode separators); @@ -517,13 +438,8 @@ let isspace,issep,isbra,issymb,isalpha,isnum,isalnum = do_list (fun c -> Array.set ctable (charcode c) 16) (explode alphas); do_list (fun c -> Array.set ctable (charcode c) 32) (explode nums); let isspace c = Array.get ctable (charcode c) = 1 - and issep c = Array.get ctable (charcode c) = 2 - and isbra c = Array.get ctable (charcode c) = 4 - and issymb c = Array.get ctable (charcode c) = 8 - and isalpha c = Array.get ctable (charcode c) = 16 - and isnum c = Array.get ctable (charcode c) = 32 - and isalnum c = Array.get ctable (charcode c) >= 16 in - isspace,issep,isbra,issymb,isalpha,isnum,isalnum;; + and isnum c = Array.get ctable (charcode c) = 32 in + isspace,isnum;; let parser_or parser1 parser2 input = try parser1 input @@ -566,9 +482,6 @@ let rec atleast n prs i = (if n <= 0 then many prs else prs ++ atleast (n - 1) prs >> (fun (h,t) -> h::t)) i;; -let finished input = - if input = [] then 0,input else failwith "Unparsed input";; - (* ------------------------------------------------------------------------- *) let temp_path = Filename.get_temp_dir_name ();; @@ -589,7 +502,7 @@ let strings_of_file filename = (Pervasives.close_in fd; data);; let string_of_file filename = - end_itlist (fun s t -> s^"\n"^t) (strings_of_file filename);; + String.concat "\n" (strings_of_file filename);; let file_of_string filename s = let fd = Pervasives.open_out filename in diff --git a/plugins/micromega/sos_lib.mli b/plugins/micromega/sos_lib.mli new file mode 100644 index 000000000..8b53b8151 --- /dev/null +++ b/plugins/micromega/sos_lib.mli @@ -0,0 +1,79 @@ +(************************************************************************) +(* * The Coq Proof Assistant / The Coq Development Team *) +(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *) +(* 'b) -> ('c -> 'a) -> 'c -> 'b + +val num_1 : Num.num +val pow10 : int -> Num.num +val pow2 : int -> Num.num + +val implode : string list -> string +val explode : string -> string list + +val funpow : int -> ('a -> 'a) -> 'a -> 'a +val tryfind : ('a -> 'b) -> 'a list -> 'b + +type ('a,'b) func = + | Empty + | Leaf of int * ('a*'b) list + | Branch of int * int * ('a,'b) func * ('a,'b) func + +val undefined : ('a, 'b) func +val is_undefined : ('a, 'b) func -> bool +val (|->) : 'a -> 'b -> ('a, 'b) func -> ('a, 'b) func +val (|=>) : 'a -> 'b -> ('a, 'b) func +val choose : ('a, 'b) func -> 'a * 'b +val combine : ('a -> 'a -> 'a) -> ('a -> bool) -> ('b, 'a) func -> ('b, 'a) func -> ('b, 'a) func +val (--) : int -> int -> int list + +val tryapplyd : ('a, 'b) func -> 'a -> 'b -> 'b +val apply : ('a, 'b) func -> 'a -> 'b + +val foldl : ('a -> 'b -> 'c -> 'a) -> 'a -> ('b, 'c) func -> 'a +val foldr : ('a -> 'b -> 'c -> 'c) -> ('a, 'b) func -> 'c -> 'c +val mapf : ('a -> 'b) -> ('c, 'a) func -> ('c, 'b) func + +val undefine : 'a -> ('a, 'b) func -> ('a, 'b) func + +val dom : ('a, 'b) func -> 'a list +val graph : ('a, 'b) func -> ('a * 'b) list + +val union : 'a list -> 'a list -> 'a list +val subtract : 'a list -> 'a list -> 'a list +val sort : ('a -> 'a -> bool) -> 'a list -> 'a list +val setify : 'a list -> 'a list +val increasing : ('a -> 'b) -> 'a -> 'a -> bool +val allpairs : ('a -> 'b -> 'c) -> 'a list -> 'b list -> 'c list + +val gcd_num : Num.num -> Num.num -> Num.num +val lcm_num : Num.num -> Num.num -> Num.num +val numerator : Num.num -> Num.num +val denominator : Num.num -> Num.num +val end_itlist : ('a -> 'a -> 'a) -> 'a list -> 'a + +val (>>) : ('a -> 'b * 'c) -> ('b -> 'd) -> 'a -> 'd * 'c +val (++) : ('a -> 'b * 'c) -> ('c -> 'd * 'e) -> 'a -> ('b * 'd) * 'e + +val a : 'a -> 'a list -> 'a * 'a list +val many : ('a -> 'b * 'a) -> 'a -> 'b list * 'a +val some : ('a -> bool) -> 'a list -> 'a * 'a list +val possibly : ('a -> 'b * 'a) -> 'a -> 'b list * 'a +val isspace : string -> bool +val parser_or : ('a -> 'b) -> ('a -> 'b) -> 'a -> 'b +val isnum : string -> bool +val atleast : int -> ('a -> 'b * 'a) -> 'a -> 'b list * 'a +val listof : ('a -> 'b * 'c) -> ('c -> 'd * 'a) -> string -> 'a -> 'b list * 'c + +val temp_path : string +val string_of_file : string -> string +val file_of_string : string -> string -> unit + +val deepen_until : int -> (int -> 'a) -> int -> 'a +exception TooDeep -- cgit v1.2.3