diff options
author | Stephane Glondu <steph@glondu.net> | 2012-01-12 16:02:20 +0100 |
---|---|---|
committer | Stephane Glondu <steph@glondu.net> | 2012-01-12 16:02:20 +0100 |
commit | 97fefe1fcca363a1317e066e7f4b99b9c1e9987b (patch) | |
tree | 97ec6b7d831cc5fb66328b0c63a11db1cbb2f158 /plugins/micromega/certificate.ml | |
parent | 300293c119981054c95182a90c829058530a6b6f (diff) |
Imported Upstream version 8.4~betaupstream/8.4_beta
Diffstat (limited to 'plugins/micromega/certificate.ml')
-rw-r--r-- | plugins/micromega/certificate.ml | 1244 |
1 files changed, 829 insertions, 415 deletions
diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml index bcab73ec..540d1b9c 100644 --- a/plugins/micromega/certificate.ml +++ b/plugins/micromega/certificate.ml @@ -1,6 +1,6 @@ (************************************************************************) (* v * The Coq Proof Assistant / The Coq Development Team *) -(* <O___,, * INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2011 *) +(* <O___,, * INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010 *) (* \VV/ **************************************************************) (* // * This file is distributed under the terms of the *) (* * GNU Lesser General Public License Version 2.1 *) @@ -15,153 +15,18 @@ (* We take as input a list of polynomials [p1...pn] and return an unfeasibility certificate polynomial. *) -(*open Micromega.Polynomial*) +type var = int + + + open Big_int open Num -open Sos_lib +open Polynomial module Mc = Micromega module Ml2C = Mutils.CamlToCoq module C2Ml = Mutils.CoqToCaml -let (<+>) = add_num -let (<->) = minus_num -let (<*>) = mult_num - -type var = Mc.positive - -module Monomial : -sig - type t - val const : t - val var : var -> t - val find : var -> t -> int - val mult : var -> t -> t - val prod : t -> t -> t - val compare : t -> t -> int - val pp : out_channel -> t -> unit - val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a -end - = -struct - (* A monomial is represented by a multiset of variables *) - module Map = Map.Make(struct type t = var let compare = Pervasives.compare end) - open Map - - type t = int Map.t - - (* The monomial that corresponds to a constant *) - let const = Map.empty - - (* The monomial 'x' *) - let var x = Map.add x 1 Map.empty - - (* Get the degre of a variable in a monomial *) - let find x m = try find x m with Not_found -> 0 - - (* Multiply a monomial by a variable *) - let mult x m = add x ( (find x m) + 1) m - - (* Product of monomials *) - let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2 - - (* Total ordering of monomials *) - let compare m1 m2 = Map.compare Pervasives.compare m1 m2 - - let pp o m = Map.iter (fun k v -> - if v = 1 then Printf.fprintf o "x%i." (C2Ml.index k) - else Printf.fprintf o "x%i^%i." (C2Ml.index k) v) m - - let fold = fold - -end - - -module Poly : - (* A polynomial is a map of monomials *) - (* - This is probably a naive implementation - (expected to be fast enough - Coq is probably the bottleneck) - *The new ring contribution is using a sparse Horner representation. - *) -sig - type t - val get : Monomial.t -> t -> num - val variable : var -> t - val add : Monomial.t -> num -> t -> t - val constant : num -> t - val mult : Monomial.t -> num -> t -> t - val product : t -> t -> t - val addition : t -> t -> t - val uminus : t -> t - val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a - val pp : out_channel -> t -> unit - val compare : t -> t -> int - val is_null : t -> bool -end = -struct - (*normalisation bug : 0*x ... *) - module P = Map.Make(Monomial) - open P - - type t = num P.t - - let pp o p = P.iter (fun k v -> - if compare_num v (Int 0) <> 0 - then - if Monomial.compare Monomial.const k = 0 - then Printf.fprintf o "%s " (string_of_num v) - else Printf.fprintf o "%s*%a " (string_of_num v) Monomial.pp k) p - - (* Get the coefficient of monomial mn *) - let get : Monomial.t -> t -> num = - fun mn p -> try find mn p with Not_found -> (Int 0) - - - (* The polynomial 1.x *) - let variable : var -> t = - fun x -> add (Monomial.var x) (Int 1) empty - - (*The constant polynomial *) - let constant : num -> t = - fun c -> add (Monomial.const) c empty - - (* The addition of a monomial *) - - let add : Monomial.t -> num -> t -> t = - fun mn v p -> - let vl = (get mn p) <+> v in - add mn vl p - - - (** Design choice: empty is not a polynomial - I do not remember why .... - **) - - (* The product by a monomial *) - let mult : Monomial.t -> num -> t -> t = - fun mn v p -> - fold (fun mn' v' res -> P.add (Monomial.prod mn mn') (v<*>v') res) p empty - - - let addition : t -> t -> t = - fun p1 p2 -> fold (fun mn v p -> add mn v p) p1 p2 - - - let product : t -> t -> t = - fun p1 p2 -> - fold (fun mn v res -> addition (mult mn v p2) res ) p1 empty - - - let uminus : t -> t = - fun p -> map (fun v -> minus_num v) p - - let fold = P.fold - - let is_null p = fold (fun mn vl b -> b & sign_num vl = 0) p true - - let compare = compare compare_num -end open Mutils type 'a number_spec = { @@ -178,10 +43,10 @@ let z_spec = { number_to_num = (fun x -> Big_int (C2Ml.z_big_int x)); zero = Mc.Z0; unit = Mc.Zpos Mc.XH; - mult = Mc.zmult; + mult = Mc.Z.mul; eqb = Mc.zeq_bool } - + let q_spec = { bigint_to_number = (fun x -> {Mc.qnum = Ml2C.bigint x; Mc.qden = Mc.XH}); @@ -195,56 +60,58 @@ let q_spec = { let r_spec = z_spec - - let dev_form n_spec p = - let rec dev_form p = + let rec dev_form p = match p with | Mc.PEc z -> Poly.constant (n_spec.number_to_num z) - | Mc.PEX v -> Poly.variable v - | Mc.PEmul(p1,p2) -> + | Mc.PEX v -> Poly.variable (C2Ml.positive v) + | Mc.PEmul(p1,p2) -> let p1 = dev_form p1 in let p2 = dev_form p2 in - Poly.product p1 p2 + Poly.product p1 p2 | Mc.PEadd(p1,p2) -> Poly.addition (dev_form p1) (dev_form p2) | Mc.PEopp p -> Poly.uminus (dev_form p) | Mc.PEsub(p1,p2) -> Poly.addition (dev_form p1) (Poly.uminus (dev_form p2)) - | Mc.PEpow(p,n) -> + | Mc.PEpow(p,n) -> let p = dev_form p in let n = C2Ml.n n in - let rec pow n = - if n = 0 + let rec pow n = + if n = 0 then Poly.constant (n_spec.number_to_num n_spec.unit) else Poly.product p (pow (n-1)) in pow n in dev_form p -let monomial_to_polynomial mn = - Monomial.fold - (fun v i acc -> - let mn = if i = 1 then Mc.PEX v else Mc.PEpow (Mc.PEX v ,Ml2C.n i) in - if acc = Mc.PEc (Mc.Zpos Mc.XH) - then mn - else Mc.PEmul(mn,acc)) - mn - (Mc.PEc (Mc.Zpos Mc.XH)) +let monomial_to_polynomial mn = + Monomial.fold + (fun v i acc -> + let v = Ml2C.positive v in + let mn = if i = 1 then Mc.PEX v else Mc.PEpow (Mc.PEX v ,Ml2C.n i) in + if acc = Mc.PEc (Mc.Zpos Mc.XH) + then mn + else Mc.PEmul(mn,acc)) + mn + (Mc.PEc (Mc.Zpos Mc.XH)) + -let list_to_polynomial vars l = + +let list_to_polynomial vars l = assert (List.for_all (fun x -> ceiling_num x =/ x) l); let var x = monomial_to_polynomial (List.nth vars x) in + let rec xtopoly p i = function | [] -> p - | c::l -> if c =/ (Int 0) then xtopoly p (i+1) l + | c::l -> if c =/ (Int 0) then xtopoly p (i+1) l else let c = Mc.PEc (Ml2C.bigint (numerator c)) in - let mn = + let mn = if c = Mc.PEc (Mc.Zpos Mc.XH) then var i else Mc.PEmul (c,var i) in let p' = if p = Mc.PEc Mc.Z0 then mn else Mc.PEadd (mn, p) in xtopoly p' (i+1) l in - + xtopoly (Mc.PEc Mc.Z0) 0 l let rec fixpoint f x = @@ -252,61 +119,54 @@ let rec fixpoint f x = if y' = x then y' else fixpoint f y' - - - - - - - -let rec_simpl_cone n_spec e = - let simpl_cone = +let rec_simpl_cone n_spec e = + let simpl_cone = Mc.simpl_cone n_spec.zero n_spec.unit n_spec.mult n_spec.eqb in let rec rec_simpl_cone = function - | Mc.PsatzMulE(t1, t2) -> + | Mc.PsatzMulE(t1, t2) -> simpl_cone (Mc.PsatzMulE (rec_simpl_cone t1, rec_simpl_cone t2)) - | Mc.PsatzAdd(t1,t2) -> + | Mc.PsatzAdd(t1,t2) -> simpl_cone (Mc.PsatzAdd (rec_simpl_cone t1, rec_simpl_cone t2)) | x -> simpl_cone x in rec_simpl_cone e - - + + let simplify_cone n_spec c = fixpoint (rec_simpl_cone n_spec) c - -type cone_prod = - Const of cone - | Ideal of cone *cone - | Mult of cone * cone + +type cone_prod = + Const of cone + | Ideal of cone *cone + | Mult of cone * cone | Other of cone and cone = Mc.zWitness let factorise_linear_cone c = - - let rec cone_list c l = + + let rec cone_list c l = match c with | Mc.PsatzAdd (x,r) -> cone_list r (x::l) | _ -> c :: l in - + let factorise c1 c2 = match c1 , c2 with - | Mc.PsatzMulC(x,y) , Mc.PsatzMulC(x',y') -> + | Mc.PsatzMulC(x,y) , Mc.PsatzMulC(x',y') -> if x = x' then Some (Mc.PsatzMulC(x, Mc.PsatzAdd(y,y'))) else None - | Mc.PsatzMulE(x,y) , Mc.PsatzMulE(x',y') -> + | Mc.PsatzMulE(x,y) , Mc.PsatzMulE(x',y') -> if x = x' then Some (Mc.PsatzMulE(x, Mc.PsatzAdd(y,y'))) else None | _ -> None in - + let rec rebuild_cone l pending = match l with | [] -> (match pending with | None -> Mc.PsatzZ | Some p -> p ) - | e::l -> + | e::l -> (match pending with - | None -> rebuild_cone l (Some e) + | None -> rebuild_cone l (Some e) | Some p -> (match factorise p e with | None -> Mc.PsatzAdd(p, rebuild_cone l (Some e)) | Some f -> rebuild_cone l (Some f) ) @@ -316,15 +176,15 @@ let factorise_linear_cone c = -(* The binding with Fourier might be a bit obsolete +(* The binding with Fourier might be a bit obsolete -- how does it handle equalities ? *) (* Certificates are elements of the cone such that P = 0 *) (* To begin with, we search for certificates of the form: - a1.p1 + ... an.pn + b1.q1 +... + bn.qn + c = 0 + a1.p1 + ... an.pn + b1.q1 +... + bn.qn + c = 0 where pi >= 0 qi > 0 - ai >= 0 + ai >= 0 bi >= 0 Sum bi + c >= 1 This is a linear problem: each monomial is considered as a variable. @@ -334,216 +194,209 @@ 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 Vect = Fourier.Vect*) -(*open Fourier.Cstr*) (* fold_left followed by a rev ! *) -let constrain_monomial mn l = +let constrain_monomial mn l = let coeffs = List.fold_left (fun acc p -> (Poly.get mn p)::acc) [] l in if mn = Monomial.const - then - { coeffs = Vect.from_list ((Big_int unit_big_int):: (List.rev coeffs)) ; - op = Eq ; + then + { coeffs = Vect.from_list ((Big_int unit_big_int):: (List.rev coeffs)) ; + op = Eq ; cst = Big_int zero_big_int } else - { coeffs = Vect.from_list ((Big_int zero_big_int):: (List.rev coeffs)) ; - op = Eq ; + { coeffs = Vect.from_list ((Big_int zero_big_int):: (List.rev coeffs)) ; + op = Eq ; cst = Big_int zero_big_int } - -let positivity l = - let rec xpositivity i l = + +let positivity l = + let rec xpositivity i l = match l with | [] -> [] | (_,Mc.Equal)::l -> xpositivity (i+1) l - | (_,_)::l -> - {coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ; - op = Ge ; + | (_,_)::l -> + {coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ; + op = Ge ; cst = Int 0 } :: (xpositivity (i+1) l) in xpositivity 0 l let string_of_op = function - | Mc.Strict -> "> 0" - | Mc.NonStrict -> ">= 0" + | Mc.Strict -> "> 0" + | Mc.NonStrict -> ">= 0" | Mc.Equal -> "= 0" | Mc.NonEqual -> "<> 0" -(* If the certificate includes at least one strict inequality, +(* If the certificate includes at least one strict inequality, the obtained polynomial can also be 0 *) let build_linear_system l = - (* Gather the monomials: HINT add up of the polynomials *) + (* Gather the monomials: HINT add up of the polynomials ==> This does not work anymore *) let l' = List.map fst l in - let monomials = - List.fold_left (fun acc p -> Poly.addition p acc) (Poly.constant (Int 0)) l' + + let module MonSet = Set.Make(Monomial) in + + let monomials = + List.fold_left (fun acc p -> + Poly.fold (fun m _ acc -> MonSet.add m acc) p acc) + (MonSet.singleton Monomial.const) l' in (* For each monomial, compute a constraint *) - let s0 = - Poly.fold (fun mn _ res -> (constrain_monomial mn l')::res) monomials [] in + let s0 = + MonSet.fold (fun mn res -> (constrain_monomial mn l')::res) monomials [] in (* I need at least something strictly positive *) let strict = { coeffs = Vect.from_list ((Big_int unit_big_int):: - (List.map (fun (x,y) -> - match y with Mc.Strict -> - Big_int unit_big_int + (List.map (fun (x,y) -> + match y with Mc.Strict -> + Big_int unit_big_int | _ -> Big_int zero_big_int) l)); op = Ge ; cst = Big_int unit_big_int } in (* Add the positivity constraint *) - {coeffs = Vect.from_list ([Big_int unit_big_int]) ; - op = Ge ; + {coeffs = Vect.from_list ([Big_int unit_big_int]) ; + op = Ge ; cst = Big_int zero_big_int}::(strict::(positivity l)@s0) let big_int_to_z = Ml2C.bigint - -(* For Q, this is a pity that the certificate has been scaled + +(* For Q, this is a pity that the certificate has been scaled -- at a lower layer, certificates are using nums... *) -let make_certificate n_spec (cert,li) = +let make_certificate n_spec (cert,li) = let bint_to_cst = n_spec.bigint_to_number in match cert with | [] -> failwith "empty_certificate" - | e::cert' -> - let cst = match compare_big_int e zero_big_int with + | e::cert' -> +(* let cst = match compare_big_int e zero_big_int with | 0 -> Mc.PsatzZ - | 1 -> Mc.PsatzC (bint_to_cst e) - | _ -> failwith "positivity error" - in + | 1 -> Mc.PsatzC (bint_to_cst e) + | _ -> failwith "positivity error" + in *) let rec scalar_product cert l = match cert with | [] -> Mc.PsatzZ - | c::cert -> match l with - | [] -> failwith "make_certificate(1)" - | i::l -> - let r = scalar_product cert l in - match compare_big_int c zero_big_int with - | -1 -> Mc.PsatzAdd ( - Mc.PsatzMulC (Mc.Pc ( bint_to_cst c), Mc.PsatzIn (Ml2C.nat i)), - r) - | 0 -> r - | _ -> Mc.PsatzAdd ( - Mc.PsatzMulE (Mc.PsatzC (bint_to_cst c), Mc.PsatzIn (Ml2C.nat i)), - r) in - - ((factorise_linear_cone - (simplify_cone n_spec (Mc.PsatzAdd (cst, scalar_product cert' li))))) + | c::cert -> + match l with + | [] -> failwith "make_certificate(1)" + | i::l -> + let r = scalar_product cert l in + match compare_big_int c zero_big_int with + | -1 -> Mc.PsatzAdd ( + Mc.PsatzMulC (Mc.Pc ( bint_to_cst c), Mc.PsatzIn (Ml2C.nat i)), + r) + | 0 -> r + | _ -> Mc.PsatzAdd ( + Mc.PsatzMulE (Mc.PsatzC (bint_to_cst c), Mc.PsatzIn (Ml2C.nat i)), + r) in + (factorise_linear_cone + (simplify_cone n_spec (scalar_product cert' li))) exception Found of Monomial.t exception Strict -let primal l = +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 + Poly.fold (fun mn vl (map,vect) -> + if mn = Monomial.const then (map,vect) - else + 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 (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) = +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 | Inr _ -> None - | Inl cert -> Some (rats_to_ints (Vect.to_list cert)) + | Inl cert -> Some (rats_to_ints (Vect.to_list cert)) (* should not use rats_to_ints *) - with x -> - if debug - then (Printf.printf "raw certificate %s" (Printexc.to_string x); + with x -> + if debug + then (Printf.printf "raw certificate %s" (Printexc.to_string x); flush stdout) ; None -let raw_certificate l = - try +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 ; + | 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 ; + if debug then Printf.printf "CProof : %a" Vect.pp_vect cert ; Some (rats_to_ints (Vect.to_list cert)) | Inl _ -> None - with Strict -> + with Strict -> (* Fourier elimination should handle > *) - dual_raw_certificate l + dual_raw_certificate l -let simple_linear_prover (*to_constant*) l = +let simple_linear_prover l = let (lc,li) = List.split l in match raw_certificate lc with | None -> None (* No certificate *) - | Some cert -> (* make_certificate to_constant*)Some (cert,li) + | Some cert -> Some (cert,li) + -let linear_prover n_spec l = - let li = List.combine l (interval 0 (List.length l -1)) in - let (l1,l') = List.partition - (fun (x,_) -> if snd x = Mc.NonEqual then true else false) li in - let l' = List.map - (fun ((x,y),i) -> match y with - Mc.NonEqual -> failwith "cannot happen" - | y -> ((dev_form n_spec x, y),i)) l' in - simple_linear_prover (*n_spec*) l' +let linear_prover n_spec l = + let build_system n_spec l = + let li = List.combine l (interval 0 (List.length l -1)) in + let (l1,l') = List.partition + (fun (x,_) -> if snd x = Mc.NonEqual then true else false) li in + List.map + (fun ((x,y),i) -> match y with + Mc.NonEqual -> failwith "cannot happen" + | y -> ((dev_form n_spec x, y),i)) l' in + let l' = build_system n_spec l in + simple_linear_prover (*n_spec*) l' let linear_prover n_spec l = try linear_prover n_spec l with x -> (print_string (Printexc.to_string x); None) -let linear_prover_with_cert spec l = +let linear_prover_with_cert spec l = match linear_prover spec l with | None -> None | Some cert -> Some (make_certificate spec cert) -(* zprover.... *) - -(* I need to gather the set of variables ---> - Then go for fold - Once I have an interval, I need a certificate : 2 other fourier elims. - (I could probably get the certificate directly - as it is done in the fourier contrib.) -*) let make_linear_system l = let l' = List.map fst l in - let monomials = List.fold_left (fun acc p -> Poly.addition p acc) + let monomials = List.fold_left (fun acc p -> Poly.addition p acc) (Poly.constant (Int 0)) l' in - let monomials = Poly.fold + let monomials = Poly.fold (fun mn _ l -> if mn = Monomial.const then l else mn::l) monomials [] in - (List.map (fun (c,op) -> - {coeffs = Vect.from_list (List.map (fun mn -> (Poly.get mn c)) monomials) ; - op = op ; + (List.map (fun (c,op) -> + {coeffs = Vect.from_list (List.map (fun mn -> (Poly.get mn c)) monomials) ; + op = op ; cst = minus_num ( (Poly.get Monomial.const c))}) l ,monomials) @@ -552,130 +405,66 @@ let pplus x y = Mc.PEadd(x,y) let pmult x y = Mc.PEmul(x,y) let pconst x = Mc.PEc x let popp x = Mc.PEopp x - + let debug = false - + (* keep track of enumerated vectors *) -let rec mem p x l = +let rec mem p x l = match l with [] -> false | e::l -> if p x e then true else mem p x l -let rec remove_assoc p x l = +let rec remove_assoc p x l = match l with [] -> [] | e::l -> if p x (fst e) then - remove_assoc p x l else e::(remove_assoc p x l) + remove_assoc p x l else e::(remove_assoc p x l) let eq x y = Vect.compare x y = 0 let remove e l = List.fold_left (fun l x -> if eq x e then l else x::l) [] l -(* The prover is (probably) incomplete -- +(* The prover is (probably) incomplete -- only searching for naive cutting planes *) -let candidates sys = - let ll = List.fold_right ( - fun (e,k) r -> - match k with - | Mc.NonStrict -> (dev_form z_spec e , Ge)::r - | Mc.Equal -> (dev_form z_spec e , Eq)::r - (* we already know the bound -- don't compute it again *) - | _ -> failwith "Cannot happen candidates") sys [] in - - let (sys,var_mn) = make_linear_system ll in - let vars = mapi (fun _ i -> Vect.set i (Int 1) Vect.null) var_mn in - (List.fold_left (fun l cstr -> - let gcd = Big_int (Vect.gcd cstr.coeffs) in - if gcd =/ (Int 1) && cstr.op = Eq - then l - else (Vect.mul (Int 1 // gcd) cstr.coeffs)::l) [] sys) @ vars - - - - -let rec xzlinear_prover planes sys = - match linear_prover z_spec sys with - | Some prf -> Some (Mc.RatProof (make_certificate z_spec prf,Mc.DoneProof)) - | None -> (* find the candidate with the smallest range *) - (* Grrr - linear_prover is also calling 'make_linear_system' *) - let ll = List.fold_right (fun (e,k) r -> match k with - Mc.NonEqual -> r - | k -> (dev_form z_spec e , - match k with - Mc.NonStrict -> Ge - | Mc.Equal -> Eq - | Mc.Strict | Mc.NonEqual -> failwith "Cannot happen") :: r) sys [] in - let (ll,var) = make_linear_system ll in - let candidates = List.fold_left (fun acc vect -> - match Fourier.optimise vect ll with - | None -> acc - | Some i -> -(* Printf.printf "%s in %s\n" (Vect.string vect) (string_of_intrvl i) ; *) - flush stdout ; - (vect,i) ::acc) [] planes in - - let smallest_interval = - match List.fold_left (fun (x1,i1) (x2,i2) -> - if Itv.smaller_itv i1 i2 - then (x1,i1) else (x2,i2)) (Vect.null,(None,None)) candidates - with - | (x,(Some i, Some j)) -> Some(i,x,j) - | x -> None (* This might be a cutting plane *) - in - match smallest_interval with - | Some (lb,e,ub) -> - let (lbn,lbd) = - (Ml2C.bigint (sub_big_int (numerator lb) unit_big_int), - Ml2C.bigint (denominator lb)) in - let (ubn,ubd) = - (Ml2C.bigint (add_big_int unit_big_int (numerator ub)) , - Ml2C.bigint (denominator ub)) in - let expr = list_to_polynomial var (Vect.to_list e) in - (match - (*x <= ub -> x > ub *) - linear_prover z_spec - ((pplus (pmult (pconst ubd) expr) (popp (pconst ubn)), - Mc.NonStrict) :: sys), - (* lb <= x -> lb > x *) - linear_prover z_spec - ((pplus (popp (pmult (pconst lbd) expr)) (pconst lbn), - Mc.NonStrict)::sys) - with - | Some cub , Some clb -> - (match zlinear_enum (remove e planes) expr - (ceiling_num lb) (floor_num ub) sys - with - | None -> None - | Some prf -> - let bound_proof (c,l) = make_certificate z_spec (List.tl c , List.tl (List.map (fun x -> x -1) l)) in - - Some (Mc.EnumProof((*Ml2C.q lb,expr,Ml2C.q ub,*) bound_proof clb, bound_proof cub,prf))) - | _ -> None - ) - | _ -> None -and zlinear_enum planes expr clb cub l = - if clb >/ cub - then Some [] - else - let pexpr = pplus (popp (pconst (Ml2C.bigint (numerator clb)))) expr in - let sys' = (pexpr, Mc.Equal)::l in - (*let enum = *) - match xzlinear_prover planes sys' with - | None -> if debug then print_string "zlp?"; None - | Some prf -> if debug then print_string "zlp!"; - match zlinear_enum planes expr (clb +/ (Int 1)) cub l with - | None -> None - | Some prfl -> Some (prf :: prfl) +let develop_constraint z_spec (e,k) = + match k with + | Mc.NonStrict -> (dev_form z_spec e , Ge) + | Mc.Equal -> (dev_form z_spec e , Eq) + | _ -> assert false + + +let op_of_op_compat = function + | Ge -> Mc.NonStrict + | Eq -> Mc.Equal + + +let integer_vector coeffs = + let vars , coeffs = List.split coeffs in + List.combine vars (List.map (fun x -> Big_int x) (rats_to_ints coeffs)) + +let integer_cstr {coeffs = coeffs ; op = op ; cst = cst } = + let vars , coeffs = List.split coeffs in + match rats_to_ints (cst::coeffs) with + | cst :: coeffs -> + { + coeffs = List.combine vars (List.map (fun x -> Big_int x) coeffs) ; + op = op ; cst = Big_int cst} + | _ -> assert false + + +let pexpr_of_cstr_compat var cstr = + let {coeffs = coeffs ; op = op ; cst = cst } = integer_cstr cstr in + try + let expr = list_to_polynomial var (Vect.to_list coeffs) in + let d = Ml2C.bigint (denominator cst) in + let n = Ml2C.bigint (numerator cst) in + (pplus (pmult (pconst d) expr) (popp (pconst n)), op_of_op_compat op) + with Failure _ -> failwith "pexpr_of_cstr_compat" + + -let zlinear_prover sys = - let candidates = candidates sys in - (* Printf.printf "candidates %d" (List.length candidates) ; *) - (*let t0 = Sys.time () in*) - let res = xzlinear_prover candidates sys in - (*Printf.printf "Time prover : %f" (Sys.time () -. t0) ;*) res open Sos_types -open Mutils -let rec scale_term t = +let rec scale_term t = match t with | Zero -> unit_big_int , Zero | Const n -> (denominator n) , Const (Big_int (numerator n)) @@ -708,7 +497,7 @@ let get_index_of_ith_match f i l = match l with | [] -> failwith "bad index" | e::l -> if f e - then + then (if j = i then res else get (j+1) (res+1) l ) else get j (res+1) l in get 0 0 l @@ -722,19 +511,19 @@ let rec scale_certificate pos = match pos with | Rational_eq n -> (denominator n) , Rational_eq (Big_int (numerator n)) | Rational_le n -> (denominator n) , Rational_le (Big_int (numerator n)) | Rational_lt n -> (denominator n) , Rational_lt (Big_int (numerator n)) - | Square t -> let s,t' = scale_term t in + | Square t -> let s,t' = scale_term t in mult_big_int s s , Square t' | Eqmul (t, y) -> let s1,y1 = scale_term t and s2,y2 = scale_certificate y in mult_big_int s1 s2 , Eqmul (y1,y2) - | Sum (y, z) -> let s1,y1 = scale_certificate y + | Sum (y, z) -> let s1,y1 = scale_certificate y and s2,y2 = scale_certificate z in let g = gcd_big_int s1 s2 in let s1' = div_big_int s1 g in let s2' = div_big_int s2 g in - mult_big_int g (mult_big_int s1' s2'), + mult_big_int g (mult_big_int s1' s2'), Sum (Product(Rational_le (Big_int s2'), y1), Product (Rational_le (Big_int s1'), y2)) - | Product (y, z) -> + | Product (y, z) -> let s1,y1 = scale_certificate y and s2,y2 = scale_certificate z in mult_big_int s1 s2 , Product (y1,y2) @@ -743,7 +532,7 @@ open Micromega let rec term_to_q_expr = function | Const n -> PEc (Ml2C.q n) | Zero -> PEc ( Ml2C.q (Int 0)) - | Var s -> PEX (Ml2C.index + | Var s -> PEX (Ml2C.index (int_of_string (String.sub s 1 (String.length s - 1)))) | Mul(p1,p2) -> PEmul(term_to_q_expr p1, term_to_q_expr p2) | Add(p1,p2) -> PEadd(term_to_q_expr p1, term_to_q_expr p2) @@ -755,20 +544,20 @@ open Micromega let term_to_q_pol e = Mc.norm_aux (Ml2C.q (Int 0)) (Ml2C.q (Int 1)) Mc.qplus Mc.qmult Mc.qminus Mc.qopp Mc.qeq_bool (term_to_q_expr e) - let rec product l = + let rec product l = match l with | [] -> Mc.PsatzZ | [i] -> Mc.PsatzIn (Ml2C.nat i) | i ::l -> Mc.PsatzMulE(Mc.PsatzIn (Ml2C.nat i), product l) -let q_cert_of_pos pos = +let q_cert_of_pos pos = let rec _cert_of_pos = function Axiom_eq i -> Mc.PsatzIn (Ml2C.nat i) | Axiom_le i -> Mc.PsatzIn (Ml2C.nat i) | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i) | Monoid l -> product l - | Rational_eq n | Rational_le n | Rational_lt n -> + | Rational_eq n | Rational_le n | Rational_lt n -> if compare_num n (Int 0) = 0 then Mc.PsatzZ else Mc.PsatzC (Ml2C.q n) | Square t -> Mc.PsatzSquare (term_to_q_pol t) @@ -781,7 +570,7 @@ let q_cert_of_pos pos = let rec term_to_z_expr = function | Const n -> PEc (Ml2C.bigint (big_int_of_num n)) | Zero -> PEc ( Z0) - | Var s -> PEX (Ml2C.index + | Var s -> PEX (Ml2C.index (int_of_string (String.sub s 1 (String.length s - 1)))) | Mul(p1,p2) -> PEmul(term_to_z_expr p1, term_to_z_expr p2) | Add(p1,p2) -> PEadd(term_to_z_expr p1, term_to_z_expr p2) @@ -790,24 +579,649 @@ let q_cert_of_pos pos = | Sub(t1,t2) -> PEsub (term_to_z_expr t1, term_to_z_expr t2) | _ -> failwith "term_to_z_expr: not implemented" - let term_to_z_pol e = Mc.norm_aux (Ml2C.z 0) (Ml2C.z 1) Mc.zplus Mc.zmult Mc.zminus Mc.zopp Mc.zeq_bool (term_to_z_expr e) + let term_to_z_pol e = Mc.norm_aux (Ml2C.z 0) (Ml2C.z 1) Mc.Z.add Mc.Z.mul Mc.Z.sub Mc.Z.opp Mc.zeq_bool (term_to_z_expr e) -let z_cert_of_pos pos = +let z_cert_of_pos pos = let s,pos = (scale_certificate pos) in let rec _cert_of_pos = function Axiom_eq i -> Mc.PsatzIn (Ml2C.nat i) | Axiom_le i -> Mc.PsatzIn (Ml2C.nat i) | Axiom_lt i -> Mc.PsatzIn (Ml2C.nat i) | Monoid l -> product l - | Rational_eq n | Rational_le n | Rational_lt n -> + | Rational_eq n | Rational_le n | Rational_lt n -> if compare_num n (Int 0) = 0 then Mc.PsatzZ else Mc.PsatzC (Ml2C.bigint (big_int_of_num n)) | Square t -> Mc.PsatzSquare (term_to_z_pol t) - | Eqmul (t, y) -> Mc.PsatzMulC(term_to_z_pol t, _cert_of_pos y) + | Eqmul (t, y) -> + let is_unit = + match t with + | Const n -> n =/ Int 1 + | _ -> false in + if is_unit + then _cert_of_pos y + else Mc.PsatzMulC(term_to_z_pol t, _cert_of_pos y) | Sum (y, z) -> Mc.PsatzAdd (_cert_of_pos y, _cert_of_pos z) | Product (y, z) -> Mc.PsatzMulE (_cert_of_pos y, _cert_of_pos z) in simplify_cone z_spec (_cert_of_pos pos) +(** All constraints (initial or derived) have an index and have a justification i.e., proof. + Given a constraint, all the coefficients are always integers. +*) +open Mutils +open Mfourier +open Num +open Big_int +open Polynomial + +(*module Mc = Micromega*) +(*module Ml2C = Mutils.CamlToCoq +module C2Ml = Mutils.CoqToCaml +*) +let debug = false + + + +module Env = +struct + + type t = int list + + let id_of_hyp hyp l = + let rec xid_of_hyp i l = + match l with + | [] -> failwith "id_of_hyp" + | hyp'::l -> if hyp = hyp' then i else xid_of_hyp (i+1) l in + xid_of_hyp 0 l + +end + + +let coq_poly_of_linpol (p,c) = + + let pol_of_mon m = + Monomial.fold (fun x v p -> Mc.PEmul(Mc.PEpow(Mc.PEX(Ml2C.positive x),Ml2C.n v),p)) m (Mc.PEc (Mc.Zpos Mc.XH)) in + + List.fold_left (fun acc (x,v) -> + let mn = LinPoly.MonT.retrieve x in + Mc.PEadd(Mc.PEmul(Mc.PEc (Ml2C.bigint (numerator v)), pol_of_mon mn),acc)) (Mc.PEc (Ml2C.bigint (numerator c))) p + + + + +let rec cmpl_prf_rule env = function + | Hyp i | Def i -> Mc.PsatzIn (Ml2C.nat (Env.id_of_hyp i env)) + | Cst i -> Mc.PsatzC (Ml2C.bigint i) + | Zero -> Mc.PsatzZ + | MulPrf(p1,p2) -> Mc.PsatzMulE(cmpl_prf_rule env p1, cmpl_prf_rule env p2) + | AddPrf(p1,p2) -> Mc.PsatzAdd(cmpl_prf_rule env p1 , cmpl_prf_rule env p2) + | MulC(lp,p) -> let lp = Mc.norm0 (coq_poly_of_linpol lp) in + Mc.PsatzMulC(lp,cmpl_prf_rule env p) + | Square lp -> Mc.PsatzSquare (Mc.norm0 (coq_poly_of_linpol lp)) + | _ -> failwith "Cuts should already be compiled" + + +let rec cmpl_proof env = function + | Done -> Mc.DoneProof + | Step(i,p,prf) -> + begin + match p with + | CutPrf p' -> + Mc.CutProof(cmpl_prf_rule env p', cmpl_proof (i::env) prf) + | _ -> Mc.RatProof(cmpl_prf_rule env p,cmpl_proof (i::env) prf) + end + | Enum(i,p1,_,p2,l) -> + Mc.EnumProof(cmpl_prf_rule env p1,cmpl_prf_rule env p2,List.map (cmpl_proof (i::env)) l) + + +let compile_proof env prf = + let id = 1 + proof_max_id prf in + let _,prf = normalise_proof id prf in + if debug then Printf.fprintf stdout "compiled proof %a\n" output_proof prf; + cmpl_proof env prf + +type prf_sys = (cstr_compat * prf_rule) list + + +let xlinear_prover sys = + match Fourier.find_point sys with + | Inr prf -> + if debug then Printf.printf "AProof : %a\n" pp_proof prf ; + let cert = (*List.map (fun (x,n) -> x+1,n)*) (fst (List.hd (Proof.mk_proof sys prf))) in + if debug then Printf.printf "CProof : %a" Vect.pp_vect cert ; + Some (rats_to_ints (Vect.to_list cert)) + | Inl _ -> None + + +let output_num o n = output_string o (string_of_num n) +let output_bigint o n = output_string o (string_of_big_int n) + +let proof_of_farkas prf cert = +(* Printf.printf "\nproof_of_farkas %a , %a \n" (pp_list output_prf_rule) prf (pp_list output_bigint) cert ; *) + let rec mk_farkas acc prf cert = + match prf, cert with + | _ , [] -> acc + | [] , _ -> failwith "proof_of_farkas : not enough hyps" + | p::prf,c::cert -> + mk_farkas (add_proof (mul_proof c p) acc) prf cert in + let res = mk_farkas Zero prf cert in + (*Printf.printf "==> %a" output_prf_rule res ; *) + res + + +let linear_prover sys = + let (sysi,prfi) = List.split sys in + match xlinear_prover sysi with + | None -> None + | Some cert -> Some (proof_of_farkas prfi cert) + +let linear_prover = + if debug + then + fun sys -> + Printf.printf "<linear_prover"; flush stdout ; + let res = linear_prover sys in + Printf.printf ">"; flush stdout ; + res + else linear_prover + + + + +(** A single constraint can be unsat for the following reasons: + - 0 >= c for c a negative constant + - 0 = c for c a non-zero constant + - e = c when the coeffs of e are all integers and c is rational +*) + +type checksat = + | Tauto (* Tautology *) + | Unsat of prf_rule (* Unsatisfiable *) + | Cut of cstr_compat * prf_rule (* Cutting plane *) + | Normalise of cstr_compat * prf_rule (* coefficients are relatively prime *) + + +(** [check_sat] + - detects constraints that are not satisfiable; + - normalises constraints and generate cuts. +*) + +let check_sat (cstr,prf) = + let {coeffs=coeffs ; op=op ; cst=cst} = cstr in + match coeffs with + | [] -> + if eval_op op (Int 0) cst then Tauto else Unsat prf + | _ -> + let gcdi = (gcd_list (List.map snd coeffs)) in + let gcd = Big_int gcdi in + if eq_num gcd (Int 1) + then Normalise(cstr,prf) + else + if sign_num (mod_num cst gcd) = 0 + then (* We can really normalise *) + begin + assert (sign_num gcd >=1 ) ; + let cstr = { + coeffs = List.map (fun (x,v) -> (x, v // gcd)) coeffs; + op = op ; cst = cst // gcd + } in + Normalise(cstr,Gcd(gcdi,prf)) + (* Normalise(cstr,CutPrf prf)*) + end + else + match op with + | Eq -> Unsat (CutPrf prf) + | Ge -> + let cstr = { + coeffs = List.map (fun (x,v) -> (x, v // gcd)) coeffs; + op = op ; cst = ceiling_num (cst // gcd) + } in Cut(cstr,CutPrf prf) + + +(** Proof generating pivoting over variable v *) +let pivot v (c1,p1) (c2,p2) = + let {coeffs = v1 ; op = op1 ; cst = n1} = c1 + and {coeffs = v2 ; op = op2 ; cst = n2} = c2 in + + + + (* Could factorise gcd... *) + let xpivot cv1 cv2 = + ( + {coeffs = Vect.add (Vect.mul cv1 v1) (Vect.mul cv2 v2) ; + op = Proof.add_op op1 op2 ; + cst = n1 */ cv1 +/ n2 */ cv2 }, + + AddPrf(mul_proof (numerator cv1) p1,mul_proof (numerator cv2) p2)) in + + match Vect.get v v1 , Vect.get v v2 with + | None , _ | _ , None -> None + | Some a , Some b -> + if (sign_num a) * (sign_num b) = -1 + then + let cv1 = abs_num b + and cv2 = abs_num a in + Some (xpivot cv1 cv2) + else + if op1 = Eq + then + let cv1 = minus_num (b */ (Int (sign_num a))) + and cv2 = abs_num a in + Some (xpivot cv1 cv2) + else if op2 = Eq + then + let cv1 = abs_num b + and cv2 = minus_num (a */ (Int (sign_num b))) in + Some (xpivot cv1 cv2) + else None (* op2 could be Eq ... this might happen *) + +exception FoundProof of prf_rule + +let rec simpl_sys sys = + List.fold_left (fun acc (c,p) -> + match check_sat (c,p) with + | Tauto -> acc + | Unsat prf -> raise (FoundProof prf) + | Cut(c,p) -> (c,p)::acc + | Normalise (c,p) -> (c,p)::acc) [] sys + + +(** [ext_gcd a b] is the extended Euclid algorithm. + [ext_gcd a b = (x,y,g)] iff [ax+by=g] + Source: http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm +*) +let rec ext_gcd a b = + if sign_big_int b = 0 + then (unit_big_int,zero_big_int) + else + let (q,r) = quomod_big_int a b in + let (s,t) = ext_gcd b r in + (t, sub_big_int s (mult_big_int q t)) + + +let pp_ext_gcd a b = + let a' = big_int_of_int a in + let b' = big_int_of_int b in + + let (x,y) = ext_gcd a' b' in + Printf.fprintf stdout "%s * %s + %s * %s = %s\n" + (string_of_big_int x) (string_of_big_int a') + (string_of_big_int y) (string_of_big_int b') + (string_of_big_int (add_big_int (mult_big_int x a') (mult_big_int y b'))) + +exception Result of (int * (proof * cstr_compat)) + +let split_equations psys = + List.partition (fun (c,p) -> c.op = Eq) + + +let extract_coprime (c1,p1) (c2,p2) = + let rec exist2 vect1 vect2 = + match vect1 , vect2 with + | _ , [] | [], _ -> None + | (v1,n1)::vect1' , (v2, n2) :: vect2' -> + if v1 = v2 + then + if compare_big_int (gcd_big_int (numerator n1) (numerator n2)) unit_big_int = 0 + then Some (v1,n1,n2) + else + exist2 vect1' vect2' + else + if v1 < v2 + then exist2 vect1' vect2 + else exist2 vect1 vect2' in + + if c1.op = Eq && c2.op = Eq + then exist2 c1.coeffs c2.coeffs + else None + +let extract2 pred l = + let rec xextract2 rl l = + match l with + | [] -> (None,rl) (* Did not find *) + | e::l -> + match extract (pred e) l with + | None,_ -> xextract2 (e::rl) l + | Some (r,e'),l' -> Some (r,e,e'), List.rev_append rl l' in + + xextract2 [] l + + +let extract_coprime_equation psys = + extract2 extract_coprime psys + + +let apply_and_normalise f psys = + List.fold_left (fun acc pc' -> + match f pc' with + | None -> pc'::acc + | Some pc' -> + match check_sat pc' with + | Tauto -> acc + | Unsat prf -> raise (FoundProof prf) + | Cut(c,p) -> (c,p)::acc + | Normalise (c,p) -> (c,p)::acc + ) [] psys + + + + +let pivot_sys v pc psys = apply_and_normalise (pivot v pc) psys + + +let reduce_coprime psys = + let oeq,sys = extract_coprime_equation psys in + match oeq with + | None -> None (* Nothing to do *) + | Some((v,n1,n2),(c1,p1),(c2,p2) ) -> + let (l1,l2) = ext_gcd (numerator n1) (numerator n2) in + let l1' = Big_int l1 and l2' = Big_int l2 in + let cstr = + {coeffs = Vect.add (Vect.mul l1' c1.coeffs) (Vect.mul l2' c2.coeffs); + op = Eq ; + cst = (l1' */ c1.cst) +/ (l2' */ c2.cst) + } in + let prf = add_proof (mul_proof (numerator l1') p1) (mul_proof (numerator l2') p2) in + + Some (pivot_sys v (cstr,prf) ((c1,p1)::sys)) + +(** If there is an equation [eq] of the form 1.x + e = c, do a pivot over x with equation [eq] *) +let reduce_unary psys = + let is_unary_equation (cstr,prf) = + if cstr.op = Eq + then + try + Some (fst (List.find (fun (_,n) -> n =/ (Int 1) || n=/ (Int (-1))) cstr.coeffs)) + with Not_found -> None + else None in + + let (oeq,sys) = extract is_unary_equation psys in + match oeq with + | None -> None (* Nothing to do *) + | Some(v,pc) -> + Some(pivot_sys v pc sys) + +let reduce_non_lin_unary psys = + + let is_unary_equation (cstr,prf) = + if cstr.op = Eq + then + try + let x = fst (List.find (fun (x,n) -> (n =/ (Int 1) || n=/ (Int (-1))) && Monomial.is_var (LinPoly.MonT.retrieve x) ) cstr.coeffs) in + let x' = LinPoly.MonT.retrieve x in + if List.for_all (fun (y,_) -> y = x || snd (Monomial.div (LinPoly.MonT.retrieve y) x') = 0) cstr.coeffs + then Some x + else None + with Not_found -> None + else None in + + + let (oeq,sys) = extract is_unary_equation psys in + match oeq with + | None -> None (* Nothing to do *) + | Some(v,pc) -> + Some(apply_and_normalise (LinPoly.pivot_eq v pc) sys) + +let reduce_var_change psys = + + let rec rel_prime vect = + match vect with + | [] -> None + | (x,v)::vect -> + let v = numerator v in + try + let (x',v') = List.find (fun (_,v') -> + let v' = numerator v' in + eq_big_int (gcd_big_int v v') unit_big_int) vect in + Some ((x,v),(x',numerator v')) + with Not_found -> rel_prime vect in + + let rel_prime (cstr,prf) = if cstr.op = Eq then rel_prime cstr.coeffs else None in + + let (oeq,sys) = extract rel_prime psys in + + match oeq with + | None -> None + | Some(((x,v),(x',v')),(c,p)) -> + let (l1,l2) = ext_gcd v v' in + let l1,l2 = Big_int l1 , Big_int l2 in + + let get v vect = + match Vect.get v vect with + | None -> Int 0 + | Some n -> n in + + let pivot_eq (c',p') = + let {coeffs = coeffs ; op = op ; cst = cst} = c' in + let vx = get x coeffs in + let vx' = get x' coeffs in + let m = minus_num (vx */ l1 +/ vx' */ l2) in + Some ({coeffs = + Vect.add (Vect.mul m c.coeffs) coeffs ; op = op ; cst = m */ c.cst +/ cst} , + AddPrf(MulC(([], m),p),p')) in + + Some (apply_and_normalise pivot_eq sys) + + + + + let reduce_pivot psys = + let is_equation (cstr,prf) = + if cstr.op = Eq + then + try + Some (fst (List.hd cstr.coeffs)) + with Not_found -> None + else None in + let (oeq,sys) = extract is_equation psys in + match oeq with + | None -> None (* Nothing to do *) + | Some(v,pc) -> + if debug then + Printf.printf "Bad news : loss of completeness %a=%s" Vect.pp_vect (fst pc).coeffs (string_of_num (fst pc).cst); + Some(pivot_sys v pc sys) + + + + + + let iterate_until_stable f x = + let rec iter x = + match f x with + | None -> x + | Some x' -> iter x' in + iter x + + let rec app_funs l x = + match l with + | [] -> None + | f::fl -> + match f x with + | None -> app_funs fl x + | Some x' -> Some x' + + let reduction_equations psys = + iterate_until_stable (app_funs + [reduce_unary ; reduce_coprime ; + reduce_var_change (*; reduce_pivot*)]) psys + + let reduction_non_lin_equations psys = + iterate_until_stable (app_funs + [reduce_non_lin_unary (*; reduce_coprime ; + reduce_var_change ; reduce_pivot *)]) psys + + + + + (** [get_bound sys] returns upon success an interval (lb,e,ub) with proofs *) + let get_bound sys = + let is_small (v,i) = + match Itv.range i with + | None -> false + | Some i -> i <=/ (Int 1) in + + let select_best (x1,i1) (x2,i2) = + if Itv.smaller_itv i1 i2 + then (x1,i1) else (x2,i2) in + + (* For lia, there are no equations => these precautions are not needed *) + (* For nlia, there are equations => do not enumerate over equations! *) + let all_planes sys = + let (eq,ineq) = List.partition (fun c -> c.op = Eq) sys in + match eq with + | [] -> List.rev_map (fun c -> c.coeffs) ineq + | _ -> + List.fold_left (fun acc c -> + if List.exists (fun c' -> Vect.equal c.coeffs c'.coeffs) eq + then acc else c.coeffs ::acc) [] ineq in + + let smallest_interval = + List.fold_left + (fun acc vect -> + if is_small acc + then acc + else + match Fourier.optimise vect sys with + | None -> acc + | Some i -> + if debug then Printf.printf "Found a new bound %a" Vect.pp_vect vect ; + select_best (vect,i) acc) (Vect.null, (None,None)) (all_planes sys) in + let smallest_interval = + match smallest_interval + with + | (x,(Some i, Some j)) -> Some(i,x,j) + | x -> None (* This should not be possible *) + in + match smallest_interval with + | Some (lb,e,ub) -> + let (lbn,lbd) = (sub_big_int (numerator lb) unit_big_int, denominator lb) in + let (ubn,ubd) = (add_big_int unit_big_int (numerator ub) , denominator ub) in + (match + (* x <= ub -> x > ub *) + xlinear_prover ({coeffs = Vect.mul (Big_int ubd) e ; op = Ge ; cst = Big_int ubn} :: sys), + (* lb <= x -> lb > x *) + xlinear_prover + ({coeffs = Vect.mul (minus_num (Big_int lbd)) e ; op = Ge ; cst = minus_num (Big_int lbn)} :: sys) + with + | Some cub , Some clb -> Some(List.tl clb,(lb,e,ub), List.tl cub) + | _ -> failwith "Interval without proof" + ) + | None -> None + + + let check_sys sys = + List.for_all (fun (c,p) -> List.for_all (fun (_,n) -> sign_num n <> 0) c.coeffs) sys + + + let xlia reduction_equations sys = + + let rec enum_proof (id:int) (sys:prf_sys) : proof option = + if debug then (Printf.printf "enum_proof\n" ; flush stdout) ; + assert (check_sys sys) ; + + let nsys,prf = List.split sys in + match get_bound nsys with + | None -> None (* Is the systeme really unbounded ? *) + | Some(prf1,(lb,e,ub),prf2) -> + if debug then Printf.printf "Found interval: %a in [%s;%s] -> " Vect.pp_vect e (string_of_num lb) (string_of_num ub) ; + (match start_enum id e (ceiling_num lb) (floor_num ub) sys + with + | Some prfl -> + Some(Enum(id,proof_of_farkas prf prf1,e, proof_of_farkas prf prf2,prfl)) + | None -> None + ) + + and start_enum id e clb cub sys = + if clb >/ cub + then Some [] + else + let eq = {coeffs = e ; op = Eq ; cst = clb} in + match aux_lia (id+1) ((eq, Def id) :: sys) with + | None -> None + | Some prf -> + match start_enum id e (clb +/ (Int 1)) cub sys with + | None -> None + | Some l -> Some (prf::l) + + and aux_lia (id:int) (sys:prf_sys) : proof option = + assert (check_sys sys) ; + if debug then Printf.printf "xlia: %a \n" (pp_list (fun o (c,_) -> output_cstr o c)) sys ; + try + let sys = reduction_equations sys in + if debug then + Printf.printf "after reduction: %a \n" (pp_list (fun o (c,_) -> output_cstr o c)) sys ; + match linear_prover sys with + | Some prf -> Some (Step(id,prf,Done)) + | None -> enum_proof id sys + with FoundProof prf -> + (* [reduction_equations] can find a proof *) + Some(Step(id,prf,Done)) in + + (* let sys' = List.map (fun (p,o) -> Mc.norm0 p , o) sys in*) + let id = List.length sys in + let orpf = + try + let sys = simpl_sys sys in + aux_lia id sys + with FoundProof pr -> Some(Step(id,pr,Done)) in + match orpf with + | None -> None + | Some prf -> + (*Printf.printf "direct proof %a\n" output_proof prf ; *) + let env = mapi (fun _ i -> i) sys in + let prf = compile_proof env prf in + (*try + if Mc.zChecker sys' prf then Some prf else + raise Certificate.BadCertificate + with Failure s -> (Printf.printf "%s" s ; Some prf) + *) Some prf + + + let cstr_compat_of_poly (p,o) = + let (v,c) = LinPoly.linpol_of_pol p in + {coeffs = v ; op = o ; cst = minus_num c } + + + let lia sys = + LinPoly.MonT.clear (); + let sys = List.map (develop_constraint z_spec) sys in + let (sys:cstr_compat list) = List.map cstr_compat_of_poly sys in + let sys = mapi (fun c i -> (c,Hyp i)) sys in + xlia reduction_equations sys + + + let nlia sys = + LinPoly.MonT.clear (); + let sys = List.map (develop_constraint z_spec) sys in + let sys = mapi (fun c i -> (c,Hyp i)) sys in + + let is_linear = List.for_all (fun ((p,_),_) -> Poly.is_linear p) sys in + + let module MonMap = Map.Make(Monomial) in + + let collect_square = + List.fold_left (fun acc ((p,_),_) -> Poly.fold + (fun m _ acc -> + match Monomial.sqrt m with + | None -> acc + | Some s -> MonMap.add s m acc) p acc) MonMap.empty sys in + let sys = MonMap.fold (fun s m acc -> + let s = LinPoly.linpol_of_pol (Poly.add s (Int 1) (Poly.constant (Int 0))) in + let m = Poly.add m (Int 1) (Poly.constant (Int 0)) in + ((m, Ge), (Square s))::acc) collect_square sys in + +(* List.iter (fun ((p,_),_) -> Printf.printf "square %a\n" Poly.pp p) gen_square*) + + let sys = + if is_linear then sys + else sys @ (all_sym_pairs (fun ((c,o),p) ((c',o'),p') -> + ((Poly.product c c',opMult o o'), MulPrf(p,p'))) sys) in + + let sys = List.map (fun (c,p) -> cstr_compat_of_poly c,p) sys in + assert (check_sys sys) ; + xlia (if is_linear then reduction_equations else reduction_non_lin_equations) sys + + + (* Local Variables: *) (* coding: utf-8 *) (* End: *) |