diff options
Diffstat (limited to 'plugins/micromega/polynomial.ml')
| -rw-r--r-- | plugins/micromega/polynomial.ml | 1336 |
1 files changed, 785 insertions, 551 deletions
diff --git a/plugins/micromega/polynomial.ml b/plugins/micromega/polynomial.ml index 1d18a26f33..5f31b6f145 100644 --- a/plugins/micromega/polynomial.ml +++ b/plugins/micromega/polynomial.ml @@ -10,7 +10,7 @@ (* *) (* Micromega: A reflexive tactic using the Positivstellensatz *) (* *) -(* Frédéric Besson (Irisa/Inria) 2006-20011 *) +(* Frédéric Besson (Irisa/Inria) 2006-20018 *) (* *) (************************************************************************) @@ -18,6 +18,10 @@ open Num module Utils = Mutils open Utils +module Mc = Micromega + +let max_nb_cstr = ref max_int + type var = int let debug = false @@ -25,652 +29,882 @@ let debug = false let (<+>) = add_num let (<*>) = mult_num - module Monomial : sig - type t - val const : t - val is_const : t -> bool - val var : var -> t - val is_var : t -> bool - val prod : t -> t -> t - val exp : t -> int -> t - val div : t -> t -> t * int - val compare : t -> t -> int - val pp : out_channel -> t -> unit - val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a - val sqrt : t -> t option + type t + val const : t + val is_const : t -> bool + val var : var -> t + val is_var : t -> bool + val get_var : t -> var option + val prod : t -> t -> t + val exp : t -> int -> t + val div : t -> t -> t * int + val compare : t -> t -> int + val pp : out_channel -> t -> unit + val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a + val sqrt : t -> t option + val variables : t -> ISet.t end - = -struct - (* A monomial is represented by a multiset of variables *) - module Map = Map.Make(Int) - open Map - - type t = int Map.t - - let pp o m = Map.iter - (fun k v -> - if v = 1 then Printf.fprintf o "x%i." k - else Printf.fprintf o "x%i^%i." k v) m - - - (* The monomial that corresponds to a constant *) - let const = Map.empty - - let sum_degree m = Map.fold (fun _ n s -> s + n) m 0 - - (* Total ordering of monomials *) - let compare: t -> t -> int = - fun m1 m2 -> - let s1 = sum_degree m1 - and s2 = sum_degree m2 in - if Int.equal s1 s2 then Map.compare Int.compare m1 m2 - else Int.compare s1 s2 - - let is_const m = (m = Map.empty) - - (* The monomial 'x' *) - let var x = Map.add x 1 Map.empty - - let is_var m = - try - not (Map.fold (fun _ i fk -> - if fk = true (* first key *) - then - if i = 1 then false - else raise Not_found - else raise Not_found) m true) - with Not_found -> false - - let sqrt m = - if is_const m then None - else - try - Some (Map.fold (fun v i acc -> - let i' = i / 2 in - if i mod 2 = 0 - then add v i' m - else raise Not_found) m const) - with Not_found -> None - - (* Get the degre of a variable in a monomial *) - let find x m = try find x m with Not_found -> 0 - - (* Product of monomials *) - let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2 - - - let exp m n = - let rec exp acc n = - if n = 0 then acc - else exp (prod acc m) (n - 1) in - - exp const n - - - (* [div m1 m2 = mr,n] such that mr * (m2)^n = m1 *) - let div m1 m2 = - let n = fold (fun x i n -> let i' = find x m1 in - let nx = i' / i in - min n nx) m2 max_int in - - let mr = fold (fun x i' m -> - let i = find x m2 in - let ir = i' - i * n in - if ir = 0 then m - else add x ir m) m1 empty in - (mr,n) - - - let fold = fold + = struct + (* A monomial is represented by a multiset of variables *) + module Map = Map.Make(Int) + open Map + + type t = int Map.t + + let is_singleton m = + try + let (k,v) = choose m in + let (l,e,r) = split k m in + if is_empty l && is_empty r + then Some(k,v) else None + with Not_found -> None + + let pp o m = + let pp_elt o (k,v)= + if v = 1 then Printf.fprintf o "x%i" k + else Printf.fprintf o "x%i^%i" k v in + + let rec pp_list o l = + match l with + [] -> () + | [e] -> pp_elt o e + | e::l -> Printf.fprintf o "%a*%a" pp_elt e pp_list l in + + pp_list o (Map.bindings m) + + + + (* The monomial that corresponds to a constant *) + let const = Map.empty + + let sum_degree m = Map.fold (fun _ n s -> s + n) m 0 + + (* Total ordering of monomials *) + let compare: t -> t -> int = + fun m1 m2 -> + let s1 = sum_degree m1 + and s2 = sum_degree m2 in + if Int.equal s1 s2 then Map.compare Int.compare m1 m2 + else Int.compare s1 s2 + + let is_const m = (m = Map.empty) + + (* The monomial 'x' *) + let var x = Map.add x 1 Map.empty + + let is_var m = + match is_singleton m with + | None -> false + | Some (_,i) -> i = 1 + + let get_var m = + match is_singleton m with + | None -> None + | Some (k,i) -> if i = 1 then Some k else None + + + let sqrt m = + if is_const m then None + else + try + Some (Map.fold (fun v i acc -> + let i' = i / 2 in + if i mod 2 = 0 + then add v i' acc + else raise Not_found) m const) + with Not_found -> None + + + (* Get the degre of a variable in a monomial *) + let find x m = try find x m with Not_found -> 0 + + (* Product of monomials *) + let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2 + + let exp m n = + let rec exp acc n = + if n = 0 then acc + else exp (prod acc m) (n - 1) in + + exp const n + + (* [div m1 m2 = mr,n] such that mr * (m2)^n = m1 *) + let div m1 m2 = + let n = fold (fun x i n -> let i' = find x m1 in + let nx = i' / i in + min n nx) m2 max_int in + + let mr = fold (fun x i' m -> + let i = find x m2 in + let ir = i' - i * n in + if ir = 0 then m + else add x ir m) m1 empty in + (mr,n) + + + let variables m = fold (fun v i acc -> ISet.add v acc) m ISet.empty + + let fold = fold end +module MonMap = + struct + include Map.Make(Monomial) + + let union f = merge + (fun x v1 v2 -> + match v1 , v2 with + | None , None -> None + | Some v , None | None , Some v -> Some v + | Some v1 , Some v2 -> f x v1 v2) + end + +let pp_mon o (m, i) = + if Monomial.is_const m + then if eq_num (Int 0) i then () + else Printf.fprintf o "%s" (string_of_num i) + else + match i with + | Int 1 -> Monomial.pp o m + | Int -1 -> Printf.fprintf o "-%a" Monomial.pp m + | Int 0 -> () + | _ -> Printf.fprintf o "%s*%a" (string_of_num i) Monomial.pp m + + + module Poly : - (* A polynomial is a map of monomials *) - (* - This is probably a naive implementation +(* 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 product : t -> t -> t - val addition : t -> t -> t - val uminus : t -> t - val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a - val is_linear : t -> bool -end = -struct - (*normalisation bug : 0*x ... *) - module P = Map.Make(Monomial) - open P - - type t = num P.t - - (* 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 -> + type t + val pp : out_channel -> t -> unit + val get : Monomial.t -> t -> num + val variable : var -> t + val add : Monomial.t -> num -> t -> t + val constant : num -> 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 is_linear : t -> bool + val variables : t -> ISet.t + val factorise : var -> t -> t * t +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 mn i -> Printf.fprintf o "%a + " pp_mon (mn, i)) 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 -> if sign_num v = 0 then p else let vl = (get mn p) <+> v in - if sign_num vl = 0 then - remove mn p - else add mn vl p + if sign_num vl = 0 then + remove mn p + else add mn vl p - (** Design choice: empty is not a polynomial - I do not remember why .... - **) + (** 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 -> - if sign_num v = 0 + (* The product by a monomial *) + let mult : Monomial.t -> num -> t -> t = + fun mn v p -> + if sign_num v = 0 then constant (Int 0) else 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 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 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 uminus : t -> t = + fun p -> map (fun v -> minus_num v) p - let is_linear p = P.fold (fun m _ acc -> acc && (Monomial.is_const m || Monomial.is_var m)) p true + let fold = P.fold -(* let is_linear p = - let res = is_linear p in - Printf.printf "is_linear %a = %b\n" pp p res ; res -*) -end + let is_linear p = P.fold (fun m _ acc -> acc && (Monomial.is_const m || Monomial.is_var m)) p true + let variables p = P.fold (fun m _ acc -> ISet.union (Monomial.variables m) acc) p ISet.empty + + let factorise x p = + let x = Monomial.var x in + P.fold (fun m v (px,cx) -> + let (m1,i) = Monomial.div m x in + if i = 0 + then (px, add m v cx) + else + let mx = Monomial.prod m1 (Monomial.exp x (i-1)) in + (add mx v px,cx) ) p (constant (Int 0) , constant (Int 0)) + +end -module Vect = - struct - (** [t] is the type of vectors. - A vector [(x1,v1) ; ... ; (xn,vn)] is such that: - - variables indexes are ordered (x1 <c ... < xn - - values are all non-zero - *) - type var = int - type t = (var * num) list - -(** [equal v1 v2 = true] if the vectors are syntactically equal. *) - - let rec equal v1 v2 = - match v1 , v2 with - | [] , [] -> true - | [] , _ -> false - | _::_ , [] -> false - | (i1,n1)::v1 , (i2,n2)::v2 -> - (Int.equal 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 ) - - - 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 - - xfrom_list 0 l - - let zero_num = Int 0 - - - 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 Int.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 Int.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 mul z t = - match z with - | Int 0 -> [] - | Int 1 -> t - | _ -> List.map (fun (i,n) -> (i, mult_num z n)) t - - - let rec add v1 v2 = - match v1 , v2 with - | (x1,n1)::v1' , (x2,n2)::v2' -> - if x1 = x2 - then - let n' = n1 +/ n2 in - if n' =/ Int 0 then add v1' v2' - else - let res = add v1' v2' in - (x1,n') ::res - else if x1 < x2 - then let res = add v1' v2 in - (x1, n1)::res - else let res = add v1 v2' in - (x2, n2)::res - | [] , [] -> [] - | [] , _ -> v2 - | _ , [] -> v1 - - - - - let compare : t -> t -> int = Mutils.Cmp.compare_list (fun x y -> Mutils.Cmp.compare_lexical - [ - (fun () -> Int.compare (fst x) (fst y)); - (fun () -> compare_num (snd x) (snd y))]) - - (** [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 Int.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 type vector = Vect.t -type cstr_compat = {coeffs : vector ; op : op ; cst : num} -and op = |Eq | Ge +type cstr = {coeffs : vector ; op : op ; cst : num} +and op = |Eq | Ge | Gt -let string_of_op = function Eq -> "=" | Ge -> ">=" +exception Strict -let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} = - Printf.fprintf o "%a %s %s" Vect.pp_vect coeffs (string_of_op op) (string_of_num cst) +let is_strict c = Pervasives.(=) c.op Gt -let opMult o1 o2 = - match o1, o2 with - | Eq , Eq -> Eq - | Eq , Ge | Ge , Eq -> Ge - | Ge , Ge -> Ge - -open Big_int - -type prf_rule = - | Hyp of int - | Def of int - | Cst of big_int - | Zero - | Square of (Vect.t * num) - | MulC of (Vect.t * num) * prf_rule - | Gcd of big_int * prf_rule - | MulPrf of prf_rule * prf_rule - | AddPrf of prf_rule * prf_rule - | CutPrf of prf_rule - -type proof = - | Done - | Step of int * prf_rule * proof - | Enum of int * prf_rule * Vect.t * prf_rule * proof list - - -let rec output_prf_rule o = function - | Hyp i -> Printf.fprintf o "Hyp %i" i - | Def i -> Printf.fprintf o "Def %i" i - | Cst c -> Printf.fprintf o "Cst %s" (string_of_big_int c) - | Zero -> Printf.fprintf o "Zero" - | Square _ -> Printf.fprintf o "( )^2" - | MulC(p,pr) -> Printf.fprintf o "P * %a" output_prf_rule pr - | MulPrf(p1,p2) -> Printf.fprintf o "%a * %a" output_prf_rule p1 output_prf_rule p2 - | AddPrf(p1,p2) -> Printf.fprintf o "%a + %a" output_prf_rule p1 output_prf_rule p2 - | CutPrf(p) -> Printf.fprintf o "[%a]" output_prf_rule p - | Gcd(c,p) -> Printf.fprintf o "(%a)/%s" output_prf_rule p (string_of_big_int c) - -let rec output_proof o = function - | Done -> Printf.fprintf o "." - | Step(i,p,pf) -> Printf.fprintf o "%i:= %a ; %a" i output_prf_rule p output_proof pf - | Enum(i,p1,v,p2,pl) -> Printf.fprintf o "%i{%a<=%a<=%a}%a" i - output_prf_rule p1 Vect.pp_vect v output_prf_rule p2 - (pp_list output_proof) pl - -let rec pr_rule_max_id = function - | Hyp i | Def i -> i - | Cst _ | Zero | Square _ -> -1 - | MulC(_,p) | CutPrf p | Gcd(_,p) -> pr_rule_max_id p - | MulPrf(p1,p2)| AddPrf(p1,p2) -> max (pr_rule_max_id p1) (pr_rule_max_id p2) - -let rec proof_max_id = function - | Done -> -1 - | Step(i,pr,prf) -> max i (max (pr_rule_max_id pr) (proof_max_id prf)) - | Enum(i,p1,_,p2,l) -> - let m = max (pr_rule_max_id p1) (pr_rule_max_id p2) in - List.fold_left (fun i prf -> max i (proof_max_id prf)) (max i m) l - -let rec pr_rule_def_cut id = function - | MulC(p,prf) -> - let (bds,id',prf') = pr_rule_def_cut id prf in - (bds, id', MulC(p,prf')) - | MulPrf(p1,p2) -> - let (bds1,id,p1) = pr_rule_def_cut id p1 in - let (bds2,id,p2) = pr_rule_def_cut id p2 in - (bds2@bds1,id,MulPrf(p1,p2)) - | AddPrf(p1,p2) -> - let (bds1,id,p1) = pr_rule_def_cut id p1 in - let (bds2,id,p2) = pr_rule_def_cut id p2 in - (bds2@bds1,id,AddPrf(p1,p2)) - | CutPrf p -> - let (bds,id,p) = pr_rule_def_cut id p in - ((id,p)::bds,id+1,Def id) - | Gcd(c,p) -> - let (bds,id,p) = pr_rule_def_cut id p in - ((id,p)::bds,id+1,Def id) - | Square _|Cst _|Def _|Hyp _|Zero as x -> ([],id,x) - - -(* Do not define top-level cuts *) -let pr_rule_def_cut id = function - | CutPrf p -> - let (bds,ids,p') = pr_rule_def_cut id p in - bds,ids, CutPrf p' - | p -> pr_rule_def_cut id p - - -let rec implicit_cut p = - match p with - | CutPrf p -> implicit_cut p - | _ -> p - - -let rec normalise_proof id prf = - match prf with - | Done -> (id,Done) - | Step(i,Gcd(c,p),Done) -> normalise_proof id (Step(i,p,Done)) - | Step(i,p,prf) -> - let bds,id,p' = pr_rule_def_cut id p in - let (id,prf) = normalise_proof id prf in - let prf = List.fold_left (fun acc (i,p) -> Step(i, CutPrf p,acc)) - (Step(i,p',prf)) bds in - - (id,prf) - | Enum(i,p1,v,p2,pl) -> - (* Why do I have top-level cuts ? *) -(* let p1 = implicit_cut p1 in - let p2 = implicit_cut p2 in - let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in - (List.fold_left max 0 ids , - Enum(i,p1,v,p2,prfs)) -*) +let eval_op = function + | Eq -> (=/) + | Ge -> (>=/) + | Gt -> (>/) - let bds1,id,p1' = pr_rule_def_cut id (implicit_cut p1) in - let bds2,id,p2' = pr_rule_def_cut id (implicit_cut p2) in - let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in - (List.fold_left max 0 ids , - List.fold_left (fun acc (i,p) -> Step(i, CutPrf p,acc)) - (Enum(i,p1',v,p2',prfs)) (bds2@bds1)) +let string_of_op = function Eq -> "=" | Ge -> ">=" | Gt -> ">" -let normalise_proof id prf = - let res = normalise_proof id prf in - if debug then Printf.printf "normalise_proof %a -> %a" output_proof prf output_proof (snd res) ; - res +let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} = + Printf.fprintf o "%a %s %s" Vect.pp coeffs (string_of_op op) (string_of_num cst) +let opMult o1 o2 = + match o1, o2 with + | Eq , _ | _ , Eq -> Eq + | Ge , _ | _ , Ge -> Ge + | Gt , Gt -> Gt -let add_proof x y = - match x, y with - | Zero , p | p , Zero -> p - | _ -> AddPrf(x,y) +let opAdd o1 o2 = + match o1, o2 with + | Eq , x | x , Eq -> x + | Gt , x | x , Gt -> Gt + | Ge , Ge -> Ge -let mul_proof c p = - match sign_big_int c with - | 0 -> Zero (* This is likely to be a bug *) - | -1 -> MulC(([],Big_int c),p) (* [p] should represent an equality *) - | 1 -> - if eq_big_int c unit_big_int - then p - else MulPrf(Cst c,p) - | _ -> assert false -let mul_proof_ext (p,c) prf = - match p with - | [] -> mul_proof (numerator c) prf - | _ -> MulC((p,c),prf) - +module LinPoly = struct + (** A linear polynomial a0 + a1.x1 + ... + an.xn + By convention, the constant a0 is the coefficient of the variable 0. + *) -module LinPoly = -struct - type t = Vect.t * num + type t = Vect.t - module MonT = - struct + module MonT = struct module MonoMap = Map.Make(Monomial) module IntMap = Map.Make(Int) - + (** A hash table might be preferable but requires a hash function. *) let (index_of_monomial : int MonoMap.t ref) = ref (MonoMap.empty) let (monomial_of_index : Monomial.t IntMap.t ref) = ref (IntMap.empty) let fresh = ref 0 - let clear () = + let clear () = index_of_monomial := MonoMap.empty; - monomial_of_index := IntMap.empty ; + monomial_of_index := IntMap.empty ; fresh := 0 - let register m = + let register m = try - MonoMap.find m !index_of_monomial - with Not_found -> - begin - let res = !fresh in - index_of_monomial := MonoMap.add m res !index_of_monomial ; - monomial_of_index := IntMap.add res m !monomial_of_index ; - incr fresh ; res - end + MonoMap.find m !index_of_monomial + with Not_found -> + begin + let res = !fresh in + index_of_monomial := MonoMap.add m res !index_of_monomial ; + monomial_of_index := IntMap.add res m !monomial_of_index ; + incr fresh ; res + end let retrieve i = IntMap.find i !monomial_of_index + let _ = register Monomial.const - end + end - let normalise (v,c) = - (List.sort (fun x y -> Int.compare (fst x) (fst y)) v , c) + let var v = Vect.set (MonT.register (Monomial.var v)) (Int 1) Vect.null + let of_monomial m = + let v = MonT.register m in + Vect.set v (Int 1) Vect.null - let output_mon o (x,v) = - Printf.fprintf o "%s.%a +" (string_of_num v) Monomial.pp (MonT.retrieve x) + let linpol_of_pol p = + Poly.fold + (fun mon num vct -> + let vr = MonT.register mon in + Vect.set vr num vct) p Vect.null + let pol_of_linpol v = + Vect.fold (fun p vr n -> Poly.add (MonT.retrieve vr) n p) (Poly.constant (Int 0)) v + let coq_poly_of_linpol cst p = - let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} = - Printf.fprintf o "%a %s %s" (pp_list output_mon) coeffs (string_of_op op) (string_of_num cst) + let pol_of_mon m = + Monomial.fold (fun x v p -> Mc.PEmul(Mc.PEpow(Mc.PEX(CamlToCoq.positive x),CamlToCoq.n v),p)) m (Mc.PEc (cst (Int 1))) in + Vect.fold (fun acc x v -> + let mn = MonT.retrieve x in + Mc.PEadd(Mc.PEmul(Mc.PEc (cst v), pol_of_mon mn),acc)) (Mc.PEc (cst (Int 0))) p + let pp_var o vr = + try + Monomial.pp o (MonT.retrieve vr) (* this is a non-linear variable *) + with Not_found -> Printf.fprintf o "v%i" vr + + + let pp o p = Vect.pp_gen pp_var o p + + let constant c = + if sign_num c = 0 + then Vect.null + else Vect.set 0 c Vect.null + + + let is_linear p = + Vect.for_all (fun v _ -> + let mn = (MonT.retrieve v) in + Monomial.is_var mn || Monomial.is_const mn) p + + + let factorise x p = + let (px,cx) = Poly.factorise x (pol_of_linpol p) in + (linpol_of_pol px, linpol_of_pol cx) + + + let is_linear_for x p = + let (a,b) = factorise x p in + Vect.is_constant a + + let search_linear p l = + + Vect.find (fun x v -> + if p v + then + let x' = MonT.retrieve x in + match Monomial.get_var x' with + | None -> None + | Some x -> if is_linear_for x l + then Some x + else None + else None) l + + + let search_all_linear p l = + Vect.fold (fun acc x v -> + if p v + then + let x' = MonT.retrieve x in + match Monomial.get_var x' with + | None -> acc + | Some x -> + if is_linear_for x l + then x::acc + else acc + else acc) [] l + + + let product p1 p2 = + linpol_of_pol (Poly.product (pol_of_linpol p1) (pol_of_linpol p2)) + + let addition p1 p2 = Vect.add p1 p2 + + let variables p = Vect.fold + (fun acc v _ -> + ISet.union (Monomial.variables (MonT.retrieve v)) acc) ISet.empty p + + + let pp_goal typ o l = + let vars = List.fold_left (fun acc p -> ISet.union acc (variables (fst p))) ISet.empty l in + let pp_vars o i = ISet.iter (fun v -> Printf.fprintf o "(x%i : %s) " v typ) vars in + + Printf.fprintf o "forall %a\n" pp_vars vars ; + List.iteri (fun i (p,op) -> Printf.fprintf o "(H%i : %a %s 0)\n" i pp p (string_of_op op)) l; + Printf.fprintf o ", False\n" - let linpol_of_pol p = - let (v,c) = - Poly.fold - (fun mon num (vct,cst) -> - if Monomial.is_const mon then (vct,num) - else - let vr = MonT.register mon in - ((vr,num)::vct,cst)) p ([], Int 0) in - normalise (v,c) - let mult v m (vect,c) = - if Monomial.is_const m - then - (Vect.mul v vect, v <*> c) - else - if sign_num v <> 0 - then - let hd = - if sign_num c <> 0 - then [MonT.register m,v <*> c] - else [] in - - let vect = hd @ (List.map (fun (x,n) -> - let x = MonT.retrieve x in - let x_m = MonT.register (Monomial.prod m x) in - (x_m, v <*> n)) vect ) in - normalise (vect , Int 0) - else ([],Int 0) - let mult v m (vect,c) = - let (vect',c') = mult v m (vect,c) in - if debug then - Printf.printf "mult %s %a (%a,%s) -> (%a,%s)\n" (string_of_num v) Monomial.pp m - (pp_list output_mon) vect (string_of_num c) - (pp_list output_mon) vect' (string_of_num c') ; - (vect',c') + let collect_square p = + Vect.fold (fun acc v _ -> + let m = (MonT.retrieve v) in + match Monomial.sqrt m with + | None -> acc + | Some s -> MonMap.add s m acc + ) MonMap.empty p - let make_lin_pol v mon = - if Monomial.is_const mon - then [] , v - else [MonT.register mon, v],Int 0 - +end + +let output_nlin_cstr o {coeffs = coeffs ; op = op ; cst = cst} = + let p = LinPoly.pol_of_linpol coeffs in + + Printf.fprintf o "%a %s %s" Poly.pp p (string_of_op op) (string_of_num cst) + + +module ProofFormat = struct + open Big_int + + type prf_rule = + | Annot of string * prf_rule + | Hyp of int + | Def of int + | Cst of Num.num + | Zero + | Square of Vect.t + | MulC of Vect.t * prf_rule + | Gcd of Big_int.big_int * prf_rule + | MulPrf of prf_rule * prf_rule + | AddPrf of prf_rule * prf_rule + | CutPrf of prf_rule + + type proof = + | Done + | Step of int * prf_rule * proof + | Enum of int * prf_rule * Vect.t * prf_rule * proof list + + + let rec output_prf_rule o = function + | Annot(s,p) -> Printf.fprintf o "(%a)@%s" output_prf_rule p s + | Hyp i -> Printf.fprintf o "Hyp %i" i + | Def i -> Printf.fprintf o "Def %i" i + | Cst c -> Printf.fprintf o "Cst %s" (string_of_num c) + | Zero -> Printf.fprintf o "Zero" + | Square s -> Printf.fprintf o "(%a)^2" Poly.pp (LinPoly.pol_of_linpol s) + | MulC(p,pr) -> Printf.fprintf o "(%a) * %a" Poly.pp (LinPoly.pol_of_linpol p) output_prf_rule pr + | MulPrf(p1,p2) -> Printf.fprintf o "%a * %a" output_prf_rule p1 output_prf_rule p2 + | AddPrf(p1,p2) -> Printf.fprintf o "%a + %a" output_prf_rule p1 output_prf_rule p2 + | CutPrf(p) -> Printf.fprintf o "[%a]" output_prf_rule p + | Gcd(c,p) -> Printf.fprintf o "(%a)/%s" output_prf_rule p (string_of_big_int c) + + let rec output_proof o = function + | Done -> Printf.fprintf o "." + | Step(i,p,pf) -> Printf.fprintf o "%i:= %a ; %a" i output_prf_rule p output_proof pf + | Enum(i,p1,v,p2,pl) -> Printf.fprintf o "%i{%a<=%a<=%a}%a" i + output_prf_rule p1 Vect.pp v output_prf_rule p2 + (pp_list ";" output_proof) pl + + let rec pr_rule_max_id = function + | Annot(_,p) -> pr_rule_max_id p + | Hyp i | Def i -> i + | Cst _ | Zero | Square _ -> -1 + | MulC(_,p) | CutPrf p | Gcd(_,p) -> pr_rule_max_id p + | MulPrf(p1,p2)| AddPrf(p1,p2) -> max (pr_rule_max_id p1) (pr_rule_max_id p2) + + let rec proof_max_id = function + | Done -> -1 + | Step(i,pr,prf) -> max i (max (pr_rule_max_id pr) (proof_max_id prf)) + | Enum(i,p1,_,p2,l) -> + let m = max (pr_rule_max_id p1) (pr_rule_max_id p2) in + List.fold_left (fun i prf -> max i (proof_max_id prf)) (max i m) l + + + let rec pr_rule_def_cut id = function + | Annot(_,p) -> pr_rule_def_cut id p + | MulC(p,prf) -> + let (bds,id',prf') = pr_rule_def_cut id prf in + (bds, id', MulC(p,prf')) + | MulPrf(p1,p2) -> + let (bds1,id,p1) = pr_rule_def_cut id p1 in + let (bds2,id,p2) = pr_rule_def_cut id p2 in + (bds2@bds1,id,MulPrf(p1,p2)) + | AddPrf(p1,p2) -> + let (bds1,id,p1) = pr_rule_def_cut id p1 in + let (bds2,id,p2) = pr_rule_def_cut id p2 in + (bds2@bds1,id,AddPrf(p1,p2)) + | CutPrf p -> + let (bds,id,p) = pr_rule_def_cut id p in + ((id,p)::bds,id+1,Def id) + | Gcd(c,p) -> + let (bds,id,p) = pr_rule_def_cut id p in + ((id,p)::bds,id+1,Def id) + | Square _|Cst _|Def _|Hyp _|Zero as x -> ([],id,x) + + + (* Do not define top-level cuts *) + let pr_rule_def_cut id = function + | CutPrf p -> + let (bds,ids,p') = pr_rule_def_cut id p in + bds,ids, CutPrf p' + | p -> pr_rule_def_cut id p + + + let rec implicit_cut p = + match p with + | CutPrf p -> implicit_cut p + | _ -> p + + + let rec pr_rule_collect_hyps pr = + match pr with + | Annot(_,pr) -> pr_rule_collect_hyps pr + | Hyp i | Def i -> ISet.add i ISet.empty + | Cst _ | Zero | Square _ -> ISet.empty + | MulC(_,pr) | Gcd(_,pr)| CutPrf pr -> pr_rule_collect_hyps pr + | MulPrf(p1,p2) | AddPrf(p1,p2) -> ISet.union (pr_rule_collect_hyps p1) (pr_rule_collect_hyps p2) + + let simplify_proof p = + let rec simplify_proof p = + match p with + | Done -> (Done, ISet.empty) + | Step(i,pr,Done) -> (p, ISet.add i (pr_rule_collect_hyps pr)) + | Step(i,pr,prf) -> + let (prf',hyps) = simplify_proof prf in + if not (ISet.mem i hyps) + then (prf',hyps) + else + (Step(i,pr,prf'), ISet.add i (ISet.union (pr_rule_collect_hyps pr) hyps)) + | Enum(i,p1,v,p2,pl) -> + let (pl,hl) = List.split (List.map simplify_proof pl) in + let hyps = List.fold_left ISet.union ISet.empty hl in + (Enum(i,p1,v,p2,pl),ISet.add i (ISet.union (ISet.union (pr_rule_collect_hyps p1) (pr_rule_collect_hyps p2)) hyps)) in + fst (simplify_proof p) + + + let rec normalise_proof id prf = + match prf with + | Done -> (id,Done) + | Step(i,Gcd(c,p),Done) -> normalise_proof id (Step(i,p,Done)) + | Step(i,p,prf) -> + let bds,id,p' = pr_rule_def_cut id p in + let (id,prf) = normalise_proof id prf in + let prf = List.fold_left (fun acc (i,p) -> Step(i, CutPrf p,acc)) + (Step(i,p',prf)) bds in + + (id,prf) + | Enum(i,p1,v,p2,pl) -> + (* Why do I have top-level cuts ? *) + (* let p1 = implicit_cut p1 in + let p2 = implicit_cut p2 in + let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in + (List.fold_left max 0 ids , + Enum(i,p1,v,p2,prfs)) + *) + + let bds1,id,p1' = pr_rule_def_cut id (implicit_cut p1) in + let bds2,id,p2' = pr_rule_def_cut id (implicit_cut p2) in + let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in + (List.fold_left max 0 ids , + List.fold_left (fun acc (i,p) -> Step(i, CutPrf p,acc)) + (Enum(i,p1',v,p2',prfs)) (bds2@bds1)) + + + let normalise_proof id prf = + let prf = simplify_proof prf in + let res = normalise_proof id prf in + if debug then Printf.printf "normalise_proof %a -> %a" output_proof prf output_proof (snd res) ; + res + + + + let add_proof x y = + match x, y with + | Zero , p | p , Zero -> p + | _ -> AddPrf(x,y) + + + let mul_cst_proof c p = + match sign_num c with + | 0 -> Zero (* This is likely to be a bug *) + | -1 -> MulC(LinPoly.constant c,p) (* [p] should represent an equality *) + | 1 -> + if eq_num (Int 1) c + then p + else MulPrf(Cst c,p) + | _ -> assert false + let mul_proof p1 p2 = + match p1 , p2 with + | Zero , _ | _ , Zero -> Zero + | Cst (Int 1) , p | p , Cst (Int 1) -> p + | _ , _ -> MulPrf(p1,p2) - let xpivot_eq (c,prf) x v (c',prf') = - if debug then Printf.printf "xpivot_eq {%a} %a %s {%a}\n" - output_cstr c - Monomial.pp (MonT.retrieve x) - (string_of_num v) output_cstr c' ; + let proof_of_farkas env vect = + Vect.fold (fun prf x n -> + add_proof (mul_cst_proof n (IMap.find x env)) prf) Zero vect - let {coeffs = coeffs ; op = op ; cst = cst} = c' in - let m = MonT.retrieve x in - let apply_pivot (vqn,q,n) (c',prf') = - (* Morally, we have (Vect.get (q*x^n) c'.coeffs) = vmn with n >=0 *) + module Env = struct - let cc' = abs_num v in - let cc_num = Int (- (sign_num v)) <*> vqn in - let cc_mon = Monomial.prod q (Monomial.exp m (n-1)) in + let rec string_of_int_list l = + match l with + | [] -> "" + | i::l -> Printf.sprintf "%i,%s" i (string_of_int_list l) - let (c_coeff,c_cst) = mult cc_num cc_mon (c.coeffs, minus_num c.cst) in - - let c' = {coeffs = Vect.add (Vect.mul cc' c'.coeffs) c_coeff ; op = op ; cst = (minus_num c_cst) <+> (cc' <*> c'.cst)} in - let prf' = add_proof - (mul_proof_ext (make_lin_pol cc_num cc_mon) prf) - (mul_proof (numerator cc') prf') in - if debug then Printf.printf "apply_pivot -> {%a}\n" output_cstr c' ; - (c',prf') in + let id_of_hyp hyp l = + let rec xid_of_hyp i l' = + match l' with + | [] -> failwith (Printf.sprintf "id_of_hyp %i %s" hyp (string_of_int_list l)) + | hyp'::l' -> if Pervasives.(=) hyp hyp' then i else xid_of_hyp (i+1) l' in + xid_of_hyp 0 l + end + + let cmpl_prf_rule norm (cst:num-> 'a) env prf = + let rec cmpl = + function + | Annot(s,p) -> cmpl p + | Hyp i | Def i -> Mc.PsatzIn (CamlToCoq.nat (Env.id_of_hyp i env)) + | Cst i -> Mc.PsatzC (cst i) + | Zero -> Mc.PsatzZ + | MulPrf(p1,p2) -> Mc.PsatzMulE(cmpl p1, cmpl p2) + | AddPrf(p1,p2) -> Mc.PsatzAdd(cmpl p1 , cmpl p2) + | MulC(lp,p) -> let lp = norm (LinPoly.coq_poly_of_linpol cst lp) in + Mc.PsatzMulC(lp,cmpl p) + | Square lp -> Mc.PsatzSquare (norm (LinPoly.coq_poly_of_linpol cst lp)) + | _ -> failwith "Cuts should already be compiled" in + cmpl prf + + + + + let cmpl_prf_rule_z env r = cmpl_prf_rule Mc.normZ (fun x -> CamlToCoq.bigint (numerator x)) env r + + let rec cmpl_proof env = function + | Done -> Mc.DoneProof + | Step(i,p,prf) -> + begin + match p with + | CutPrf p' -> + Mc.CutProof(cmpl_prf_rule_z env p', cmpl_proof (i::env) prf) + | _ -> Mc.RatProof(cmpl_prf_rule_z env p,cmpl_proof (i::env) prf) + end + | Enum(i,p1,_,p2,l) -> + Mc.EnumProof(cmpl_prf_rule_z env p1,cmpl_prf_rule_z 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 + cmpl_proof env prf + + let rec eval_prf_rule env = function + | Annot(s,p) -> eval_prf_rule env p + | Hyp i | Def i -> env i + | Cst n -> (Vect.set 0 n Vect.null, + match Num.compare_num n (Int 0) with + | 0 -> Ge + | 1 -> Gt + | _ -> failwith "eval_prf_rule : negative constant" + ) + | Zero -> (Vect.null, Ge) + | Square v -> (LinPoly.product v v,Ge) + | MulC(v, p) -> + let (p1,o) = eval_prf_rule env p in + begin match o with + | Eq -> (LinPoly.product v p1,Eq) + | _ -> + Printf.fprintf stdout "MulC(%a,%a) invalid 2d arg %a %s" Vect.pp v output_prf_rule p Vect.pp p1 (string_of_op o); + failwith "eval_prf_rule : not an equality" + end + | Gcd(g,p) -> let (v,op) = eval_prf_rule env p in + (Vect.div (Big_int g) v, op) + | MulPrf(p1,p2) -> + let (v1,o1) = eval_prf_rule env p1 in + let (v2,o2) = eval_prf_rule env p2 in + (LinPoly.product v1 v2, opMult o1 o2) + | AddPrf(p1,p2) -> + let (v1,o1) = eval_prf_rule env p1 in + let (v2,o2) = eval_prf_rule env p2 in + (LinPoly.addition v1 v2, opAdd o1 o2) + | CutPrf p -> eval_prf_rule env p + + + let is_unsat (p,o) = + let (c,r) = Vect.decomp_cst p in + if Vect.is_null r + then not (eval_op o c (Int 0)) + else false + + let rec eval_proof env p = + match p with + | Done -> failwith "Proof is not finished" + | Step(i, prf, rst) -> + let (p,o) = eval_prf_rule (fun i -> IMap.find i env) prf in + if is_unsat (p,o) then true + else + if Pervasives.(=) rst Done + then + begin + Printf.fprintf stdout "Last inference %a %s\n" LinPoly.pp p (string_of_op o); + false + end + else eval_proof (IMap.add i (p,o) env) rst + | Enum(i,r1,v,r2,l) -> let _ = eval_prf_rule (fun i -> IMap.find i env) r1 in + let _ = eval_prf_rule (fun i -> IMap.find i env) r2 in + (* Should check bounds *) + failwith "Not implemented" + +end + +module WithProof = struct + + type t = ((LinPoly.t * op) * ProofFormat.prf_rule) + + let annot s (p,prf) = (p, ProofFormat.Annot(s,prf)) + + let output o ((lp,op),prf) = + Printf.fprintf o "%a %s 0 by %a\n" LinPoly.pp lp (string_of_op op) ProofFormat.output_prf_rule prf - let cmp (q,n) (q',n') = - if n < n' then -1 - else if n = n' then Monomial.compare q q' - else 1 in + exception InvalidProof - - let find_pivot (c',prf') = - let (v,q,n) = List.fold_left - (fun (v,q,n) (x,v') -> - let x = MonT.retrieve x in - let (q',n') = Monomial.div x m in - if cmp (q,n) (q',n') = -1 then (v',q',n') else (v,q,n)) (Int 0, Monomial.const,0) c'.coeffs in - if n > 0 then Some (v,q,n) else None in + let zero = ((Vect.null,Eq), ProofFormat.Zero) - let rec pivot (q,n) (c',prf') = - match find_pivot (c',prf') with - | None -> (c',prf') - | Some(v,q',n') -> - if cmp (q',n') (q,n) = -1 - then pivot (q',n') (apply_pivot (v,q',n') (c',prf')) - else (c',prf') in - pivot (Monomial.const,max_int) (c',prf') + let of_cstr (c,prf) = + (Vect.set 0 (Num.minus_num (c.cst)) c.coeffs,c.op), prf + let product : t -> t -> t = fun ((p1,o1),prf1) ((p2,o2),prf2) -> + ((LinPoly.product p1 p2 , opMult o1 o2), ProofFormat.mul_proof prf1 prf2) - let pivot_eq x (c,prf) = - match Vect.get x c.coeffs with - | None -> (fun x -> None) - | Some v -> fun cp' -> Some (xpivot_eq (c,prf) x v cp') + let addition : t -> t -> t = fun ((p1,o1),prf1) ((p2,o2),prf2) -> + ((Vect.add p1 p2, opAdd o1 o2), ProofFormat.add_proof prf1 prf2) + let mult p ((p1,o1),prf1) = + match o1 with + | Eq -> ((LinPoly.product p p1,o1), ProofFormat.MulC(p, prf1)) + | Gt| Ge -> let (n,r) = Vect.decomp_cst p in + if Vect.is_null r && n >/ Int 0 + then ((LinPoly.product p p1, o1), ProofFormat.mul_cst_proof n prf1) + else raise InvalidProof + + + let cutting_plane ((p,o),prf) = + let (c,p') = Vect.decomp_cst p in + let g = (Vect.gcd p') in + if (Big_int.eq_big_int Big_int.unit_big_int g) || c =/ Int 0 || + not (Big_int.eq_big_int (denominator c) Big_int.unit_big_int) + then None (* Nothing to do *) + else + let c1 = c // (Big_int g) in + let c1' = Num.floor_num c1 in + if c1 =/ c1' + then None + else + match o with + | Eq -> Some ((Vect.set 0 (Int (-1)) Vect.null,Eq), ProofFormat.Gcd(g,prf)) + | Gt -> failwith "cutting_plane ignore strict constraints" + | Ge -> + (* This is a non-trivial common divisor *) + Some ((Vect.set 0 c1' (Vect.div (Big_int g) p),o),ProofFormat.Gcd(g, prf)) + + + let construct_sign p = + let (c,p') = Vect.decomp_cst p in + if Vect.is_null p' + then + Some (begin match sign_num c with + | 0 -> (true, Eq, ProofFormat.Zero) + | 1 -> (true,Gt, ProofFormat.Cst c) + | _ (*-1*) -> (false,Gt, ProofFormat.Cst (minus_num c)) + end) + else None + + + let get_sign l p = + match construct_sign p with + | None -> begin + try + let ((p',o),prf) = + List.find (fun ((p',o),prf) -> Vect.equal p p') l in + Some (true,o,prf) + with Not_found -> + let p = Vect.uminus p in + try + let ((p',o),prf) = List.find (fun ((p',o),prf) -> Vect.equal p p') l in + Some (false,o,prf) + with Not_found -> None + end + | Some s -> Some s + + + let mult_sign : bool -> t -> t = fun b ((p,o),prf) -> + if b then ((p,o),prf) + else ((Vect.uminus p,o),prf) + + + let rec linear_pivot sys ((lp1,op1), prf1) x ((lp2,op2), prf2) = + + (* lp1 = a1.x + b1 *) + let (a1,b1) = LinPoly.factorise x lp1 in + + (* lp2 = a2.x + b2 *) + let (a2,b2) = LinPoly.factorise x lp2 in + + if Vect.is_null a2 + then (* We are done *) + Some ((lp2,op2),prf2) + else + match op1,op2 with + | Eq , (Ge|Gt) -> begin + match get_sign sys a1 with + | None -> None (* Impossible to pivot without sign information *) + | Some(b,o,prf) -> + let sa1 = mult_sign b ((a1,o),prf) in + let sa2 = if b then (Vect.uminus a2) else a2 in + + let ((lp2,op2),prf2) = + addition (product sa1 ((lp2,op2),prf2)) + (mult sa2 ((lp1,op1),prf1)) in + linear_pivot sys ((lp1,op1),prf1) x ((lp2,op2),prf2) + + end + | Eq , Eq -> + let ((lp2,op2),prf2) = addition (mult a1 ((lp2,op2),prf2)) + (mult (Vect.uminus a2) ((lp1,op1),prf1)) in + linear_pivot sys ((lp1,op1),prf1) x ((lp2,op2),prf2) + + | (Ge | Gt) , (Ge| Gt) -> begin + match get_sign sys a1 , get_sign sys a2 with + | Some(b1,o1,p1) , Some(b2,o2,p2) -> + if b1 <> b2 + then + let ((lp2,op2),prf2) = + addition (product (mult_sign b1 ((a1,o1), p1)) ((lp2,op2),prf2)) + (product (mult_sign b2 ((a2,o2), p2)) ((lp1,op1),prf1)) in + linear_pivot sys ((lp1,op1),prf1) x ((lp2,op2),prf2) + else None + | _ -> None + end + | (Ge|Gt) , Eq -> failwith "pivot: equality as second argument" end + + +(* Local Variables: *) +(* coding: utf-8 *) +(* End: *) |