diff options
Diffstat (limited to 'contrib/micromega/certificate.ml')
-rw-r--r-- | contrib/micromega/certificate.ml | 618 |
1 files changed, 618 insertions, 0 deletions
diff --git a/contrib/micromega/certificate.ml b/contrib/micromega/certificate.ml new file mode 100644 index 00000000..88e882e6 --- /dev/null +++ b/contrib/micromega/certificate.ml @@ -0,0 +1,618 @@ +(************************************************************************) +(* 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 *) +(* *) +(************************************************************************) + +(* We take as input a list of polynomials [p1...pn] and return an unfeasibility + certificate polynomial. *) + +(*open Micromega.Polynomial*) +open Big_int +open Num + +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 +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 compare = compare compare_num +end + +open Mutils +type 'a number_spec = { + bigint_to_number : big_int -> 'a; + number_to_num : 'a -> num; + zero : 'a; + unit : 'a; + mult : 'a -> 'a -> 'a; + eqb : 'a -> 'a -> Mc.bool +} + +let z_spec = { + bigint_to_number = Ml2C.bigint ; + number_to_num = (fun x -> Big_int (C2Ml.z_big_int x)); + zero = Mc.Z0; + unit = Mc.Zpos Mc.XH; + mult = Mc.zmult; + eqb = Mc.zeq_bool +} + + +let q_spec = { + bigint_to_number = (fun x -> {Mc.qnum = Ml2C.bigint x; Mc.qden = Mc.XH}); + number_to_num = C2Ml.q_to_num; + zero = {Mc.qnum = Mc.Z0;Mc.qden = Mc.XH}; + unit = {Mc.qnum = (Mc.Zpos Mc.XH) ; Mc.qden = Mc.XH}; + mult = Mc.qmult; + eqb = Mc.qeq_bool +} + +let r_spec = z_spec + + + + +let dev_form n_spec 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) -> + let p1 = dev_form p1 in + let p2 = dev_form p2 in + 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) -> + let p = dev_form p in + let n = C2Ml.n n in + 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 list_to_polynomial vars l = + assert (List.for_all (fun x -> ceiling_num x =/ x) l); + let var x = monomial_to_polynomial (List.nth vars x) in + let rec xtopoly p i = function + | [] -> p + | c::l -> if c =/ (Int 0) then xtopoly p (i+1) l + else let c = Mc.PEc (Ml2C.bigint (numerator c)) in + let mn = + if 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 = + let y' = f x in + if y' = x then y' + else fixpoint f y' + + + + + + + + +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.S_Mult(t1, t2) -> + simpl_cone (Mc.S_Mult (rec_simpl_cone t1, rec_simpl_cone t2)) + | Mc.S_Add(t1,t2) -> + simpl_cone (Mc.S_Add (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 + | Other of cone +and cone = Mc.zWitness + + + +let factorise_linear_cone c = + + let rec cone_list c l = + match c with + | Mc.S_Add (x,r) -> cone_list r (x::l) + | _ -> c :: l in + + let factorise c1 c2 = + match c1 , c2 with + | Mc.S_Ideal(x,y) , Mc.S_Ideal(x',y') -> + if x = x' then Some (Mc.S_Ideal(x, Mc.S_Add(y,y'))) else None + | Mc.S_Mult(x,y) , Mc.S_Mult(x',y') -> + if x = x' then Some (Mc.S_Mult(x, Mc.S_Add(y,y'))) else None + | _ -> None in + + let rec rebuild_cone l pending = + match l with + | [] -> (match pending with + | None -> Mc.S_Z + | Some p -> p + ) + | e::l -> + (match pending with + | None -> rebuild_cone l (Some e) + | Some p -> (match factorise p e with + | None -> Mc.S_Add(p, rebuild_cone l (Some e)) + | Some f -> rebuild_cone l (Some f) ) + ) in + + (rebuild_cone (List.sort Pervasives.compare (cone_list c [])) None) + + + +(* 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 + where pi >= 0 qi > 0 + ai >= 0 + bi >= 0 + Sum bi + c >= 1 + This is a linear problem: each monomial is considered as a variable. + Hence, we can use fourier. + + The variable c is at index 0 +*) + +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 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 ; + cst = Big_int zero_big_int } + else + { 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 = + match l with + | [] -> [] + | (_,Mc.Equal)::l -> xpositivity (i+1) l + | (_,_)::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.Equal -> "= 0" + | Mc.NonEqual -> "<> 0" + + + +(* 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 *) + let l' = List.map fst l in + let monomials = + List.fold_left (fun acc p -> Poly.addition p acc) (Poly.constant (Int 0)) l' + in (* For each monomial, compute a constraint *) + let s0 = + Poly.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 + | _ -> 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 ; + cst = Big_int zero_big_int}::(strict::(positivity l)@s0) + + +let big_int_to_z = Ml2C.bigint + +(* For Q, this is a pity that the certificate has been scaled + -- at a lower layer, certificates are using nums... *) +let make_certificate n_spec cert li = + let bint_to_cst = n_spec.bigint_to_number in + match cert with + | [] -> None + | e::cert' -> + let cst = match compare_big_int e zero_big_int with + | 0 -> Mc.S_Z + | 1 -> Mc.S_Pos (bint_to_cst e) + | _ -> failwith "positivity error" + in + let rec scalar_product cert l = + match cert with + | [] -> Mc.S_Z + | c::cert -> match l with + | [] -> failwith "make_certificate(1)" + | i::l -> + let r = scalar_product cert l in + match compare_big_int c zero_big_int with + | -1 -> Mc.S_Add ( + Mc.S_Ideal (Mc.PEc ( bint_to_cst c), Mc.S_In (Ml2C.nat i)), + r) + | 0 -> r + | _ -> Mc.S_Add ( + Mc.S_Mult (Mc.S_Pos (bint_to_cst c), Mc.S_In (Ml2C.nat i)), + r) in + + Some ((factorise_linear_cone + (simplify_cone n_spec (Mc.S_Add (cst, scalar_product cert' li))))) + + +exception Found of Monomial.t + +let raw_certificate l = + let sys = build_linear_system l in + try + match Fourier.find_point sys with + | None -> None + | Some 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); + flush stdout) ; + None + + +let simple_linear_prover to_constant l = + let (lc,li) = List.split l in + match raw_certificate lc with + | None -> None (* No certificate *) + | Some cert -> make_certificate to_constant 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 (c,i) -> let (Mc.Pair(x,y)) = c in + match y with + Mc.NonEqual -> failwith "cannot happen" + | y -> ((dev_form n_spec x, y),i)) l' in + + 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) + +(* 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) + (Poly.constant (Int 0)) l' in + 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 ; + cst = minus_num ( (Poly.get Monomial.const c))}) 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 +let popp x = Mc.PEopp x + +let debug = false + +(* keep track of enumerated vectors *) +let rec mem p x l = + match l with [] -> false | e::l -> if p x e then true else mem p x l + +let rec remove_assoc p x l = + match l with [] -> [] | e::l -> if p x (fst e) then + remove_assoc p x l else e::(remove_assoc p x l) + +let eq x y = Vect.compare x y = 0 + +(* Beurk... this code is a shame *) + +let rec zlinear_prover sys = xzlinear_prover [] sys + +and xzlinear_prover enum l : (Mc.proofTerm option) = + match linear_prover z_spec l with + | Some prf -> Some (Mc.RatProof prf) + | None -> + let ll = List.fold_right (fun (Mc.Pair(e,k)) r -> match k with + Mc.NonEqual -> r + | k -> (dev_form z_spec e , + match k with + | Mc.Strict | Mc.NonStrict -> Ge + (* Loss of precision -- weakness of fourier*) + | Mc.Equal -> Eq + | Mc.NonEqual -> failwith "Cannot happen") :: r) l [] in + + let (sys,var) = make_linear_system ll in + let res = + match Fourier.find_Q_interval sys with + | Some(i,x,j) -> if i =/ j + then Some(i,Vect.set x (Int 1) Vect.null,i) else None + | None -> None in + let res = match res with + | None -> + begin + let candidates = List.fold_right + (fun cstr acc -> + let gcd = Big_int (Vect.gcd cstr.coeffs) in + let vect = Vect.mul (Int 1 // gcd) cstr.coeffs in + if mem eq vect enum then acc + else ((vect,Fourier.optimise vect sys)::acc)) sys [] in + let candidates = List.fold_left (fun l (x,i) -> + match i with + None -> (x,Empty)::l + | Some i -> (x,i)::l) [] (candidates) in + match List.fold_left (fun (x1,i1) (x2,i2) -> + if smaller_itv i1 i2 + then (x1,i1) else (x2,i2)) (Vect.null,Itv(None,None)) candidates + with + | (i,Empty) -> None + | (x,Itv(Some i, Some j)) -> Some(i,x,j) + | (x,Point n) -> Some(n,x,n) + | x -> match Fourier.find_Q_interval sys with + | None -> None + | Some(i,x,j) -> + if i =/ j + then Some(i,Vect.set x (Int 1) Vect.null,i) + else None + end + | _ -> res in + + match res with + | Some (lb,e,ub) -> + let (lbn,lbd) = + (Ml2C.bigint (sub_big_int (numerator lb) unit_big_int), + Ml2C.bigint (denominator lb)) in + let (ubn,ubd) = + (Ml2C.bigint (add_big_int unit_big_int (numerator ub)) , + Ml2C.bigint (denominator ub)) in + let expr = list_to_polynomial var (Vect.to_list e) in + (match + (*x <= ub -> x > ub *) + linear_prover z_spec + (Mc.Pair(pplus (pmult (pconst ubd) expr) (popp (pconst ubn)), + Mc.NonStrict) :: l), + (* lb <= x -> lb > x *) + linear_prover z_spec + (Mc.Pair( pplus (popp (pmult (pconst lbd) expr)) (pconst lbn) , + Mc.NonStrict)::l) + with + | Some cub , Some clb -> + (match zlinear_enum (e::enum) expr + (ceiling_num lb) (floor_num ub) l + with + | None -> None + | Some prf -> + Some (Mc.EnumProof(Ml2C.q lb,expr,Ml2C.q ub,clb,cub,prf))) + | _ -> None + ) + | _ -> None +and xzlinear_enum enum expr clb cub l = + if clb >/ cub + then Some Mc.Nil + else + let pexpr = pplus (popp (pconst (Ml2C.bigint (numerator clb)))) expr in + let sys' = (Mc.Pair(pexpr, Mc.Equal))::l in + match xzlinear_prover enum sys' with + | None -> if debug then print_string "zlp?"; None + | Some prf -> if debug then print_string "zlp!"; + match zlinear_enum enum expr (clb +/ (Int 1)) cub l with + | None -> None + | Some prfl -> Some (Mc.Cons(prf,prfl)) + +and zlinear_enum enum expr clb cub l = + let res = xzlinear_enum enum expr clb cub l in + if debug then Printf.printf "zlinear_enum %s %s -> %s\n" + (string_of_num clb) + (string_of_num cub) + (match res with + | None -> "None" + | Some r -> "Some") ; res + |