|
(************************************************************************)
(* * The Coq Proof Assistant / The Coq Development Team *)
(* v * INRIA, CNRS and contributors - Copyright 1999-2018 *)
(* <O___,, * (see CREDITS file for the list of authors) *)
(* \VV/ **************************************************************)
(* // * This file is distributed under the terms of the *)
(* * GNU Lesser General Public License Version 2.1 *)
(* * (see LICENSE file for the text of the license) *)
(************************************************************************)
(* *)
(* Micromega: A reflexive tactic using the Positivstellensatz *)
(* *)
(* Frédéric Besson (Irisa/Inria) 2006-20018 *)
(* *)
(************************************************************************)
open Num
module Utils = Mutils
open Utils
module Mc = Micromega
let max_nb_cstr = ref max_int
type var = int
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 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 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
(expected to be fast enough - Coq is probably the bottleneck)
*The new ring contribution is using a sparse Horner representation.
*)
sig
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 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
(** 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
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 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 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
type vector = Vect.t
type cstr = {coeffs : vector ; op : op ; cst : num}
and op = |Eq | Ge | Gt
exception Strict
let is_strict c = Pervasives.(=) c.op Gt
let eval_op = function
| Eq -> (=/)
| Ge -> (>=/)
| Gt -> (>/)
let string_of_op = function Eq -> "=" | Ge -> ">=" | Gt -> ">"
let output_cstr o { coeffs ; op ; 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 opAdd o1 o2 =
match o1, o2 with
| Eq , x | x , Eq -> x
| Gt , x | x , Gt -> Gt
| Ge , Ge -> Ge
module LinPoly = struct
(** A linear polynomial a0 + a1.x1 + ... + an.xn
By convention, the constant a0 is the coefficient of the variable 0.
*)
type t = Vect.t
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 () =
index_of_monomial := MonoMap.empty;
monomial_of_index := IntMap.empty ;
fresh := 0
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
let retrieve i = IntMap.find i !monomial_of_index
let _ = register Monomial.const
end
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 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 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 is_variable p =
let ((x,v),r) = Vect.decomp_fst p in
if Vect.is_null r && v >/ Int 0
then Monomial.get_var (MonT.retrieve x)
else None
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_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 min_list (l:int list) =
match l with
| [] -> None
| e::l -> Some (List.fold_left Pervasives.min e l)
let search_linear p l =
min_list (search_all_linear p 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 of_vect v =
Vect.fold (fun acc v vl -> addition (product (var v) (constant vl)) acc) Vect.null v
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 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
end
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_size = function
| Annot(_,p) -> pr_size p
| Zero| Square _ -> Int 0
| Hyp _ -> Int 1
| Def _ -> Int 1
| Cst n -> n
| Gcd(i, p) -> pr_size p // (Big_int i)
| MulPrf(p1,p2) | AddPrf(p1,p2) -> pr_size p1 +/ pr_size p2
| CutPrf p -> pr_size p
| MulC(v, p) -> pr_size p
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
module OrdPrfRule =
struct
type t = prf_rule
let id_of_constr = function
| Annot _ -> 0
| Hyp _ -> 1
| Def _ -> 2
| Cst _ -> 3
| Zero -> 4
| Square _ -> 5
| MulC _ -> 6
| Gcd _ -> 7
| MulPrf _ -> 8
| AddPrf _ -> 9
| CutPrf _ -> 10
let cmp_pair c1 c2 (x1,x2) (y1,y2) =
match c1 x1 y1 with
| 0 -> c2 x2 y2
| i -> i
let rec compare p1 p2 =
match p1, p2 with
| Annot(s1,p1) , Annot(s2,p2) -> if s1 = s2 then compare p1 p2
else Pervasives.compare s1 s2
| Hyp i , Hyp j -> Pervasives.compare i j
| Def i , Def j -> Pervasives.compare i j
| Cst n , Cst m -> Num.compare_num n m
| Zero , Zero -> 0
| Square v1 , Square v2 -> Vect.compare v1 v2
| MulC(v1,p1) , MulC(v2,p2) -> cmp_pair Vect.compare compare (v1,p1) (v2,p2)
| Gcd(b1,p1) , Gcd(b2,p2) -> cmp_pair Big_int.compare_big_int compare (b1,p1) (b2,p2)
| MulPrf(p1,q1) , MulPrf(p2,q2) -> cmp_pair compare compare (p1,q1) (p2,q2)
| AddPrf(p1,q1) , MulPrf(p2,q2) -> cmp_pair compare compare (p1,q1) (p2,q2)
| CutPrf p , CutPrf p' -> compare p p'
| _ , _ -> Pervasives.compare (id_of_constr p1) (id_of_constr p2)
end
let add_proof x y =
match x, y with
| Zero , p | p , Zero -> p
| _ -> AddPrf(x,y)
let rec mul_cst_proof c p =
match p with
| Annot(s,p) -> Annot(s,mul_cst_proof c p)
| MulC(v,p') -> MulC(Vect.mul c v,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 sMulC v p =
let (c,v') = Vect.decomp_cst v in
if Vect.is_null v' then mul_cst_proof c p
else MulC(v,p)
let mul_proof p1 p2 =
match p1 , p2 with
| Zero , _ | _ , Zero -> Zero
| Cst c , p | p , Cst c -> mul_cst_proof c p
| _ , _ ->
MulPrf(p1,p2)
module PrfRuleMap = Map.Make(OrdPrfRule)
let prf_rule_of_map m =
PrfRuleMap.fold (fun k v acc -> add_proof (sMulC v k) acc) m Zero
let rec dev_prf_rule p =
match p with
| Annot(s,p) -> dev_prf_rule p
| Hyp _ | Def _ | Cst _ | Zero | Square _ -> PrfRuleMap.singleton p (LinPoly.constant (Int 1))
| MulC(v,p) -> PrfRuleMap.map (fun v1 -> LinPoly.product v v1) (dev_prf_rule p)
| AddPrf(p1,p2) -> PrfRuleMap.merge (fun k o1 o2 ->
match o1 , o2 with
| None , None -> None
| None , Some v | Some v, None -> Some v
| Some v1 , Some v2 -> Some (LinPoly.addition v1 v2)) (dev_prf_rule p1) (dev_prf_rule p2)
| MulPrf(p1, p2) ->
begin
let p1' = dev_prf_rule p1 in
let p2' = dev_prf_rule p2 in
let p1'' = prf_rule_of_map p1' in
let p2'' = prf_rule_of_map p2' in
match p1'' with
| Cst c -> PrfRuleMap.map (fun v1 -> Vect.mul c v1) p2'
| _ -> PrfRuleMap.singleton (MulPrf(p1'',p2'')) (LinPoly.constant (Int 1))
end
| _ -> PrfRuleMap.singleton p (LinPoly.constant (Int 1))
let simplify_prf_rule p =
prf_rule_of_map (dev_prf_rule p)
(*
let mul_proof p1 p2 =
let res = mul_proof p1 p2 in
Printf.printf "mul_proof %a %a = %a\n"
output_prf_rule p1 output_prf_rule p2 output_prf_rule res; res
let add_proof p1 p2 =
let res = add_proof p1 p2 in
Printf.printf "add_proof %a %a = %a\n"
output_prf_rule p1 output_prf_rule p2 output_prf_rule res; res
let sMulC v p =
let res = sMulC v p in
Printf.printf "sMulC %a %a = %a\n" Vect.pp v output_prf_rule p output_prf_rule res ;
res
let mul_cst_proof c p =
let res = mul_cst_proof c p in
Printf.printf "mul_cst_proof %s %a = %a\n" (Num.string_of_num c) output_prf_rule p output_prf_rule res ;
res
*)
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
module Env = struct
let rec string_of_int_list l =
match l with
| [] -> ""
| i::l -> Printf.sprintf "%i,%s" i (string_of_int_list l)
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 output_sys o l =
List.iter (Printf.fprintf o "%a\n" output) l
exception InvalidProof
let zero = ((Vect.null,Eq), ProofFormat.Zero)
let const n = ((LinPoly.constant n,Ge), ProofFormat.Cst n)
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 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.sMulC 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"
let linear_pivot sys ((lp1,op1), prf1) x ((lp2,op2), prf2) =
match linear_pivot sys ((lp1,op1), prf1) x ((lp2,op2), prf2) with
| None -> None
| Some (c,p) -> Some(c, ProofFormat.simplify_prf_rule p)
let is_substitution strict ((p,o),prf) =
let pred v = if strict then v =/ Int 1 || v =/ Int (-1) else true in
match o with
| Eq -> LinPoly.search_linear pred p
| _ -> None
let subst1 sys0 =
let (oeq,sys') = extract (is_substitution true) sys0 in
match oeq with
| None -> sys0
| Some(v,pc) ->
match simplify (linear_pivot sys0 pc v) sys' with
| None -> sys0
| Some sys' -> sys'
let subst sys0 =
let elim sys =
let (oeq,sys') = extract (is_substitution true) sys in
match oeq with
| None -> None
| Some(v,pc) -> simplify (linear_pivot sys0 pc v) sys' in
iterate_until_stable elim sys0
let saturate_subst b sys0 =
let select = is_substitution b in
let gen (v,pc) ((c,op),prf) =
if ISet.mem v (LinPoly.variables c)
then linear_pivot sys0 pc v ((c,op),prf)
else None
in
saturate select gen sys0
end
(* Local Variables: *)
(* coding: utf-8 *)
(* End: *)
¤ Dauer der Verarbeitung: 0.31 Sekunden
(vorverarbeitet)
¤
|
Haftungshinweis
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung ist noch experimentell.
|
