diff options
| -rw-r--r-- | Makefile.common | 5 | ||||
| -rw-r--r-- | contrib/micromega/certificate.ml | 80 | ||||
| -rw-r--r-- | contrib/micromega/coq_micromega.ml | 537 | ||||
| -rw-r--r-- | contrib/micromega/mfourier.ml | 1579 | ||||
| -rw-r--r-- | contrib/micromega/mutils.ml | 40 | ||||
| -rw-r--r-- | contrib/micromega/vector.ml | 674 | ||||
| -rw-r--r-- | test-suite/micromega/zomicron.v | 3 |
7 files changed, 1450 insertions, 1468 deletions
diff --git a/Makefile.common b/Makefile.common index af380dd63..6b7da1a57 100644 --- a/Makefile.common +++ b/Makefile.common @@ -236,7 +236,7 @@ ROMEGACMO:=\ contrib/romega/refl_omega.cmo contrib/romega/g_romega.cmo MICROMEGACMO:=\ - contrib/micromega/mutils.cmo contrib/micromega/vector.cmo \ + contrib/micromega/mutils.cmo \ contrib/micromega/micromega.cmo contrib/micromega/mfourier.cmo \ contrib/micromega/certificate.cmo \ contrib/micromega/coq_micromega.cmo contrib/micromega/g_micromega.cmo @@ -390,8 +390,7 @@ PARSERCMX:= $(PARSERREQUIRESCMX) $(PARSERCODE:.cmo=.cmx) INTERFACERC:= contrib/interface/vernacrc CSDPCERTCMO:= contrib/micromega/mutils.cmo contrib/micromega/micromega.cmo \ - contrib/micromega/vector.cmo contrib/micromega/mfourier.cmo \ - contrib/micromega/certificate.cmo \ + contrib/micromega/mfourier.cmo contrib/micromega/certificate.cmo \ contrib/micromega/sos.cmo contrib/micromega/csdpcert.cmo CSDPCERTCMX:= $(CSDPCERTCMO:.cmo=.cmx) diff --git a/contrib/micromega/certificate.ml b/contrib/micromega/certificate.ml index f4efcd083..a297c1633 100644 --- a/contrib/micromega/certificate.ml +++ b/contrib/micromega/certificate.ml @@ -76,7 +76,7 @@ struct end -module Poly : +module Poly : (* A polynomial is a map of monomials *) (* This is probably a naive implementation @@ -96,6 +96,7 @@ sig 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 end = struct (*normalisation bug : 0*x ... *) @@ -156,6 +157,8 @@ 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 end @@ -332,10 +335,10 @@ let factorise_linear_cone c = open Mfourier (*module Fourier = Fourier(Vector.VList)(SysSet(Vector.VList))*) (*module Fourier = Fourier(Vector.VSparse)(SysSetAlt(Vector.VSparse))*) -module Fourier = Mfourier.Fourier(Vector.VSparse)(*(SysSetAlt(Vector.VMap))*) +(*module Fourier = Mfourier.Fourier(Vector.VSparse)(*(SysSetAlt(Vector.VMap))*)*) -module Vect = Fourier.Vect -open Fourier.Cstr +(*module Vect = Fourier.Vect*) +(*open Fourier.Cstr*) (* fold_left followed by a rev ! *) @@ -388,8 +391,8 @@ let build_linear_system l = let strict = { coeffs = Vect.from_list ((Big_int unit_big_int):: (List.map (fun (x,y) -> - match y with Mc.Strict -> - Big_int unit_big_int + match y with Mc.Strict -> + Big_int unit_big_int | _ -> Big_int zero_big_int) l)); op = Ge ; cst = Big_int unit_big_int } in (* Add the positivity constraint *) @@ -433,13 +436,41 @@ let make_certificate n_spec cert li = exception Found of Monomial.t - -let raw_certificate l = + +exception Strict + +let primal l = + let vr = ref 0 in + let module Mmn = Map.Make(Monomial) in + + let vect_of_poly map p = + Poly.fold (fun mn vl (map,vect) -> + if mn = Monomial.const + then (map,vect) + else + let (mn,m) = try (Mmn.find mn map,map) with Not_found -> let res = (!vr, Mmn.add mn !vr map) in incr vr ; res in + (m,if sign_num vl = 0 then vect else (mn,vl)::vect)) p (map,[]) in + + let op_op = function Mc.NonStrict -> Ge |Mc.Equal -> Eq | _ -> raise Strict in + + let cmp x y = Pervasives.compare (fst x) (fst y) in + + snd (List.fold_right (fun (p,op) (map,l) -> + let (mp,vect) = vect_of_poly map p in + let cstr = {coeffs = List.sort cmp vect; op = op_op op ; cst = minus_num (Poly.get Monomial.const p)} in + + (mp,cstr::l)) l (Mmn.empty,[])) + +let dual_raw_certificate (l: (Poly.t * Mc.op1) list) = +(* List.iter (fun (p,op) -> Printf.fprintf stdout "%a %s 0\n" Poly.pp p (string_of_op op) ) l ; *) + + let sys = build_linear_system l in - try + + try match Fourier.find_point sys with - | None -> None - | Some cert -> Some (rats_to_ints (Vect.to_list cert)) + | Inr _ -> None + | Inl cert -> Some (rats_to_ints (Vect.to_list cert)) (* should not use rats_to_ints *) with x -> if debug @@ -448,6 +479,21 @@ let raw_certificate l = None +let raw_certificate l = + try + let p = primal l in + match Fourier.find_point p with + | Inr prf -> + if debug then Printf.printf "AProof : %a\n" pp_proof prf ; + let cert = List.map (fun (x,n) -> x+1,n) (fst (List.hd (Proof.mk_proof p prf))) in + if debug then Printf.printf "CProof : %a" Vect.pp_vect cert ; + Some (rats_to_ints (Vect.to_list cert)) + | Inl _ -> None + with Strict -> + (* Fourier elimination should handle > *) + dual_raw_certificate l + + let simple_linear_prover to_constant l = let (lc,li) = List.split l in match raw_certificate lc with @@ -495,7 +541,6 @@ let make_linear_system l = ,monomials) -open Interval let pplus x y = Mc.PEadd(x,y) let pmult x y = Mc.PEmul(x,y) let pconst x = Mc.PEc x @@ -560,11 +605,10 @@ let rec xzlinear_prover planes sys = let smallest_interval = match List.fold_left (fun (x1,i1) (x2,i2) -> - if smaller_itv i1 i2 - then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates + if Itv.smaller_itv i1 i2 + then (x1,i1) else (x2,i2)) (Vect.null,(None,None)) candidates with - | (x,Itv(Some i, Some j)) -> Some(i,x,j) - | (x,Point n) -> Some(n,x,n) + | (x,(Some i, Some j)) -> Some(i,x,j) | x -> None (* This might be a cutting plane *) in match smallest_interval with @@ -613,7 +657,9 @@ and zlinear_enum planes expr clb cub l = let zlinear_prover sys = let candidates = candidates sys in (* Printf.printf "candidates %d" (List.length candidates) ; *) - xzlinear_prover candidates sys + (*let t0 = Sys.time () in*) + let res = xzlinear_prover candidates sys in + (*Printf.printf "Time prover : %f" (Sys.time () -. t0) ;*) res open Sos diff --git a/contrib/micromega/coq_micromega.ml b/contrib/micromega/coq_micromega.ml index 02ff61a19..847ae7baa 100644 --- a/contrib/micromega/coq_micromega.ml +++ b/contrib/micromega/coq_micromega.ml @@ -23,43 +23,68 @@ let time str f x = flush stdout); res -type ('a,'b) formula = + +type tag = int +type 'cst atom = 'cst Micromega.formula + +type 'cst formula = | TT | FF - | X of 'b - | A of 'a * Names.name - | C of ('a,'b) formula * ('a,'b) formula * Names.name - | D of ('a,'b) formula * ('a,'b) formula * Names.name - | N of ('a,'b) formula * Names.name - | I of ('a,'b) formula * ('a,'b) formula * Names.name + | X of Term.constr + | A of 'cst atom * tag * Term.constr + | C of 'cst formula * 'cst formula + | D of 'cst formula * 'cst formula + | N of 'cst formula + | I of 'cst formula * Names.identifier option * 'cst formula + + +let rec pp_formula o f = + match f with + | TT -> output_string o "tt" + | FF -> output_string o "ff" + | X c -> output_string o "X " + | A(_,t,_) -> Printf.fprintf o "A(%i)" t + | C(f1,f2) -> Printf.fprintf o "C(%a,%a)" pp_formula f1 pp_formula f2 + | D(f1,f2) -> Printf.fprintf o "D(%a,%a)" pp_formula f1 pp_formula f2 + | I(f1,n,f2) -> Printf.fprintf o "I(%a%s,%a)" + pp_formula f1 + (match n with + | Some id -> Names.string_of_id id + | None -> "") pp_formula f2 + | N(f) -> Printf.fprintf o "N(%a)" pp_formula f + +let rec ids_of_formula f = + match f with + | I(f1,Some id,f2) -> id::(ids_of_formula f2) + | _ -> [] + +(* obsolete *) + +let tag_formula t f = f (* to be removed *) +(* obsolete *) + +module ISet = Set.Make(struct type t = int let compare : int -> int -> int = Pervasives.compare end) -let none = Names.Anonymous +type 'cst clause = ('cst Micromega.nFormula * tag) list -let tag_formula t f = - match f with - | A(x,_) -> A(x,t) - | C(x,y,_) -> C(x,y,t) - | D(x,y,_) -> D(x,y,t) - | N(x,_) -> N(x,t) - | I(x,y,_) -> I(x,y,t) - | _ -> f +type 'cst cnf = ('cst clause) list -let tt = [] -let ff = [ [] ] +let tt : 'cst cnf = [] +let ff : 'cst cnf = [ [] ] -type ('constant,'contr) sentence = - ('constant Micromega.formula, 'contr) formula +type 'cst mc_cnf = ('cst Micromega.nFormula) Micromega.list Micromega.list -let cnf negate normalise f = - let negate a = - CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> x) cl) (negate a) in +let cnf (negate: 'cst atom -> 'cst mc_cnf) (normalise:'cst atom -> 'cst mc_cnf) (f:'cst formula) = + let negate a t = + CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> (x,t)) cl) (negate a) in - let normalise a = - CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> x) cl) (normalise a) in + let normalise a t = + CoqToCaml.list (fun cl -> CoqToCaml.list (fun x -> (x,t)) cl) (normalise a) in let and_cnf x y = x @ y in - let or_clause_cnf t f = List.map (fun x -> t@x ) f in + + let or_clause_cnf t f = List.map (fun x -> t@x) f in let rec or_cnf f f' = match f with @@ -71,19 +96,17 @@ let cnf negate normalise f = | TT -> if pol then tt else ff (* ?? *) | FF -> if pol then ff else tt (* ?? *) | X p -> if pol then ff else ff (* ?? *) - | A(x,t) -> if pol then normalise x else negate x - | N(e,t) -> xcnf (not pol) e - | C(e1,e2,t) -> + | A(x,t,_) -> if pol then normalise x t else negate x t + | N(e) -> xcnf (not pol) e + | C(e1,e2) -> (if pol then and_cnf else or_cnf) (xcnf pol e1) (xcnf pol e2) - | D(e1,e2,t) -> + | D(e1,e2) -> (if pol then or_cnf else and_cnf) (xcnf pol e1) (xcnf pol e2) - | I(e1,e2,t) -> + | I(e1,_,e2) -> (if pol then or_cnf else and_cnf) (xcnf (not pol) e1) (xcnf pol e2) in xcnf true f - - module M = struct open Coqlib @@ -415,11 +438,11 @@ let parse_q term = Printf.fprintf o "%s%a%s" op _pp l cl - let pp_var = pp_positive let dump_var = dump_positive - let rec pp_expr o e = + 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 @@ -427,7 +450,8 @@ let parse_q term = | 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 + | 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 = @@ -476,9 +500,9 @@ let parse_q term = | Mc.S_In n -> Printf.fprintf o "(S_In %a)%%nat" pp_nat n | Mc.S_Ideal(e,c) -> - Printf.fprintf o "(S_Ideal %a %a)" pp_expr e pp_cone c + Printf.fprintf o "(S_Ideal %a %a)" (pp_expr pp_z) e pp_cone c | Mc.S_Square e -> - Printf.fprintf o "(S_Square %a)" pp_expr e + Printf.fprintf o "(S_Square %a)" (pp_expr pp_z) e | Mc.S_Monoid l -> Printf.fprintf o "(S_Monoid %a)" (pp_list "[" "]" pp_nat) l | Mc.S_Add(e1,e2) -> @@ -514,8 +538,8 @@ let parse_q term = - let pp_cstr o {Mc.flhs = l ; Mc.fop = op ; Mc.frhs = r } = - Printf.fprintf o"(%a %a %a)" pp_expr l pp_op op pp_expr r + 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} = Term.mkApp(Lazy.force coq_Build, @@ -743,21 +767,21 @@ let parse_rexpr = | 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) + | 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 is_prop t = match t with | Names.Anonymous -> true (* Not quite right *) | Names.Name x -> false - let mkC f1 f2 = C(f1,f2,none) - let mkD f1 f2 = D(f1,f2,none) - let mkIff f1 f2 = C(I(f1,f2,none),I(f2,f2,none),none) - let mkI f1 f2 = I(f1,f2,none) + 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)) + let mkI f1 f2 = I(f1,None,f2) let mkformula_binary g term f1 f2 = match f1 , f2 with @@ -765,34 +789,36 @@ let parse_rexpr = | _ -> g f1 f2 let parse_formula parse_atom env term = - let parse_atom env t = try let (at,env) = parse_atom env t in (A(at,none), env) with _ -> (X(t),env) in - let rec xparse_formula env term = + let parse_atom env tg t = try let (at,env) = parse_atom env t in + (A(at,tg,t), env,tg+1) with _ -> (X(t),env,tg) in + + let rec xparse_formula env tg term = match kind_of_term term with | App(l,rst) -> (match rst with | [|a;b|] when l = Lazy.force coq_and -> - let f,env = xparse_formula env a in - let g,env = xparse_formula env b in - mkformula_binary mkC term f g,env + let f,env,tg = xparse_formula env tg a in + let g,env, tg = xparse_formula env tg b in + mkformula_binary mkC term f g,env,tg | [|a;b|] when l = Lazy.force coq_or -> - let f,env = xparse_formula env a in - let g,env = xparse_formula env b in - mkformula_binary mkD term f g,env + let f,env,tg = xparse_formula env tg a in + let g,env,tg = xparse_formula env tg b in + mkformula_binary mkD term f g,env,tg | [|a|] when l = Lazy.force coq_not -> - let (f,env) = xparse_formula env a in (N(f,none), env) + let (f,env,tg) = xparse_formula env tg a in (N(f), env,tg) | [|a;b|] when l = Lazy.force coq_iff -> - let f,env = xparse_formula env a in - let g,env = xparse_formula env b in - mkformula_binary mkIff term f g,env - | _ -> parse_atom env term) + let f,env,tg = xparse_formula env tg a in + let g,env,tg = xparse_formula env tg b in + mkformula_binary mkIff term f g,env,tg + | _ -> parse_atom env tg term) | Prod(typ,a,b) when not (Termops.dependent (mkRel 1) b) -> - let f,env = xparse_formula env a in - let g,env = xparse_formula env b in - mkformula_binary mkI term f g,env - | _ when term = Lazy.force coq_True -> (TT,env) - | _ when term = Lazy.force coq_False -> (FF,env) - | _ -> X(term),env in + let f,env,tg = xparse_formula env tg a in + let g,env,tg = xparse_formula env tg b in + mkformula_binary mkI term f g,env,tg + | _ when term = Lazy.force coq_True -> (TT,env,tg) + | _ when term = Lazy.force coq_False -> (FF,env,tg) + | _ -> X(term),env,tg in xparse_formula env term let coq_TT = lazy @@ -828,11 +854,11 @@ let parse_rexpr = match f with | TT -> mkApp(Lazy.force coq_TT,[| typ|]) | FF -> mkApp(Lazy.force coq_FF,[| typ|]) - | C(x,y,_) -> mkApp(Lazy.force coq_And,[| typ ; xdump x ; xdump y|]) - | D(x,y,_) -> mkApp(Lazy.force coq_Or,[| typ ; xdump x ; xdump y|]) - | I(x,y,_) -> mkApp(Lazy.force coq_Impl,[| typ ; xdump x ; xdump y|]) - | N(x,_) -> mkApp(Lazy.force coq_Neg,[| typ ; xdump x|]) - | A(x,_) -> mkApp(Lazy.force coq_Atom,[| typ ; dump_atom x|]) + | C(x,y) -> mkApp(Lazy.force coq_And,[| typ ; xdump x ; xdump y|]) + | D(x,y) -> mkApp(Lazy.force coq_Or,[| typ ; xdump x ; xdump y|]) + | I(x,_,y) -> mkApp(Lazy.force coq_Impl,[| typ ; xdump x ; xdump y|]) + | N(x) -> mkApp(Lazy.force coq_Neg,[| typ ; xdump x|]) + | A(x,_,_) -> mkApp(Lazy.force coq_Atom,[| typ ; dump_atom x|]) | X(t) -> mkApp(Lazy.force coq_X,[| typ ; t|]) in xdump f @@ -978,22 +1004,22 @@ let pp_q o q = Printf.fprintf o "%a/%a" pp_z q.Micromega.qnum pp_positive q.Micr let rec pp_proof_term o = function - | Micromega.RatProof cone -> Printf.fprintf o "R[%a]" (pp_cone pp_z) cone + | Micromega.RatProof cone -> Printf.fprintf o "R[%a]" (pp_cone pp_z) cone | Micromega.CutProof(e,q,_,p) -> failwith "not implemented" | Micromega.EnumProof(q1,e1,q2,c1,c2,rst) -> Printf.fprintf o "EP[%a,%a,%a,%a,%a,%a]" - pp_q q1 pp_expr e1 pp_q q2 (pp_cone pp_z) c1 (pp_cone pp_z) c2 + pp_q q1 (pp_expr pp_z) e1 pp_q q2 (pp_cone pp_z) c1 (pp_cone pp_z) c2 (pp_list "[" "]" pp_proof_term) rst -let rec parse_hyps parse_arith env hyps = +let rec parse_hyps parse_arith env tg hyps = match hyps with - | [] -> ([],env) + | [] -> ([],env,tg) | (i,t)::l -> - let (lhyps,env) = parse_hyps parse_arith env l in + let (lhyps,env,tg) = parse_hyps parse_arith env tg l in try - let (c,env) = parse_formula parse_arith env t in - ((i,c)::lhyps, env) - with _ -> (lhyps,env) + let (c,env,tg) = parse_formula parse_arith env tg t in + ((i,c)::lhyps, env,tg) + with _ -> (lhyps,env,tg) (*(if debug then Printf.printf "parse_arith : %s\n" x);*) @@ -1001,18 +1027,18 @@ exception ParseError let parse_goal parse_arith env hyps term = (* try*) - let (f,env) = parse_formula parse_arith env term in - let (lhyps,env) = parse_hyps parse_arith env hyps in + let (f,env,tg) = parse_formula parse_arith env 0 term in + let (lhyps,env,tg) = parse_hyps parse_arith env tg hyps in (lhyps,f,env) (* with Failure x -> raise ParseError*) -type ('a, 'b) domain_spec = { +type ('d, 'prf) domain_spec = { typ : Term.constr; (* Z, Q , R *) coeff : Term.constr ; (* Z, Q *) - dump_coeff : 'a -> Term.constr ; + dump_coeff : 'd -> Term.constr ; proof_typ : Term.constr ; - dump_proof : 'b -> Term.constr + dump_proof : 'prf -> Term.constr } let zz_domain_spec = lazy { @@ -1040,6 +1066,34 @@ let rz_domain_spec = lazy { } +let abstract_formula hyps f = + + let rec xabs f = + match f with + | X c -> X c + | A(a,t,term) -> if ISet.mem t hyps then A(a,t,term) else X(term) + | C(f1,f2) -> + (match xabs f1 , xabs f2 with + | X a1 , X a2 -> X (Term.mkApp(Lazy.force coq_and, [|a1;a2|])) + | f1 , f2 -> C(f1,f2) ) + | D(f1,f2) -> + (match xabs f1 , xabs f2 with + | X a1 , X a2 -> X (Term.mkApp(Lazy.force coq_or, [|a1;a2|])) + | f1 , f2 -> D(f1,f2) ) + | N(f) -> + (match xabs f with + | X a -> X (Term.mkApp(Lazy.force coq_not, [|a|])) + | f -> N f) + | I(f1,hyp,f2) -> + (match xabs f1 , hyp, xabs f2 with + | X a1 , Some _ , af2 -> af2 + | X a1 , None , X a2 -> X (Term.mkArrow a1 a2) + | af1 , _ , af2 -> I(af1,hyp,af2) + ) + | x -> x in + xabs f + + let micromega_order_change spec cert cert_typ env ff gl = @@ -1062,42 +1116,27 @@ let micromega_order_change spec cert cert_typ env ff gl = gl -let detect_duplicates cnf wit = - let cnf = CoqToCaml.list (fun x -> x) cnf in - let wit = CoqToCaml.list (fun x -> x) wit in - - let rec xdup cnf wit = - match wit with - | [] -> [] - | w :: wit -> - let sg = sig_of_cone w in - match cnf with - | [] -> [] - | e::cnf -> - let (dups,cnf) = (List.partition (fun x -> same_proof sg e x) cnf) in - dups@(xdup cnf wit) in - xdup cnf wit +type ('a,'prf) prover = { + name : string ; (* name of the prover *) + prover : 'a list -> 'prf option; (* the prover itself *) + hyps : 'prf -> ISet.t; (* extract the indexes of the hypotheses really used in the proof *) + compact : 'prf -> (int -> int) -> 'prf; (* remap the hyp indexes according to function *) + pp_prf : out_channel -> 'prf -> unit ;(* pretting printing of proof *) + pp_f : out_channel -> 'a -> unit +} -let find_witness prover polys1 = - try_any prover polys1 +let find_witness provers polys1 = + + let provers = List.map (fun p -> + (fun l -> + match p.prover l with + | None -> None + | Some prf -> Some(prf,p)) , p.name) provers in + + try_any provers (List.map fst polys1) -let witness_list_with_tags prover l = - - let rec xwitness_list l = - match l with - | [] -> Some([]) - | e::l -> - match find_witness prover (List.map fst e) with - | None -> None - | Some w -> - (match xwitness_list l with - | None -> None - | Some l -> Some (w::l) - ) in - xwitness_list l - -let witness_list_without_tags prover l = +let witness_list_tags prover l = let rec xwitness_list l = match l with @@ -1127,10 +1166,43 @@ let witness_list prover l = xwitness_list l +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 "]" + +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 remap i = + (*Printf.printf "nth (len:%i) %i\n" (List.length old_cl) i;*) + let formula = try fst (List.nth old_cl i) with Failure _ -> failwith "bad old index" in + List.assoc formula new_cl in + let res = prover.compact prf remap in + + if debug then + begin + Printf.printf "\ncompact_proof : %a %a %a -> %a\n" + (pp_ml_list prover.pp_f) (List.map fst old_cl) + prover.pp_prf prf + (pp_ml_list prover.pp_f) (List.map fst new_cl) + prover.pp_prf res ; + flush stdout + end + ; + res in + + (* This seems overkill when List.length cnf_ff = List.length cnf_ff'. Isn't there a 1-1-mapping *) + let cnf_res = List.combine cnf_ff res in + let l = + List.map (fun x -> let (o,p) = List.find (fun lp -> is_sublist x (fst lp)) cnf_res in (o,p,x)) cnf_ff' in + + List.map (fun (o,p,n) -> compact_proof o p n) l -let is_singleton = function [] -> true | [e] -> true | _ -> false - let micromega_tauto negate normalise spec prover env polys1 polys2 gl = let spec = Lazy.force spec in @@ -1139,34 +1211,63 @@ let micromega_tauto negate normalise spec prover env polys1 polys2 gl = (fun (id,f) (cc,ids) -> match f with X _ -> (cc,ids) - | _ -> (I(tag_formula (Names.Name id) f,cc,none), id::ids)) + | _ -> (I(f,Some id,cc), id::ids)) polys1 (polys2,[]) in let cnf_ff = cnf negate normalise ff in if debug then - (Pp.pp (Pp.str "Formula....\n") ; - let formula_typ = (Term.mkApp( Lazy.force coq_Cstr,[| spec.coeff|])) in - let ff = dump_formula formula_typ - (dump_cstr spec.typ spec.dump_coeff) ff in - Pp.pp (Printer.prterm ff) ; Pp.pp_flush ()) ; - - match witness_list_without_tags prover cnf_ff with + begin + Pp.pp (Pp.str "Formula....\n") ; + let formula_typ = (Term.mkApp( Lazy.force coq_Cstr,[| spec.coeff|])) in + let ff = dump_formula formula_typ + (dump_cstr spec.typ spec.dump_coeff) ff in + Pp.pp (Printer.prterm ff) ; Pp.pp_flush () + end; + + match witness_list_tags prover cnf_ff with | None -> Tacticals.tclFAIL 0 (Pp.str "Cannot find witness") gl | Some res -> (*Printf.printf "\nList %i" (List.length res); *) - let (ff,res,ids) = (ff,res,List.map Term.mkVar ids) in - let res' = dump_ml_list (spec.proof_typ) spec.dump_proof res in - (Tacticals.tclTHENSEQ - [ - Tactics.generalize ids; - micromega_order_change spec res' - (Term.mkApp(Lazy.force coq_list,[| spec.proof_typ|])) env ff ; - ]) gl - - -let micromega_gen parse_arith negate normalise spec prover gl = - let concl = Tacmach.pf_concl gl in - let hyps = Tacmach.pf_hyps_types gl in + let hyps = List.fold_left (fun s (cl,(prf,p)) -> + let tags = ISet.fold (fun i s -> try ISet.add (snd (List.nth cl i)) s with Invalid_argument _ -> s) (p.hyps prf) ISet.empty in + ISet.union s tags) ISet.empty (List.combine cnf_ff res) in + + if debug then (Printf.printf "TForm : %a\n" pp_formula ff ; flush stdout; + Printf.printf "Hyps : %a\n" (fun o s -> ISet.fold (fun i _ -> Printf.fprintf o "%i " i) s ()) hyps) ; + + let ff' = abstract_formula hyps ff in + + if debug then + begin + Pp.pp (Pp.str "\nAFormula\n") ; + let formula_typ = (Term.mkApp( Lazy.force coq_Cstr,[| spec.coeff|])) in + let ff' = dump_formula formula_typ + (dump_cstr spec.typ spec.dump_coeff) ff' in + Pp.pp (Printer.prterm ff') ; Pp.pp_flush () + end; + + + let cnf_ff' = cnf negate normalise ff' in + let res' = compact_proofs cnf_ff res cnf_ff' in + + let (ff',res',ids) = (ff',res',List.map Term.mkVar (ids_of_formula ff')) in + + let res' = dump_ml_list (spec.proof_typ) spec.dump_proof res' in + (Tacticals.tclTHENSEQ + [ + Tactics.generalize ids; + micromega_order_change spec res' + (Term.mkApp(Lazy.force coq_list,[| spec.proof_typ|])) env ff' ; + ]) gl + + +let micromega_gen + parse_arith + (negate:'cst atom -> 'cst mc_cnf) + (normalise:'cst atom -> 'cst mc_cnf) + spec prover gl = + let concl = Tacmach.pf_concl gl in + let hyps = Tacmach.pf_hyps_types gl in try let (hyps,concl,env) = parse_goal parse_arith Env.empty hyps concl in let env = Env.elements env in @@ -1235,52 +1336,174 @@ let call_csdpcert_z provername poly = | Mc.False -> (print_string "buggy certificate" ; flush stdout) ;None +let hyps_of_cone prf = + let rec xtract e acc = + match e with + | Mc.S_Pos _ | Mc.S_Z | Mc.S_Square _ -> acc + | Mc.S_In n -> ISet.add (CoqToCaml.nat n) acc + | Mc.S_Ideal(_,c) -> xtract c acc + | Mc.S_Monoid Mc.Nil -> acc + | Mc.S_Monoid (Mc.Cons(i,l)) -> xtract (Mc.S_Monoid l) (ISet.add (CoqToCaml.nat i) acc) + | Mc.S_Add(e1,e2) | Mc.S_Mult(e1,e2) -> xtract e1 (xtract e2 acc) in + + xtract prf ISet.empty + +let compact_cone prf f = + let np n = CamlToCoq.nat (f (CoqToCaml.nat n)) in + + let rec xinterp prf = + match prf with + | Mc.S_Pos _ | Mc.S_Z | Mc.S_Square _ -> prf + | Mc.S_In n -> Mc.S_In (np n) + | Mc.S_Ideal(e,c) -> Mc.S_Ideal(e,xinterp c) + | Mc.S_Monoid l -> Mc.S_Monoid (Mc.map np l) + | Mc.S_Add(e1,e2) -> Mc.S_Add(xinterp e1,xinterp e2) + | Mc.S_Mult(e1,e2) -> Mc.S_Mult(xinterp e1,xinterp e2) in + + xinterp prf + + +let hyps_of_pt pt = + let translate x s = ISet.fold (fun i s -> ISet.add (i - x) s) s ISet.empty in + + let rec xhyps ofset pt acc = + match pt with + | Mc.RatProof p -> ISet.union acc (translate ofset (hyps_of_cone p)) + | Mc.CutProof(_,_,c,pt) -> xhyps (ofset+1) pt (ISet.union (translate (ofset+1) (hyps_of_cone c)) acc) + | Mc.EnumProof(_,_,_,c1,c2,l) -> + let l = CoqToCaml.list (fun x -> x) l in + let s = ISet.union acc (translate (ofset +1) (ISet.union (hyps_of_cone c1) (hyps_of_cone c2))) in + List.fold_left (fun s x -> xhyps (ofset + 1) x s) s l in + + xhyps 0 pt ISet.empty + + +let compact_pt pt f = + let translate ofset x = + if x < ofset then x + else (f (x-ofset) + ofset) in + + let rec compact_pt ofset pt = + match pt with + | Mc.RatProof p -> Mc.RatProof (compact_cone p (translate ofset)) + | Mc.CutProof(x,y,c,pt) -> Mc.CutProof(x,y,compact_cone c (translate (ofset+1)), compact_pt (ofset+1) pt ) + | Mc.EnumProof(a,b,c,c1,c2,l) -> Mc.EnumProof(a,b,c,compact_cone c1 (translate (ofset+1)), compact_cone c2 (translate (ofset+1)), + Mc.map (fun x -> compact_pt (ofset+1) x) l) in + compact_pt 0 pt + +(** Definition of provers *) + +let linear_prover_Z = { + name = "linear prover" ; + prover = lift_ratproof (Certificate.linear_prover Certificate.z_spec) ; + hyps = hyps_of_pt ; + compact = compact_pt ; + pp_prf = pp_proof_term; + pp_f = fun o x -> pp_expr pp_z o (fst' x) +} + +let linear_prover_Q = { + name = "linear prover"; + prover = (Certificate.linear_prover Certificate.q_spec) ; + hyps = hyps_of_cone ; + compact = compact_cone ; + pp_prf = pp_cone pp_q ; + pp_f = fun o x -> pp_expr pp_q o (fst' x) +} + +let linear_prover_R = { + name = "linear prover"; + prover = (Certificate.linear_prover Certificate.z_spec) ; + hyps = hyps_of_cone ; + compact = compact_cone ; + pp_prf = pp_cone pp_z ; + pp_f = fun o x -> pp_expr pp_z o (fst' x) +} + +let non_linear_prover_Q str o = { + name = "real nonlinear prover"; + prover = call_csdpcert_q (str, o); + hyps = hyps_of_cone; + compact = compact_cone ; + pp_prf = pp_cone pp_q ; + pp_f = fun o x -> pp_expr pp_q o (fst' x) +} + +let non_linear_prover_R str o = { + name = "real nonlinear prover"; + prover = call_csdpcert_z (str, o); + hyps = hyps_of_cone; + compact = compact_cone; + pp_prf = pp_cone pp_z; + pp_f = fun o x -> pp_expr pp_z o (fst' x) +} + +let non_linear_prover_Z str o = { + name = "real nonlinear prover"; + prover = lift_ratproof (call_csdpcert_z (str, o)); + hyps = hyps_of_pt; + compact = compact_pt; + pp_prf = pp_proof_term; + pp_f = fun o x -> pp_expr pp_z o (fst' x) +} + +let linear_Z = { + name = "lia"; + prover = Certificate.zlinear_prover ; + hyps = hyps_of_pt; + compact = compact_pt; + pp_prf = pp_proof_term; + pp_f = fun o x -> pp_expr pp_z o (fst' x) +} + + +(** Instantiation of the tactics *) let psatzl_Z gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec - [lift_ratproof - (Certificate.linear_prover Certificate.z_spec), "fourier refutation" ] gl - + [linear_prover_Z ] gl let psatzl_Q gl = micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec - [ Certificate.linear_prover Certificate.q_spec, "fourier refutation" ] gl + [ linear_prover_Q ] gl + let psatz_Q i gl = micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec - [ call_csdpcert_q ("real_nonlinear_prover", Some i), "fourier refutation" ] gl + [ non_linear_prover_Q "real_nonlinear_prover" (Some i) ] gl + let psatzl_R gl = micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec - [ Certificate.linear_prover Certificate.z_spec, "fourier refutation" ] gl - + [ linear_prover_R ] gl let psatz_R i gl = micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec - [ call_csdpcert_z ("real_nonlinear_prover", Some i), "fourier refutation" ] gl + [ non_linear_prover_R "real_nonlinear_prover" (Some i)] gl + - let psatz_Z i gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec - [lift_ratproof (call_csdpcert_z ("real_nonlinear_prover",Some i)), - "fourier refutation" ] gl + [non_linear_prover_Z "real_nonlinear_prover" (Some i) ] gl let sos_Z gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec - [lift_ratproof (call_csdpcert_z ("pure_sos", None)), "pure sos refutation"] gl + [non_linear_prover_Z "pure_sos" None] gl let sos_Q gl = micromega_gen parse_qarith Mc.cnf_negate Mc.cnf_normalise qq_domain_spec - [call_csdpcert_q ("pure_sos", None), "pure sos refutation"] gl + [non_linear_prover_Q "pure_sos" None] gl + let sos_R gl = micromega_gen parse_rarith Mc.cnf_negate Mc.cnf_normalise rz_domain_spec - [call_csdpcert_z ("pure_sos", None), "pure sos refutation"] gl + [non_linear_prover_R "pure_sos" None] gl let xlia gl = micromega_gen parse_zarith Mc.negate Mc.normalise zz_domain_spec - [Certificate.zlinear_prover, "zprover"] gl + [linear_Z] gl + diff --git a/contrib/micromega/mfourier.ml b/contrib/micromega/mfourier.ml index 415d3a3e2..c2eb572c9 100644 --- a/contrib/micromega/mfourier.ml +++ b/contrib/micromega/mfourier.ml @@ -1,667 +1,1014 @@ -(************************************************************************) -(* v * The Coq Proof Assistant / The Coq Development Team *) -(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) -(* \VV/ **************************************************************) -(* // * This file is distributed under the terms of the *) -(* * GNU Lesser General Public License Version 2.1 *) -(************************************************************************) -(* *) -(* Micromega: A reflexive tactic using the Positivstellensatz *) -(* *) -(* Frédéric Besson (Irisa/Inria) 2006-2008 *) -(* *) -(************************************************************************) - -(* Yet another implementation of Fourier *) open Num +module Utils = Mutils + +let map_option = Utils.map_option +let from_option = Utils.from_option + +let debug = false +type ('a,'b) lr = Inl of 'a | Inr of 'b + + +module Vect = + struct + (** [t] is the type of vectors. + A vector [(x1,v1) ; ... ; (xn,vn)] is such that: + - variables indexes are ordered (x1 < ... < xn + - values are all non-zero + *) + type var = int + type t = (var * num) list + +(** [equal v1 v2 = true] if the vectors are syntactically equal. + ([num] is not handled by [Pervasives.equal] *) + + let rec equal v1 v2 = + match v1 , v2 with + | [] , [] -> true + | [] , _ -> false + | _::_ , [] -> false + | (i1,n1)::v1 , (i2,n2)::v2 -> + (i1 = i2) && n1 =/ n2 && equal v1 v2 + + let hash v = + let rec hash i = function + | [] -> i + | (vr,vl)::l -> hash (i + (Hashtbl.hash (vr, float_of_num vl))) l in + Hashtbl.hash (hash 0 v ) -module Cmp = - (* How to compare pairs, lists ... *) -struct - let rec compare_lexical l = - match l with - | [] -> 0 (* Equal *) - | f::l -> - let cmp = f () in - if cmp = 0 then compare_lexical l else cmp - let rec compare_list cmp l1 l2 = - match l1 , l2 with - | [] , [] -> 0 - | [] , _ -> -1 - | _ , [] -> 1 - | e1::l1 , e2::l2 -> - let c = cmp e1 e2 in - if c = 0 then compare_list cmp l1 l2 else c + let null = [] + + let pp_vect o vect = + List.iter (fun (v,n) -> Printf.printf "%sx%i + " (string_of_num n) v) vect + + let from_list (l: num list) = + let rec xfrom_list i l = + match l with + | [] -> [] + | e::l -> + if e <>/ Int 0 + then (i,e)::(xfrom_list (i+1) l) + else xfrom_list (i+1) l in - let hash_list hash l = - let rec xhash res l = - match l with - | [] -> res - | e::l -> xhash ((hash e) lxor res) l in - xhash (Hashtbl.hash []) l + xfrom_list 0 l -end - -module Interval = -struct - (** The type of intervals. **) - type intrvl = Empty | Point of num | Itv of num option * num option + let zero_num = Int 0 + let unit_num = Int 1 + + + let to_list m = + let rec xto_list i l = + match l with + | [] -> [] + | (x,v)::l' -> + if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in + xto_list 0 m + + + let cons i v rst = if v =/ Int 0 then rst else (i,v)::rst + + let rec update i f t = + match t with + | [] -> cons i (f zero_num) [] + | (k,v)::l -> + match Pervasives.compare i k with + | 0 -> cons k (f v) l + | -1 -> cons i (f zero_num) t + | 1 -> (k,v) ::(update i f l) + | _ -> failwith "compare_num" + + let rec set i n t = + match t with + | [] -> cons i n [] + | (k,v)::l -> + match Pervasives.compare i k with + | 0 -> cons k n l + | -1 -> cons i n t + | 1 -> (k,v) :: (set i n l) + | _ -> failwith "compare_num" + + let 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 rec mul z t = + match z with + | Int 0 -> [] + | Int 1 -> t + | _ -> List.map (fun (i,n) -> (i, mult_num z n)) t + + let compare : t -> t -> int = Utils.Cmp.compare_list (fun x y -> Utils.Cmp.compare_lexical + [ + (fun () -> Pervasives.compare (fst x) (fst y)); + (fun () -> compare_num (snd x) (snd y))]) + + (** [tail v vect] returns + - [None] if [v] is not a variable of the vector [vect] + - [Some(vl,rst)] where [vl] is the value of [v] in vector [vect] + and [rst] is the remaining of the vector + We exploit that vectors are ordered lists + *) + let rec tail (v:var) (vect:t) = + match vect with + | [] -> None + | (v',vl)::vect' -> + match Pervasives.compare v' v with + | 0 -> Some (vl,vect) (* Ok, found *) + | -1 -> tail v vect' (* Might be in the tail *) + | _ -> None (* Hopeless *) + + let get v vect = + match tail v vect with + | None -> None + | Some(vl,_) -> Some vl + + + let rec fresh v = + match v with + | [] -> 1 + | [v,_] -> v + 1 + | _::v -> fresh v + + end +open Vect + +(** Implementation of intervals *) +module Itv = +struct - (** - Different intervals can denote the same set of variables e.g., - Point n && Itv (Some n, Some n) - Itv (Some x) (Some y) && Empty if x > y - see the 'belongs_to' function. - **) - - (* The set of numerics that belong to an interval *) - let belongs_to n = function - | Empty -> false - | Point x -> n =/ x - | Itv(Some x, Some y) -> x <=/ n && n <=/ y - | Itv(None,Some y) -> n <=/ y - | Itv(Some x,None) -> x <=/ n - | Itv(None,None) -> true - - let string_of_bound = function - | None -> "oo" - | Some n -> Printf.sprintf "Bd(%s)" (string_of_num n) - - let string_of_intrvl = function - | Empty -> "[]" - | Point n -> Printf.sprintf "[%s]" (string_of_num n) - | Itv(bd1,bd2) -> - Printf.sprintf "[%s,%s]" (string_of_bound bd1) (string_of_bound bd2) - - let pick_closed_to_zero = function - | Empty -> None - | Point n -> Some n - | Itv(None,None) -> Some (Int 0) - | Itv(None,Some i) -> - Some (if (Int 0) <=/ (floor_num i) then Int 0 else floor_num i) - | Itv(Some i,None) -> - Some (if i <=/ (Int 0) then Int 0 else ceiling_num i) - | Itv(Some i,Some j) -> - Some ( - if i <=/ Int 0 && Int 0 <=/ j - then Int 0 - else if ceiling_num i <=/ floor_num j - then ceiling_num i (* why not *) else i) - - type status = - | O | Qonly | Z | Q - - let interval_kind = function - | Empty -> O - | Point n -> if ceiling_num n =/ n then Z else Qonly - | Itv(None,None) -> Z - | Itv(None,Some i) -> if ceiling_num i <>/ i then Q else Z - | Itv(Some i,None) -> if ceiling_num i <>/ i then Q else Z - | Itv(Some i,Some j) -> - if ceiling_num i <>/ i or floor_num j <>/ j then Q else Z - - let empty_z = function - | Empty -> true - | Point n -> ceiling_num n <>/ n - | Itv(None,None) | Itv(None,Some _) | Itv(Some _,None) -> false - | Itv(Some i,Some j) -> ceiling_num i >/ floor_num j - - - let normalise b1 b2 = - match b1 , b2 with - | Some i , Some j -> - (match compare_num i j with - | 1 -> Empty - | 0 -> Point i - | _ -> Itv(b1,b2) - ) - | _ -> Itv(b1,b2) - - - - let min x y = - match x , y with - | None , x | x , None -> x - | Some i , Some j -> Some (min_num i j) + (** The type of intervals is *) + type interval = num option * num option + (** None models the absence of bound i.e. infinity *) + (** As a result, + - None , None -> ]-oo,+oo[ + - None , Some v -> ]-oo,v] + - Some v, None -> [v,+oo[ + - Some v, Some v' -> [v,v'] + Intervals needs to be explicitely normalised. + *) + + type who = Left | Right + + + (** if then interval [itv] is empty, [norm_itv itv] returns [None] + otherwise, it returns [Some itv] *) + + let norm_itv itv = + match itv with + | Some a , Some b -> if a <=/ b then Some itv else None + | _ -> Some itv + + (** [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 = + let (l1,r1) = i1 + and (l2,r2) = i2 in - let max x y = - match x , y with - | None , x | x , None -> x - | Some i , Some j -> Some (max_num i j) - - let inter i1 i2 = - match i1,i2 with - | Empty , _ -> Empty - | _ , Empty -> Empty - | Point n , Point m -> if n =/ m then i1 else Empty - | Point n , Itv (mn,mx) | Itv (mn,mx) , Point n-> - if (match mn with - | None -> true - | Some mn -> mn <=/ n) && - (match mx with - | None -> true - | Some mx -> n <=/ mx) then Point n else Empty - | Itv (min1,max1) , Itv (min2,max2) -> - let bmin = max min1 min2 - and bmax = min max1 max2 in - normalise bmin bmax - - (* a.x >= b*) - let bound_of_constraint (a,b) = - match compare_num a (Int 0) with - | 0 -> - if compare_num b (Int 0) = 1 - then Empty - (*actually this is a contradiction failwith "bound_of_constraint" *) - else Itv (None,None) - | 1 -> Itv (Some (div_num b a),None) - | -1 -> Itv (None, Some (div_num b a)) - | x -> failwith "bound_of_constraint(2)" - - - let bounded x = - match x with - | Itv(None,_) | Itv(_,None) -> false - | _ -> true - - - let range = function - | Empty -> Some (Int 0) - | Point n -> Some (Int (if ceiling_num n =/ n then 1 else 0)) - | Itv(None,_) | Itv(_,None)-> None - | Itv(Some i,Some j) -> Some (floor_num j -/ceiling_num i +/ (Int 1)) - - (* Returns the interval of smallest range *) - let smaller_itv i1 i2 = - match range i1 , range i2 with - | None , _ -> false - | _ , None -> true - | Some i , Some j -> i <=/ j + let inter f o1 o2 = + match o1 , o2 with + | None , None -> None + | Some _ , None -> o1 + | None , Some _ -> o2 + | Some n1 , Some n2 -> Some (f n1 n2) in + + norm_itv (inter max_num l1 l2 , inter min_num r1 r2) + + let range = function + | None,_ | _,None -> None + | Some i,Some j -> Some (floor_num j -/ceiling_num i +/ (Int 1)) + + + let smaller_itv i1 i2 = + match range i1 , range i2 with + | None , _ -> false + | _ , None -> true + | Some i , Some j -> i <=/ j + + +(** [in_bound bnd v] checks whether [v] is within the bounds [bnd] *) +let in_bound bnd v = + let (l,r) = bnd in + match l , r with + | None , None -> true + | None , Some a -> v <=/ a + | Some a , None -> a <=/ v + | Some a , Some b -> a <=/ v && v <=/ b end -open Interval +open Itv +type vector = Vect.t + +type cstr = { coeffs : vector ; bound : interval } +(** 'cstr' is the type of constraints. + {coeffs = v ; bound = (l,r) } models the constraints l <= v <= r +**) + +module ISet = Set.Make(struct type t = int let compare = Pervasives.compare end) + + +module PSet = ISet + + +module System = Hashtbl.Make(Vect) + + type proof = + | Hyp of int + | Elim of var * proof * proof + | And of proof * proof + + + +type system = { + sys : cstr_info ref System.t ; + vars : ISet.t +} +and cstr_info = { + bound : interval ; + prf : proof ; + pos : int ; + neg : int ; +} + + +(** A system of constraints has the form [{sys = s ; vars = v}]. + [s] is a hashtable mapping a normalised vector to a [cstr_info] record where + - [bound] is an interval + - [prf_idx] is the set of hypothese indexes (i.e. constraints in the initial system) used to obtain the current constraint. + In the initial system, each constraint is given an unique singleton proof_idx. + When a new constraint c is computed by a function f(c1,...,cn), its proof_idx is ISet.fold union (List.map (fun x -> x.proof_idx) [c1;...;cn] + - [pos] is the number of positive values of the vector + - [neg] is the number of negative values of the vector + ( [neg] + [pos] is therefore the length of the vector) + [v] is an upper-bound of the set of variables which appear in [s]. +*) + +(** To be thrown when a system has no solution *) +exception SystemContradiction of proof +let hyps prf = + let rec hyps prf acc = + match prf with + | Hyp 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 = + match prf with + | Hyp i -> Printf.fprintf o "H%i" i + | 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 rec pp_list f o l = + match l with + | [] -> () + | e::l -> f o e ; output_string o ";" ; pp_list f o l -(* A set of constraints *) -module Sys(V:Vector.S) (* : Vector.SystemS with module Vect = V*) = -struct - - module Vect = V - - module Cstr = Vector.Cstr(V) - open Cstr - - - module CMap = Map.Make( - struct - type t = Vect.t - let compare = Vect.compare - end) - - module CstrBag = - struct +let pp_iset o s = + output_string o "{" ; + ISet.fold (fun i _ -> Printf.fprintf o "%i " i) s (); + output_string o "}" - type mut_itv = { mutable itv : intrvl} +let pp_pset o s = + output_string o "{" ; + PSet.fold (fun i _ -> Printf.fprintf o "%i " i) s (); + output_string o "}" - type t = mut_itv CMap.t - - exception Contradiction - - let cstr_to_itv cstr = - let (n,l) = V.normalise cstr.coeffs in - if n =/ (Int 0) - then (Vect.null, bound_of_constraint (Int 0,cstr.cst)) (* Might be empty *) - else - match cstr.op with - | Eq -> let n = cstr.cst // n in (l, Point n) - | Ge -> - match compare_num n (Int 0) with - | 0 -> failwith "intrvl_of_constraint" - | 1 -> (l,Itv (Some (cstr.cst // n), None)) - | -1 -> (l, Itv(None,Some (cstr.cst // n))) - | _ -> failwith "cstr_to_itv" - - let empty = CMap.empty +let pp_info o i = pp_itv o i.bound +let pp_cstr o (vect,bnd) = + let (l,r) = bnd in + (match l with + | None -> () + | Some n -> Printf.fprintf o "%s <= " (string_of_num n)) + ; + pp_vect o vect ; + (match r with + | None -> output_string o"\n" + | Some n -> Printf.fprintf o "<=%s\n" (string_of_num n)) +let pp_system o sys= + System.iter (fun vect ibnd -> + pp_cstr o (vect,(!ibnd).bound)) sys - let is_empty = CMap.is_empty - let find_vect v bag = - try - (bag,CMap.find v bag) - with Not_found -> let x = { itv = Itv(None,None)} in (CMap.add v x bag ,x) +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) - let add (v,b) bag = - match b with - | Empty -> raise Contradiction - | Itv(None,None) -> bag - | _ -> - let (bag,intrl) = find_vect v bag in - match inter b intrl.itv with - | Empty -> raise Contradiction - | itv -> intrl.itv <- itv ; bag +(** [merge_cstr_info] takes: + - the intersection of bounds and + - the union of proofs + - [pos] and [neg] fields should be identical *) - exception Found of cstr +let merge_cstr_info i1 i2 = + let { pos = p1 ; neg = n1 ; bound = i1 ; prf = prf1 } = i1 + and { pos = p2 ; neg = n2 ; bound = i2 ; prf = prf2 } = i2 in + assert (p1 = p2 && n1 = n2) ; + match inter i1 i2 with + | None -> None (* Could directly raise a system contradiction exception *) + | Some bnd -> + Some { pos = p1 ; neg = n1 ; bound = bnd ; prf = And(prf1,prf2) } - let find_equation bag = - try - CMap.fold (fun v i () -> - match i.itv with - | Point n -> let e = {coeffs = v ; op = Eq ; cst = n} - in raise (Found e) - | _ -> () ) bag () ; None - with Found c -> Some c +(** [xadd_cstr vect cstr_info] loads an constraint into the system. + The constraint is neither redundant nor contradictory. + @raise SystemContradiction if [cstr_info] returns [None] +*) +let xadd_cstr vect cstr_info sys = + if debug && System.length sys mod 1000 = 0 then (print_string "*" ; flush stdout) ; + try + let info = System.find sys vect in + match merge_cstr_info cstr_info !info with + | None -> raise (SystemContradiction (And(cstr_info.prf, (!info).prf))) + | Some info' -> info := info' + with + | Not_found -> System.replace sys vect (ref cstr_info) + + +type cstr_ext = + | Contradiction (** The constraint is contradictory. + Typically, a [SystemContradiction] exception will be raised. *) + | Redundant (** The constrain is redundant. + Typically, the constraint will be dropped *) + | Cstr of vector * cstr_info (** Taken alone, the constraint is neither contradictory nor redundant. + Typically, it will be added to the constraint system. *) + +(** [normalise_cstr] : vector -> cstr_info -> cstr_ext *) +let normalise_cstr vect cinfo = + match norm_itv cinfo.bound with + | None -> Contradiction + | Some (l,r) -> + match vect with + | [] -> if Itv.in_bound (l,r) (Int 0) then Redundant else Contradiction + | (_,n)::_ -> Cstr( + (if n <>/ Int 1 then List.map (fun (x,nx) -> (x,nx // n)) vect else vect), + let divn x = x // n in + if 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)}) + +(** For compatibility, there an external representation of constraints *) + +type cstr_compat = {coeffs : vector ; op : op ; cst : num} +and op = |Eq | Ge + +let string_of_op = function Eq -> "=" | Ge -> ">=" + + +let eval_op = function + | Eq -> (=/) + | Ge -> (>=/) + +let count v = + let rec count n p v = + match v with + | [] -> (n,p) + | (_,vl)::v -> let sg = sign_num vl in + assert (sg <> 0) ; + if sg = 1 then count n (p+1) v else count (n+1) p v in + count 0 0 v + + +let norm_cstr {coeffs = v ; op = o ; cst = c} idx = + let (n,p) = count v in + + normalise_cstr v {pos = p ; neg = n ; bound = + (match o with + | Eq -> Some c , Some c + | Ge -> Some c , None) ; + prf = Hyp idx } + + +(** [load_system l] takes a list of constraints of type [cstr_compat] + @return a system of constraints + @raise SystemContradiction if a contradiction is found +*) +let load_system l = - let fold f bag acc = - CMap.fold (fun v itv acc -> - match itv.itv with - | Empty | Itv(None,None) -> failwith "fold Empty" - | Itv(None ,Some i) -> - f {coeffs = V.mul (Int (-1)) v ; op = Ge ; cst = minus_num i} acc - | Point n -> f {coeffs = v ; op = Eq ; cst = n} acc - | Itv(x,y) -> - (match x with - | None -> (fun x -> x) - | Some i -> f {coeffs = v ; op = Ge ; cst = i}) - (match y with - | None -> acc - | Some i -> - f {coeffs = V.mul (Int (-1)) v ; op = Ge ; cst = minus_num i} acc - ) ) bag acc - - - let remove l _ = failwith "remove:Not implemented" - - module Map = - Map.Make( - struct - type t = int - let compare : int -> int -> int = Pervasives.compare - end) - - let split f (t:t) = - let res = - fold (fun e m -> let i = f e in - Map.add i (add (cstr_to_itv e) - (try Map.find i m with - Not_found -> empty)) m) t Map.empty in - (fun i -> try Map.find i res with Not_found -> empty) - - type map = (int list * int list) Map.t + let sys = System.create 1000 in + let li = Mutils.mapi (fun e i -> (e,i)) l in + + let vars = List.fold_left (fun vrs (cstr,i) -> + match norm_cstr cstr i with + | Contradiction -> raise (SystemContradiction (Hyp i)) + | Redundant -> vrs + | Cstr(vect,info) -> + xadd_cstr vect info sys ; + List.fold_left (fun s (v,_) -> ISet.add v s) vrs cstr.coeffs) ISet.empty li in + + {sys = sys ;vars = vars} + +let system_list sys = + let { sys = s ; vars = v } = sys in + System.fold (fun k bi l -> (k, !bi)::l) s [] + + +(** [add (v1,c1) (v2,c2) ] + precondition: (c1 <>/ Int 0 && c2 <>/ Int 0) + @return a pair [(v,ln)] such that + [v] is the sum of vector [v1] divided by [c1] and vector [v2] divided by [c2] + Note that the resulting vector is not normalised. +*) + +let add (v1,c1) (v2,c2) = + assert (c1 <>/ Int 0 && c2 <>/ Int 0) ; + + let rec xadd v1 v2 = + match v1 , v2 with + | (x1,n1)::v1' , (x2,n2)::v2' -> + if x1 = x2 + then + let n' = (n1 // c1) +/ (n2 // c2) in + if n' =/ Int 0 then xadd v1' v2' + else + let res = xadd v1' v2' in + (x1,n') ::res + else if x1 < x2 + then let res = xadd v1' v2 in + (x1, n1 // c1)::res + else let res = xadd v1 v2' in + (x2, n2 // c2)::res + | [] , [] -> [] + | [] , _ -> List.map (fun (x,vl) -> (x,vl // c2)) v2 + | _ , [] -> List.map (fun (x,vl) -> (x,vl // c1)) v1 in - let status (b:t) = - let _ , map = fold (fun c ( (idx:int),(res: map)) -> - ( idx + 1, - List.fold_left (fun (res:map) (pos,s) -> - let (lp,ln) = try Map.find pos res with Not_found -> ([],[]) in - match s with - | Vect.Pos -> Map.add pos (idx::lp,ln) res - | Vect.Neg -> - Map.add pos (lp, idx::ln) res) res - (Vect.status c.coeffs))) b (0,Map.empty) in - Map.fold (fun k e res -> (k,e)::res) map [] - - - type it = num CMap.t + let res = xadd v1 v2 in + (res, count res) + +let add (v1,c1) (v2,c2) = + let res = add (v1,c1) (v2,c2) in + (* 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. *) - let iterator x = x - - let element it = failwith "element:Not implemented" +(** [split x vect info (l,m,r)] + @param v is the variable to eliminate + @param l contains constraints such that (e + a*x) // a >= c / a + @param r contains constraints such that (e + a*x) // - a >= c / -a + @param m contains constraints which do not mention [x] +*) + +let split x (vect: vector) info (l,m,r) = + match get x vect with + | None -> (* The constraint does not mention [x], store it in m *) + (l,(vect,info)::m,r) + | Some vl -> (* otherwise *) + + let cons_bound lst bd = + match bd with + | None -> lst + | Some bnd -> (vl,vect,{info with bound = Some bnd,None})::lst in + + let lb,rb = info.bound in + if sign_num vl = 1 + then (cons_bound l lb,m,cons_bound r rb) + else (* sign_num vl = -1 *) + (cons_bound l rb,m,cons_bound r lb) + + +(** [project vr sys] projects system [sys] over the set of variables [ISet.remove vr sys.vars ]. + This is a one step Fourier elimination. +*) +let project vr sys = + + let (l,m,r) = System.fold (fun vect rf l_m_r -> split vr vect !rf l_m_r) sys.sys ([],[],[]) in - end -end + let new_sys = System.create (System.length sys.sys) in + + (* Constraints in [m] belong to the projection - for those [vr] is already projected out *) + List.iter (fun (vect,info) -> System.replace new_sys vect (ref info) ) m ; + + let elim (v1,vect1,info1) (v2,vect2,info2) = + 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 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 + + List.iter(fun l_elem -> List.iter (fun r_elem -> + let (vect,info) = elim l_elem r_elem in + match normalise_cstr vect info with + | Redundant -> () + | Contradiction -> raise (SystemContradiction info.prf) + | Cstr(vect,info) -> xadd_cstr vect info new_sys) r ) l; + {sys = new_sys ; vars = ISet.remove vr sys.vars} + -module Fourier(Vect : Vector.S) = -struct - module Vect = Vect - module Sys = Sys( Vect) - module Cstr = Sys.Cstr - module Bag = Sys.CstrBag +(** [project_using_eq] performs elimination by pivoting using an equation. + This is the counter_part of the [elim] sub-function of [!project]. + @param vr is the variable to be used as pivot + @param c is the coefficient of variable [vr] in vector [vect] + @param len is the length of the equation + @param bound is the bound of the equation + @param prf is the proof of the equation +*) + +let project_using_eq vr c vect bound prf (vect',info') = + match get vr vect' with + | Some c2 -> + let c1 = if c2 >=/ Int 0 then minus_num c else c in + + let c2 = abs_num c2 in + + let (vres,(n,p)) = add (vect,c1) (vect', c2) in + + let cst = bound // c1 in + + let bndres = + let f x = cst +/ x // c2 in + let (l,r) = info'.bound in + (map_option f l , map_option 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 elim_var = project_using_eq vr c vect cst prf in - open Cstr - open Sys + let new_sys = System.create (System.length sys.sys) in - let debug = false + System.iter(fun vect iref -> + let (vect',info') = elim_var (vect,!iref) in + match normalise_cstr vect' info' with + | Redundant -> () + | Contradiction -> raise (SystemContradiction info'.prf) + | Cstr(vect,info') -> xadd_cstr vect info' new_sys) sys.sys ; + + {sys = new_sys ; vars = ISet.remove vr sys.vars} - let print_bag msg b = - print_endline msg; - CstrBag.fold (fun e () -> print_endline (Cstr.string_of_cstr e)) b () + +(** [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(struct type t = int let compare : int -> int -> int = Pervasives.compare end) - let print_bag_file file msg b = - let f = open_out file in - output_string f msg; - CstrBag.fold (fun e () -> - Printf.fprintf f "%s\n" (Cstr.string_of_cstr e)) b () +let pp_map o map = IMap.fold (fun k elt () -> Printf.fprintf o "%i -> %s\n" k (string_of_num elt)) map () +(** [eval_vect map vect] evaluates vector [vect] using the values of [map]. + If [map] binds all the variables of [vect], we get + [eval_vect map [(x1,v1);...;(xn,vn)] = (IMap.find x1 map * v1) + ... + (IMap.find xn map) * vn , []] + The function returns as second argument, a sub-vector consisting in the variables that are not in [map]. *) - (* A system with only inequations -- - *) - let partition i m = - let splitter cstr = compare_num (Vect.get i cstr.coeffs ) (Int 0) in - let split = CstrBag.split splitter m in - (split (-1) , split 0, split 1) +let eval_vect map vect = + let rec xeval_vect vect sum rst = + match vect with + | [] -> (sum,rst) + | (v,vl)::vect -> + try + let val_v = IMap.find v map in + xeval_vect vect (sum +/ (val_v */ vl)) rst + with + Not_found -> xeval_vect vect sum ((v,vl)::rst) in + xeval_vect vect (Int 0) [] + +(** [restrict_bound n sum itv] returns the interval of [x] + given that (fst itv) <= x * n + sum <= (snd itv) *) +let restrict_bound n sum (itv:interval) = + let f x = (x -/ sum) // n in + let l,r = itv in + match sign_num n with + | 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 + + +(** [bound_of_variable map v sys] computes the interval of [v] in + [sys] given a mapping [map] binding all the other variables *) +let bound_of_variable map v sys = + System.fold (fun vect iref bnd -> + let sum,rst = eval_vect map vect in + let vl = match get v rst with + | None -> Int 0 + | Some v -> v in + match inter bnd (restrict_bound vl sum (!iref).bound) with + | None -> failwith "bound_of_variable: impossible" + | Some itv -> itv) sys (None,None) + + +(** [pick_small_value bnd] picks a value being closed to zero within the interval *) +let pick_small_value bnd = + match bnd with + | None , None -> Int 0 + | None , Some i -> if (Int 0) <=/ (floor_num i) then Int 0 else floor_num i + | Some i,None -> if i <=/ (Int 0) then Int 0 else ceiling_num i + | Some i,Some j -> + if i <=/ Int 0 && Int 0 <=/ j + then Int 0 + else if ceiling_num i <=/ floor_num j + then ceiling_num i (* why not *) else i + + +(** [solution s1 sys_l = Some(sn,[(vn-1,sn-1);...; (v1,s1)]@sys_l)] + then [sn] is a system which contains only [black_v] -- if it existed in [s1] + and [sn+1] is obtained by projecting [vn] out of [sn] + @raise SystemContradiction if system [s] has no solution +*) + +let solve_sys black_v choose_eq choose_variable sys sys_l = + + let rec solve_sys sys sys_l = + if debug then Printf.printf "S #%i size %i\n" (System.length sys.sys) (size sys.sys); - (* op of the result is arbitrary Ge *) - let lin_comb n1 c1 n2 c2 = - { coeffs = Vect.lin_comb n1 c1.coeffs n2 c2.coeffs ; - op = Ge ; - cst = (n1 */ c1.cst) +/ (n2 */ c2.cst)} - - (* BUG? : operator of the result ? *) - - let combine_project i c1 c2 = - let p = Vect.get i c1.coeffs - and n = Vect.get i c2.coeffs in - assert (n </ Int 0 && p >/ Int 0) ; - let nopp = minus_num n in - let c =lin_comb nopp c1 p c2 in - let op = if c1.op = Ge || c2.op = Ge then Ge else Eq in - CstrBag.cstr_to_itv {coeffs = c.coeffs ; op = op ; cst= c.cst } - - - let project i m = - let (neg,zero,pos) = partition i m in - let project1 cpos acc = - CstrBag.fold (fun cneg res -> - CstrBag.add (combine_project i cpos cneg) res) neg acc in - (CstrBag.fold project1 pos zero) - - (* Given a vector [x1 -> v1; ... ; xn -> vn] - and a constraint {x1 ; .... xn >= c } - *) - let evaluate_constraint i map cstr = - let {coeffs = _coeffs ; op = _op ; cst = _cst} = cstr in - let vi = Vect.get i _coeffs in - let v = Vect.set i (Int 0) _coeffs in - (vi, _cst -/ Vect.dotp map v) - - - let rec bounds m itv = - match m with - | [] -> itv - | e::m -> bounds m (inter itv (bound_of_constraint e)) - - - - let compare_status (i,(lp,ln)) (i',(lp',ln')) = - let cmp = Pervasives.compare - ((List.length lp) * (List.length ln)) - ((List.length lp') * (List.length ln')) in - if cmp = 0 - then Pervasives.compare i i' - else cmp - - let cardinal m = CstrBag.fold (fun _ x -> x + 1) m 0 - - let lightest_projection l c m = - let bound = c in - if debug then (Printf.printf "l%i" bound; flush stdout) ; - let rec xlight best l = - match l with - | [] -> best - | i::l -> - let proj = (project i m) in - let cproj = cardinal proj in - (*Printf.printf " p %i " cproj; flush stdout;*) - match best with - | None -> - if cproj < bound - then Some(cproj,proj,i) - else xlight (Some(cproj,proj,i)) l - | Some (cbest,_,_) -> - if cproj < cbest - then - if cproj < bound then Some(cproj,proj,i) - else xlight (Some(cproj,proj,i)) l - else xlight best l in - match xlight None l with - | None -> None - | Some(_,p,i) -> Some (p,i) - + let eqs = choose_eq sys in + try + let (v,vect,cst,ln) = fst (List.find (fun ((v,_,_,_),_) -> v <> black_v) eqs) in + if debug then + (Printf.printf "\nE %a = %s variable %i\n" pp_vect vect (string_of_num cst) v ; + flush stdout); + let sys' = elim_var_using_eq v vect cst ln sys in + solve_sys sys' ((v,sys)::sys_l) + with Not_found -> + let vars = choose_variable sys in + try + let (v,est) = (List.find (fun (v,_) -> v <> black_v) vars) in + if debug then (Printf.printf "\nV : %i esimate %f\n" v est ; flush stdout) ; + let sys' = project v sys in + solve_sys sys' ((v,sys)::sys_l) + with Not_found -> (* we are done *) Inl (sys,sys_l) in + solve_sys sys sys_l - exception Equality of cstr - let find_equality m = Bag.find_equation m - +let solve black_v choose_eq choose_variable cstrs = - let pivot (n,v) eq ge = - assert (eq.op = Eq) ; - let res = - match - compare_num v (Int 0), - compare_num (Vect.get n ge.coeffs) (Int 0) - with - | 0 , _ -> failwith "Buggy" - | _ ,0 -> (CstrBag.cstr_to_itv ge) - | 1 , -1 -> combine_project n eq ge - | -1 , 1 -> combine_project n ge eq - | 1 , 1 -> - combine_project n ge - {coeffs = Vect.mul (Int (-1)) eq.coeffs; - op = eq.op ; - cst = minus_num eq.cst} - | -1 , -1 -> - combine_project n - {coeffs = Vect.mul (Int (-1)) eq.coeffs; - op = eq.op ; cst = minus_num eq.cst} ge - | _ -> failwith "pivot" in - res - - let check_cstr v c = - let {coeffs = _coeffs ; op = _op ; cst = _cst} = c in - let vl = Vect.dotp v _coeffs in - match _op with - | Eq -> vl =/ _cst - | Ge -> vl >= _cst - - - let forall p sys = try - CstrBag.fold (fun c () -> if p c then () else raise Not_found) sys (); true - with Not_found -> false - - - let check_sys v sys = forall (check_cstr v) sys - - let check_null_cstr c = - let {coeffs = _coeffs ; op = _op ; cst = _cst} = c in - match _op with - | Eq -> (Int 0) =/ _cst - | Ge -> (Int 0) >= _cst - - let check_null sys = forall check_null_cstr sys - - - let optimise_ge - quick_check choose choose_idx return_empty return_ge return_eq m = - let c = cardinal m in - let bound = 2 * c in - if debug then (Printf.printf "optimise_ge: %i\n" c; flush stdout); - - let rec xoptimise m = - if debug then (Printf.printf "x%i" (cardinal m) ; flush stdout); - if debug then (print_bag "xoptimise" m ; flush stdout); - if quick_check m - then return_empty m - else - match find_equality m with - | None -> xoptimise_ge m - | Some eq -> xoptimise_eq eq m - - and xoptimise_ge m = - begin - let c = cardinal m in - let l = List.map fst (List.sort compare_status (CstrBag.status m)) in - let idx = choose bound l c m in - match idx with - | None -> return_empty m - | Some (proj,i) -> - match xoptimise proj with - | None -> None - | Some mapping -> return_ge m i mapping - end - and xoptimise_eq eq m = - let l = List.map fst (Vect.status eq.coeffs) in - match choose_idx l with - | None -> (*if l = [] then None else*) return_empty m - | Some i -> - let p = (i,Vect.get i eq.coeffs) in - let m' = CstrBag.fold - (fun ge res -> CstrBag.add (pivot p eq ge) res) m CstrBag.empty in - match xoptimise ( m') with - | None -> None - | Some mapp -> return_eq m eq i mapp in - try - let res = xoptimise m in res - with CstrBag.Contradiction -> (*print_string "contradiction" ;*) None - - - - let minimise m = - let opt_zero_choose bound l c m = - if c > bound - then lightest_projection l c m - else match l with - | [] -> None - | i::_ -> Some (project i m, i) in - - let choose_idx = function [] -> None | x::l -> Some x in - - let opt_zero_return_empty m = Some Vect.null in - - - let opt_zero_return_ge m i mapping = - let (it:intrvl) = CstrBag.fold (fun cstr itv -> Interval.inter - (bound_of_constraint (evaluate_constraint i mapping cstr)) itv) m - (Itv (None, None)) in - match pick_closed_to_zero it with - | None -> print_endline "Cannot pick" ; None - | Some v -> - let res = (Vect.set i v mapping) in - if debug - then Printf.printf "xoptimise res %i [%s]" i (Vect.string res) ; - Some res in - - let opt_zero_return_eq m eq i mapp = - let (a,b) = evaluate_constraint i mapp eq in - Some (Vect.set i (div_num b a) mapp) in - - optimise_ge check_null opt_zero_choose - choose_idx opt_zero_return_empty opt_zero_return_ge opt_zero_return_eq m - - let normalise cstr = [CstrBag.cstr_to_itv cstr] - - let find_point l = - (* List.iter (fun e -> print_endline (Cstr.string_of_cstr e)) l;*) - try - let m = List.fold_left (fun sys e -> CstrBag.add (CstrBag.cstr_to_itv e) sys) - CstrBag.empty l in - match minimise m with - | None -> None - | Some res -> - if debug then Printf.printf "[%s]" (Vect.string res); - Some res - with CstrBag.Contradiction -> None + let sys = load_system cstrs in +(* Printf.printf "solve :\n %a" pp_system sys.sys ; *) + solve_sys black_v choose_eq choose_variable sys [] + with SystemContradiction prf -> Inr prf + + +(** The purpose of module [EstimateElimVar] is to try to estimate the cost of eliminating a variable. + The output is an ordered list of (variable,cost). +*) + +module EstimateElimVar = +struct + type sys_list = (vector * cstr_info) list + + let abstract_partition (v:int) (l: sys_list) = + + let rec xpart (l:sys_list) (ltl:sys_list) (n:int list) (z:int) (p:int list) = + match l with + | [] -> (ltl, n,z,p) + | (l1,info) ::rl -> + match l1 with + | [] -> xpart rl (([],info)::ltl) n (info.neg+info.pos+z) p + | (vr,vl)::rl1 -> + if v = vr + then + let cons_bound lst bd = + match bd with + | None -> lst + | Some bnd -> info.neg+info.pos::lst in + + let lb,rb = info.bound in + if sign_num vl = 1 + then xpart rl ((rl1,info)::ltl) (cons_bound n lb) z (cons_bound p rb) + else xpart rl ((rl1,info)::ltl) (cons_bound n rb) z (cons_bound p lb) + else + (* the variable is greater *) + xpart rl ((l1,info)::ltl) n (info.neg+info.pos+z) p + + in + let (sys',n,z,p) = xpart l [] [] 0 [] in + + let ln = float_of_int (List.length n) in + let sn = float_of_int (List.fold_left (+) 0 n) in + let lp = float_of_int (List.length p) in + let sp = float_of_int (List.fold_left (+) 0 p) in + (sys', float_of_int z +. lp *. sn +. ln *. sp -. lp*.ln) + + + let choose_variable sys = + let {sys = s ; vars = v} = sys in + + let sl = system_list sys in + + let evals = fst + (ISet.fold (fun v (eval,s) -> let ts,vl = abstract_partition v s in + ((v,vl)::eval, ts)) v ([],sl)) in + + List.sort (fun x y -> Pervasives.compare (snd x) (snd y) ) evals - let find_q_interval_for x m = - if debug then Printf.printf "find_q_interval_for %i\n" x ; +end +open EstimateElimVar + +(** The module [EstimateElimEq] is similar to [EstimateElimVar] but it orders equations. +*) +module EstimateElimEq = +struct + + let itv_point bnd = + match bnd with + |(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 choose bound l c m = - let rec xchoose l = + let rec unroll_until v l = match l with - | [] -> None - | i::l -> if i = x then xchoose l else Some (project i m,i) in - xchoose l in - - let rec choose_idx = function - [] -> None - | e::l -> if e = x then choose_idx l else Some e in - - let return_empty m = (* Beurk *) - (* returns the interval of x *) - Some (CstrBag.fold (fun cstr itv -> - let i = if cstr.op = Eq - then Point (cstr.cst // Vect.get x cstr.coeffs) - else if Vect.is_null (Vect.set x (Int 0) cstr.coeffs) - then bound_of_constraint (Vect.get x cstr.coeffs , cstr.cst) - else itv - in - Interval.inter i itv) m (Itv (None, None))) in - - let return_ge m i res = Some res in - - let return_eq m eq i res = Some res in - - try - optimise_ge - (fun x -> false) choose choose_idx return_empty return_ge return_eq m - with CstrBag.Contradiction -> None - - - let find_q_intervals sys = - let variables = - List.map fst (List.sort compare_status (CstrBag.status sys)) in - List.map (fun x -> (x,find_q_interval_for x sys)) variables - - let pp_option f o = function - None -> Printf.fprintf o "None" - | Some x -> Printf.fprintf o "Some %a" f x - - let optimise vect sys = - (* we have to modify the system with a dummy variable *) - let fresh = - List.fold_left (fun fr c -> Pervasives.max fr (Vect.fresh c.coeffs)) 0 sys in - assert (List.for_all (fun x -> Vect.get fresh x.coeffs =/ Int 0) sys); - let cstr = { - coeffs = Vect.set fresh (Int (-1)) vect ; - op = Eq ; - cst = (Int 0)} in - try - find_q_interval_for fresh - (List.fold_left - (fun bg c -> CstrBag.add (CstrBag.cstr_to_itv c) bg) - CstrBag.empty (cstr::sys)) - with CstrBag.Contradiction -> None - - - let optimise vect sys = - let res = optimise vect sys in - if debug - then Printf.printf "optimise %s -> %a\n" - (Vect.string vect) (pp_option (fun o x -> Printf.printf "%s" (string_of_intrvl x))) res - ; res - - let find_Q_interval sys = - try - let sys = - (List.fold_left - (fun bg c -> CstrBag.add (CstrBag.cstr_to_itv c) bg) CstrBag.empty sys) in - let candidates = - List.fold_left - (fun l (x,i) -> match i with - None -> (x,Empty)::l - | Some i -> (x,i)::l) [] (find_q_intervals sys) in - match List.fold_left - (fun (x1,i1) (x2,i2) -> - if smaller_itv i1 i2 - then (x1,i1) else (x2,i2)) (-1,Itv(None,None)) candidates - with - | (i,Empty) -> None - | (x,Itv(Some i, Some j)) -> Some(i,x,j) - | (x,Point n) -> Some(n,x,n) - | _ -> None - with CstrBag.Contradiction -> None + | [] -> (false,[]) + | (i,_)::rl -> if i = v + then (true,rl) + else if i < v then unroll_until v rl else (false,l) + + + let choose_primal_equation eqs sys_l = + + let is_primal_equation_var v = + List.fold_left (fun (nb_eq,nb_cst) (vect,info) -> + if fst (unroll_until v vect) + then if itv_point info.bound then (nb_eq + 1,nb_cst) else (nb_eq,nb_cst) + else (nb_eq,nb_cst)) (0,0) sys_l in + + let rec find_var vect = + match vect with + | [] -> None + | (i,_)::vect -> + let (nb_eq,nb_cst) = is_primal_equation_var i in + if nb_eq = 2 && nb_cst = 0 + then Some i else find_var vect in + + let rec find_eq_var eqs = + match eqs with + | [] -> None + | (vect,a,prf,ln)::l -> + match find_var vect with + | None -> find_eq_var l + | Some r -> Some (r,vect,a,prf,ln) + in + + + find_eq_var eqs + + + + + let choose_equality_var sys = + + let sys_l = system_list sys in + + let equalities = List.fold_left + (fun l (vect,info) -> + match info.bound with + | Some a , Some b -> + if a =/ b then (* This an equation *) + (vect,a,info.prf,info.neg+info.pos)::l else l + | _ -> l + ) [] sys_l in + + let rec estimate_cost v ct sysl acc tlsys = + match sysl with + | [] -> (acc,tlsys) + | (l,info)::rsys -> + let ln = info.pos + info.neg in + let (b,l) = unroll_until v l in + match b with + | true -> + if itv_point info.bound + then estimate_cost v ct rsys (acc+ln) ((l,info)::tlsys) (* this is free *) + else estimate_cost v ct rsys (acc+ln+ct) ((l,info)::tlsys) (* should be more ? *) + | false -> estimate_cost v ct rsys (acc+ln) ((l,info)::tlsys) in + + match choose_primal_equation equalities sys_l with + | None -> + let cost_eq eq const prf ln acc_costs = + + let rec cost_eq eqr sysl costs = + match eqr with + | [] -> costs + | (v,_) ::eqr -> let (cst,tlsys) = estimate_cost v (ln-1) sysl 0 [] in + cost_eq eqr tlsys (((v,eq,const,prf),cst)::costs) in + cost_eq eq sys_l acc_costs in + + let all_costs = List.fold_left (fun all_costs (vect,const,prf,ln) -> cost_eq vect const prf ln all_costs) [] equalities in + + (* pp_list (fun o ((v,eq,_,_),cst) -> Printf.fprintf o "((%i,%a),%i)\n" v pp_vect eq cst) stdout all_costs ; *) + + List.sort (fun x y -> Pervasives.compare (snd x) (snd y) ) all_costs + | Some (v,vect, const,prf,_) -> [(v,vect,const,prf),0] end +open EstimateElimEq + + + +module Fourier = +struct + + let optimise vect l = + (* We add a dummy (fresh) variable for vector *) + let fresh = + List.fold_left (fun fr c -> Pervasives.max fr (Vect.fresh c.coeffs)) 0 l in + let cstr = { + coeffs = Vect.set fresh (Int (-1)) vect ; + op = Eq ; + cst = (Int 0)} in + match solve fresh choose_equality_var choose_variable (cstr::l) with + | Inr prf -> None (* This is an unsatisfiability proof *) + | Inl (s,_) -> + try + Some (bound_of_variable IMap.empty fresh s.sys) + with + x -> Printf.printf "optimise Exception : %s" (Printexc.to_string x) ; None + + + let find_point cstrs = + + match solve max_int choose_equality_var choose_variable cstrs with + | Inr prf -> Inr prf + | Inl (_,l) -> + + let rec rebuild_solution l map = + match l with + | [] -> map + | (v,e)::l -> + let itv = bound_of_variable map v e.sys in + let map = IMap.add v (pick_small_value itv) map in + rebuild_solution l map + in + + let map = rebuild_solution l IMap.empty in + let vect = List.rev (IMap.fold (fun v i vect -> (v,i)::vect) map []) in +(* Printf.printf "SOLUTION %a" pp_vect vect ; *) + let res = Inl vect in + res + + +end + + +module Proof = +struct + + + + +(** A proof term in the sense of a ZMicromega.RatProof is a positive combination of the hypotheses which leads to a contradiction. + The proofs constructed by Fourier elimination are more like execution traces: + - certain facts are recorded but are useless + - certain inferences are implicit. + The following code implements proof reconstruction. +*) + let add x y = fst (add x y) + + + let forall_pairs f l1 l2 = + List.fold_left (fun acc e1 -> + List.fold_left (fun acc e2 -> + match f e1 e2 with + | None -> acc + | Some v -> v::acc) acc l2) [] l1 + + + let add_op x y = + match x , y with + | Eq , Eq -> Eq + | _ -> Ge + + + let pivot v (p1,c1) (p2,c2) = + let {coeffs = v1 ; op = op1 ; cst = n1} = c1 + and {coeffs = v2 ; op = op2 ; cst = n2} = c2 in + + match Vect.get v v1 , Vect.get v v2 with + | None , _ | _ , None -> None + | Some a , Some b -> + if (sign_num a) * (sign_num b) = -1 + then Some (add (p1,abs_num a) (p2,abs_num b) , + {coeffs = add (v1,abs_num a) (v2,abs_num b) ; + op = add_op op1 op2 ; + cst = n1 // (abs_num a) +/ n2 // (abs_num b) }) + else if op1 = Eq + then Some (add (p1,minus_num (a // b)) (p2,Int 1), + {coeffs = add (v1,minus_num (a// b)) (v2 ,Int 1) ; + op = add_op op1 op2; + cst = n1 // (minus_num (a// b)) +/ n2 // (Int 1)}) + else if op2 = Eq + then + Some (add (p2,minus_num (b // a)) (p1,Int 1), + {coeffs = add (v2,minus_num (b// a)) (v1 ,Int 1) ; + op = add_op op1 op2; + cst = n2 // (minus_num (b// a)) +/ n1 // (Int 1)}) + else None (* op2 could be Eq ... this might happen *) + + + let normalise_proofs l = + List.fold_left (fun acc (prf,cstr) -> + match acc with + | Inr _ -> acc (* I already found a contradiction *) + | Inl acc -> + match norm_cstr cstr 0 with + | Redundant -> Inl acc + | Contradiction -> Inr (prf,cstr) + | Cstr(v,info) -> Inl ((prf,cstr,v,info)::acc)) (Inl []) l + + + type oproof = (vector * cstr_compat * num) option + + let merge_proof (oleft:oproof) (prf,cstr,v,info) (oright:oproof) = + let (l,r) = info.bound in + + let keep p ob bd = + match ob , bd with + | None , None -> None + | None , Some b -> Some(prf,cstr,b) + | Some _ , None -> ob + | Some(prfl,cstrl,bl) , Some b -> if p bl b then Some(prf,cstr, b) else ob in + + let oleft = keep (<=/) oleft l in + let oright = keep (>=/) oright r in + (* Now, there might be a contradiction *) + match oleft , oright with + | None , _ | _ , None -> Inl (oleft,oright) + | Some(prfl,cstrl,l) , Some(prfr,cstrr,r) -> + if l <=/ r + then Inl (oleft,oright) + else (* There is a contradiction - it should show up by scaling up the vectors - any pivot should do*) + match cstrr.coeffs with + | [] -> Inr (add (prfl,Int 1) (prfr,Int 1), cstrr) (* this is wrong *) + | (v,_)::_ -> + match pivot v (prfl,cstrl) (prfr,cstrr) with + | None -> failwith "merge_proof : pivot is not possible" + | Some x -> Inr x + +let mk_proof hyps prf = + (* I am keeping list - I might have a proof for the left bound and a proof for the right bound. + If I perform aggressive elimination of redundancies, I expect the list to be of length at most 2. + For each proof list, all the vectors should be of the form a.v for different constants a. + *) + + let rec mk_proof prf = + match prf with + | Hyp i -> [ ([i, Int 1] , List.nth hyps i) ] + + | Elim(v,prf1,prf2) -> + let prfsl = mk_proof prf1 + and prfsr = mk_proof prf2 in + (* I take only the pairs for which the elimination is meaningfull *) + forall_pairs (pivot v) prfsl prfsr + | And(prf1,prf2) -> + let prfsl1 = mk_proof prf1 + and prfsl2 = mk_proof prf2 in + (* detect trivial redundancies and contradictions *) + match normalise_proofs (prfsl1@prfsl2) with + | Inr x -> [x] (* This is a contradiction - this should be the end of the proof *) + | Inl l -> (* All the vectors are the same *) + let prfs = + List.fold_left (fun acc e -> + match acc with + | Inr _ -> acc (* I have a contradiction *) + | Inl (oleft,oright) -> merge_proof oleft e oright) (Inl(None,None)) l in + match prfs with + | Inr x -> [x] + | Inl (oleft,oright) -> + match oleft , oright with + | None , None -> [] + | None , Some(prf,cstr,_) | Some(prf,cstr,_) , None -> [prf,cstr] + | Some(prf1,cstr1,_) , Some(prf2,cstr2,_) -> [prf1,cstr1;prf2,cstr2] in + + mk_proof prf + + +end diff --git a/contrib/micromega/mutils.ml b/contrib/micromega/mutils.ml index 2473608f2..03e3999fd 100644 --- a/contrib/micromega/mutils.ml +++ b/contrib/micromega/mutils.ml @@ -17,6 +17,15 @@ let debug = false let fst' (Micromega.Pair(x,y)) = x let snd' (Micromega.Pair(x,y)) = y +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 @@ -24,6 +33,37 @@ 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 mapi f l = + let rec xmap i l = + match l with + | [] -> [] + | e::l -> (f i e)::xmap (i+1) l in + xmap 0 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) + | _ -> raise (Invalid_argument "map3") + + +let rec is_sublist l1 l2 = + match l1 ,l2 with + | [] ,_ -> true + | e::l1', [] -> false + | e::l1' , e'::l2' -> + if e = e' then is_sublist l1' l2' + else is_sublist l1 l2' + + + let list_try_find f = let rec try_find_f = function | [] -> failwith "try_find" diff --git a/contrib/micromega/vector.ml b/contrib/micromega/vector.ml deleted file mode 100644 index fee4ebfc1..000000000 --- a/contrib/micromega/vector.ml +++ /dev/null @@ -1,674 +0,0 @@ -(************************************************************************) -(* v * The Coq Proof Assistant / The Coq Development Team *) -(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *) -(* \VV/ **************************************************************) -(* // * This file is distributed under the terms of the *) -(* * GNU Lesser General Public License Version 2.1 *) -(************************************************************************) -(* *) -(* Micromega: A reflexive tactic using the Positivstellensatz *) -(* *) -(* Frédéric Besson (Irisa/Inria) 2006-2008 *) -(* *) -(************************************************************************) - -open Num - -module type S = -sig - type t - - val fresh : t -> int - - val null : t - - val is_null : t -> bool - - val get : int -> t -> num - - val update : int -> (num -> num) -> t -> t - (* behaviour is undef if index < 0 -- might loop*) - - val set : int -> num -> t -> t - - (* - For efficiency... - - val get_update : int -> (num -> num) -> t -> num * t - *) - - val mul : num -> t -> t - - val uminus : t -> t - - val add : t -> t -> t - - val dotp : t -> t -> num - - val lin_comb : num -> t -> num -> t -> t - (* lin_comb n1 t1 n2 t2 = (n1 * t1) + (n2 * t2) *) - - val gcd : t -> Big_int.big_int - - val normalise : t -> num * t - - val hash : t -> int - - val compare : t -> t -> int - - type it - - val iterator : t -> it - val element : it -> (num*it) option - - val string : t -> string - - type status = Pos | Neg - - (* the result list is ordered by fst *) - val status : t -> (int * status) list - - val from_list : num list -> t - val to_list : t -> num list - -end - - -module type SystemS = -sig - - module Vect : S - - module Cstr : - sig - type kind = Eq | Ge - val string_of_kind : kind -> string - type cstr = {coeffs : Vect.t ; op : kind ; cst : num} - val string_of_cstr : cstr -> string - val compare : cstr -> cstr -> int - end - open Cstr - - - module CstrBag : - sig - type t - exception Contradiction - - val empty : t - - val is_empty : t -> bool - - val add : cstr -> t -> t - (* c can be deduced from add c t *) - - val find : (cstr -> bool) -> t -> cstr option - - val fold : (cstr -> 'a -> 'a) -> t -> 'a -> 'a - - val status : t -> (int * (int list * int list)) list - (* aggregate of vector statuses *) - - val remove : cstr -> t -> t - - (* remove_list the ith element -- it is the ith element visited by 'fold' *) - - val split : (cstr -> int) -> t -> (int -> t) - - type it - val iterator : t -> it - val element : it -> (cstr*it) option - - end - -end - -let zero_num = Int 0 -let unit_num = Int 1 - - - - -module Cstr(V:S) = -struct - type kind = Eq | Ge - let string_of_kind = function Eq -> "Eq" | Ge -> "Ge" - - type cstr = {coeffs : V.t ; op : kind ; cst : num} - - let string_of_cstr {coeffs =a ; op = b ; cst =c} = - Printf.sprintf "{coeffs = %s;op=%s;cst=%s}" (V.string a) (string_of_kind b) (string_of_num c) - - type t = cstr - let compare - {coeffs = v1 ; op = op1 ; cst = c1} - {coeffs = v2 ; op = op2 ; cst = c2} = - Mutils.Cmp.compare_lexical [ - (fun () -> V.compare v1 v2); - (fun () -> Pervasives.compare op1 op2); - (fun () -> compare_num c1 c2) - ] - - -end - - - -module VList : S with type t = num list = -struct - type t = num list - - let fresh l = failwith "not implemented" - - let null = [] - - let is_null = List.for_all ((=/) zero_num) - - let normalise l = failwith "Not implemented" - (*match l with (* Buggy : What if the first num is zero! *) - | [] -> (Int 0,[]) - | [n] -> (n,[Int 1]) - | n::l -> (n, (Int 1) :: List.map (fun x -> x // n) l) - *) - - - let get i l = try List.nth l i with _ -> zero_num - - (* This is not tail-recursive *) - let rec update i f t = - match t with - | [] -> if i = 0 then [f zero_num] else (zero_num)::(update (i-1) f []) - | e::t -> if i = 0 then (f e)::t else e::(update (i-1) f t) - - let rec set i n t = - match t with - | [] -> if i = 0 then [n] else (zero_num)::(set (i-1) n []) - | e::t -> if i = 0 then (n)::t else e::(set (i-1) n t) - - - - - let rec mul z t = - match z with - | Int 0 -> null - | Int 1 -> t - | _ -> List.map (mult_num z) t - - let uminus t = mul (Int (-1)) t - - let rec add t1 t2 = - match t1,t2 with - | [], _ -> t2 - | _ , [] -> t1 - | e1::t1,e2::t2 -> (e1 +/ e2 )::(add t1 t2) - - let dotp t1 t2 = - let rec _dotp t1 t2 acc = - match t1, t2 with - | [] , _ -> acc - | _ , [] -> acc - | e1::t1,e2::t2 -> _dotp t1 t2 (acc +/ (e1 */ e2)) in - _dotp t1 t2 zero_num - - let add_mul n t1 t2 = - match n with - | Int 0 -> t2 - | Int 1 -> add t1 t2 - | _ -> - let rec _add_mul t1 t2 = - match t1,t2 with - | [], _ -> t2 - | _ , [] -> mul n t1 - | e1::t1,e2::t2 -> ( (n */e1) +/ e2 )::(_add_mul t1 t2) in - _add_mul t1 t2 - - let lin_comb n1 t1 n2 t2 = - match n1,n2 with - | Int 0 , _ -> mul n2 t2 - | Int 1 , _ -> add_mul n2 t2 t1 - | _ , Int 0 -> mul n1 t1 - | _ , Int 1 -> add_mul n1 t1 t2 - | _ -> - let rec _lin_comb t1 t2 = - match t1,t2 with - | [], _ -> mul n2 t2 - | _ , [] -> mul n1 t1 - | e1::t1,e2::t2 -> ( (n1 */e1) +/ (n2 */ e2 ))::(_lin_comb t1 t2) in - _lin_comb t1 t2 - - (* could be computed on the fly *) - let gcd t =Mutils.gcd_list t - - - - - let hash = Mutils.Cmp.hash_list int_of_num - - let compare = Mutils.Cmp.compare_list compare_num - - type it = t - let iterator (x:t) : it = x - let element it = - match it with - | [] -> None - | e::l -> Some (e,l) - - (* TODO: Buffer! *) - let string l = List.fold_right (fun n s -> (string_of_num n)^";"^s) l "" - - type status = Pos | Neg - - let status l = - let rec xstatus i l = - match l with - | [] -> [] - | e::l -> - begin - match compare_num e (Int 0) with - | 1 -> (i,Pos):: (xstatus (i+1) l) - | 0 -> xstatus (i+1) l - | -1 -> (i,Neg) :: (xstatus (i+1) l) - | _ -> assert false - end in - xstatus 0 l - - let from_list l = l - let to_list l = l - -end - -module VMap : S = -struct - module Map = Map.Make(struct type t = int let compare (x:int) (y:int) = Pervasives.compare x y end) - - type t = num Map.t - - let null = Map.empty - - let fresh m = failwith "not implemented" - - let is_null = Map.is_empty - - let normalise m = failwith "Not implemented" - - - - let get i l = try Map.find i l with _ -> zero_num - - let update i f t = - try - let res = f (Map.find i t) in - if res =/ zero_num - then Map.remove i t - else Map.add i res t - with - Not_found -> - let res = f zero_num in - if res =/ zero_num then t else Map.add i res t - - let set i n t = - if n =/ zero_num then Map.remove i t - else Map.add i n t - - - let rec mul z t = - match z with - | Int 0 -> null - | Int 1 -> t - | _ -> Map.map (mult_num z) t - - let uminus t = mul (Int (-1)) t - - - let map2 f m1 m2 = - let res,m2' = - Map.fold (fun k e (res,m2) -> - let v = f e (get k m2) in - if v =/ zero_num - then (res,Map.remove k m2) - else (Map.add k v res,Map.remove k m2)) m1 (Map.empty,m2) in - Map.fold (fun k e res -> - let v = f zero_num e in - if v =/ zero_num - then res else Map.add k v res) m2' res - - let add t1 t2 = map2 (+/) t1 t2 - - - let dotp t1 t2 = - Map.fold (fun k e res -> - res +/ (e */ get k t2)) t1 zero_num - - - - let add_mul n t1 t2 = - match n with - | Int 0 -> t2 - | Int 1 -> add t1 t2 - | _ -> map2 (fun x y -> (n */ x) +/ y) t1 t2 - - let lin_comb n1 t1 n2 t2 = - match n1,n2 with - | Int 0 , _ -> mul n2 t2 - | Int 1 , _ -> add_mul n2 t2 t1 - | _ , Int 0 -> mul n1 t1 - | _ , Int 1 -> add_mul n1 t1 t2 - | _ -> map2 (fun x y -> (n1 */ x) +/ (n2 */ y)) t1 t2 - - - let hash map = Map.fold (fun k e res -> k lxor (int_of_num e) lxor res) map 0 - - let compare = Map.compare compare_num - - type it = t * int - - let iterator (x:t) : it = (x,0) - - let element (mp,id) = - try - Some (Map.find id mp, (mp, id+1)) - with - Not_found -> None - - (* TODO: Buffer! *) - type status = Pos | Neg - - let status l = Map.fold (fun k e l -> - match compare_num e (Int 0) with - | 1 -> (k,Pos)::l - | 0 -> l - | -1 -> (k,Neg) :: l - | _ -> assert false) l [] - let from_list l = - let rec from_list i l map = - match l with - | [] -> map - | e::l -> from_list (i+1) l (if e <>/ Int 0 then Map.add i e map else map) in - from_list 0 l Map.empty - - let gcd m = - let res = Map.fold (fun _ e x -> Big_int.gcd_big_int x (Mutils.numerator e)) m Big_int.zero_big_int in - if Big_int.compare_big_int res Big_int.zero_big_int = 0 - then Big_int.unit_big_int else res - - - let to_list m = - let l = List.rev (Map.fold (fun k e l -> (k,e)::l) m []) in - let rec xto_list i l = - match l with - | [] -> [] - | (x,v)::l' -> if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in - xto_list 0 l - - let string l = VList.string (to_list l) - - -end - - -module VSparse : S = -struct - - type t = (int*num) list - - let null = [] - - let fresh l = List.fold_left (fun acc (i,_) -> max (i+1) acc) 0 l - - let is_null l = l = [] - - let rec is_sorted l = - match l with - | [] -> true - | [e] -> true - | (i,_)::(j,x)::l -> i < j && is_sorted ((j,x)::l) - - - let check l = (List.for_all (fun (_,n) -> compare_num n (Int 0) <> 0) l) && (is_sorted l) - - (* let get i t = - assert (check t); - try List.assoc i t with Not_found -> zero_num *) - - let rec get (i:int) t = - match t with - | [] -> zero_num - | (j,n)::t -> - match compare i j with - | 0 -> n - | 1 -> get i t - | _ -> zero_num - - let cons i v rst = if v =/ Int 0 then rst else (i,v)::rst - - let rec update i f t = - match t with - | [] -> cons i (f zero_num) [] - | (k,v)::l -> - match Pervasives.compare i k with - | 0 -> cons k (f v) l - | -1 -> cons i (f zero_num) t - | 1 -> (k,v) ::(update i f l) - | _ -> failwith "compare_num" - - let update i f t = - assert (check t); - let res = update i f t in - assert (check t) ; res - - - let rec set i n t = - match t with - | [] -> cons i n [] - | (k,v)::l -> - match Pervasives.compare i k with - | 0 -> cons k n l - | -1 -> cons i n t - | 1 -> (k,v) :: (set i n l) - | _ -> failwith "compare_num" - - - let rec map f l = - match l with - | [] -> [] - | (i,e)::l -> cons i (f e) (map f l) - - let rec mul z t = - match z with - | Int 0 -> null - | Int 1 -> t - | _ -> List.map (fun (i,n) -> (i, mult_num z n)) t - - let mul z t = - assert (check t) ; - let res = mul z t in - assert (check res) ; - res - - let uminus t = mul (Int (-1)) t - - - let normalise l = - match l with - | [] -> (Int 0,[]) - | (i,n)::_ -> (n, mul ((Int 1) // n) l) - - - let rec map2 f m1 m2 = - match m1, m2 with - | [] , [] -> [] - | l , [] -> map (fun x -> f x zero_num) l - | [] ,l -> map (f zero_num) l - | (i,e)::l1,(i',e')::l2 -> - match Pervasives.compare i i' with - | 0 -> cons i (f e e') (map2 f l1 l2) - | -1 -> cons i (f e zero_num) (map2 f l1 m2) - | 1 -> cons i' (f zero_num e') (map2 f m1 l2) - | _ -> assert false - - (* let add t1 t2 = map2 (+/) t1 t2*) - - let rec add (m1:t) (m2:t) = - match m1, m2 with - | [] , [] -> [] - | l , [] -> l - | [] ,l -> l - | (i,e)::l1,(i',e')::l2 -> - match Pervasives.compare i i' with - | 0 -> cons i ( e +/ e') (add l1 l2) - | -1 -> (i,e) :: (add l1 m2) - | 1 -> (i', e') :: (add m1 l2) - | _ -> assert false - - - - - let add t1 t2 = - assert (check t1 && check t2); - let res = add t1 t2 in - assert (check res); - res - - - let rec dotp (t1:t) (t2:t) = - match t1, t2 with - | [] , _ -> zero_num - | _ , [] -> zero_num - | (i,e)::l1 , (i',e')::l2 -> - match Pervasives.compare i i' with - | 0 -> (e */ e') +/ (dotp l1 l2) - | -1 -> dotp l1 t2 - | 1 -> dotp t1 l2 - | _ -> assert false - - let dotp t1 t2 = - assert (check t1 && check t2) ; dotp t1 t2 - - let add_mul n t1 t2 = - match n with - | Int 0 -> t2 - | Int 1 -> add t1 t2 - | _ -> map2 (fun x y -> (n */ x) +/ y) t1 t2 - - let add_mul n (t1:t) (t2:t) = - match n with - | Int 0 -> t2 - | Int 1 -> add t1 t2 - | _ -> - let rec xadd_mul m1 m2 = - match m1, m2 with - | [] , [] -> [] - | _ , [] -> mul n m1 - | [] , _ -> m2 - | (i,e)::l1,(i',e')::l2 -> - match Pervasives.compare i i' with - | 0 -> cons i ( n */ e +/ e') (xadd_mul l1 l2) - | -1 -> (i,n */ e) :: (xadd_mul l1 m2) - | 1 -> (i', e') :: (xadd_mul m1 l2) - | _ -> assert false in - xadd_mul t1 t2 - - - - - let lin_comb n1 t1 n2 t2 = - match n1,n2 with - | Int 0 , _ -> mul n2 t2 - | Int 1 , _ -> add_mul n2 t2 t1 - | _ , Int 0 -> mul n1 t1 - | _ , Int 1 -> add_mul n1 t1 t2 - | _ -> (*map2 (fun x y -> (n1 */ x) +/ (n2 */ y)) t1 t2*) - let rec xlin_comb m1 m2 = - match m1, m2 with - | [] , [] -> [] - | _ , [] -> mul n1 m1 - | [] , _ -> mul n2 m2 - | (i,e)::l1,(i',e')::l2 -> - match Pervasives.compare i i' with - | 0 -> cons i ( n1 */ e +/ n2 */ e') (xlin_comb l1 l2) - | -1 -> (i,n1 */ e) :: (xlin_comb l1 m2) - | 1 -> (i', n2 */ e') :: (xlin_comb m1 l2) - | _ -> assert false in - xlin_comb t1 t2 - - - - - - let lin_comb n1 t1 n2 t2 = - assert (check t1 && check t2); - let res = lin_comb n1 t1 n2 t2 in - assert (check res); res - - let hash = Mutils.Cmp.hash_list (fun (x,y) -> (Hashtbl.hash x) lxor (int_of_num y)) - - - let compare = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical - [ - (fun () -> Pervasives.compare (fst x) (fst y)); - (fun () -> compare_num (snd x) (snd y))]) - - (* - let compare (x:t) (y:t) = - let rec xcompare acc1 acc2 x y = - match x , y with - | [] , [] -> xcomp acc1 acc2 - | [] , _ -> -1 - | _ , [] -> 1 - | (i,n1)::l1 , (j,n2)::l2 -> - match Pervasives.compare i j with - | 0 -> xcompare (n1::acc1) (n2::acc2) l1 l2 - | c -> c - and xcomp acc1 acc2 = Mutils.Cmp.compare_list compare_num acc1 acc2 in - xcompare [] [] x y - *) - - type it = t - - let iterator (x:t) : it = x - - let element l = failwith "Not_implemented" - - (* TODO: Buffer! *) - type status = Pos | Neg - - let status l = List.map (fun (i,e) -> - match compare_num e (Int 0) with - | 1 -> i,Pos - | -1 -> i,Neg - | _ -> assert false) l - - let from_list (l: num list) = - let rec xfrom_list i l = - match l with - | [] -> [] - | e::l -> - if e <>/ Int 0 - then (i,e)::(xfrom_list (i+1) l) - else xfrom_list (i+1) l in - - let res = xfrom_list 0 l in - assert (check res) ; res - - - let gcd m = - let res = List.fold_left (fun x (i,e) -> Big_int.gcd_big_int x (Mutils.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 to_list m = - let rec xto_list i l = - match l with - | [] -> [] - | (x,v)::l' -> - if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in - xto_list 0 m - - let to_list l = - assert (check l); - to_list l - - - let string l = VList.string (to_list l) - -end diff --git a/test-suite/micromega/zomicron.v b/test-suite/micromega/zomicron.v index 2b40f6c93..3087db1ad 100644 --- a/test-suite/micromega/zomicron.v +++ b/test-suite/micromega/zomicron.v @@ -20,8 +20,9 @@ Proof. lia. Qed. -Lemma omega_nightmare : forall x y, 27 <= 11 * x + 13 * y <= 45 -> 7 * x - 9 * y = 4 -> -10 <= 7 * x - 9 * y <= 4 -> False. +Lemma omega_nightmare : forall x y, 27 <= 11 * x + 13 * y <= 45 -> -10 <= 7 * x - 9 * y <= 4 -> False. Proof. intros ; intuition auto. lia. Qed. + |