(* ========================================================================= *) (* Real quantifier elimination (using Cohen-Hormander algorithm). *) (* *) (* Copyright (c) 2003-2007, John Harrison. (See "LICENSE.txt" for details.) *) (* ========================================================================= *) (* ------------------------------------------------------------------------- *) (* Formal derivative of polynomial. *) (* ------------------------------------------------------------------------- *) let rec poly_diffn x n p = match p with Fn("+",[c; Fn("*",[y; q])]) when y = x -> Fn("+",[poly_cmul(Int n) c; Fn("*",[x; poly_diffn x (n+1) q])]) | _ -> poly_cmul(Int n) p;; let poly_diff vars p = match p with Fn("+",[c; Fn("*",[Var x; q])]) when x = hd vars -> poly_diffn (Var x) 1 q | _ -> zero;; (* ------------------------------------------------------------------------- *) (* Evaluate a quantifier-free formula given a sign matrix row for its polys. *) (* ------------------------------------------------------------------------- *) let rel_signs = ["=",[Zero]; "<=",[Zero;Negative]; ">=",[Zero;Positive]; "<",[Negative]; ">",[Positive]];; let testform pmat fm = eval fm (fun (R(a,[p;z])) -> mem (assoc p pmat) (assoc a rel_signs));; (* ------------------------------------------------------------------------- *) (* Infer sign of p(x) at points from corresponding qi(x) with pi(x) = 0 *) (* ------------------------------------------------------------------------- *) let inferpsign (pd,qd) = try let i = index Zero pd in el i qd :: pd with Failure _ -> Nonzero :: pd;; (* ------------------------------------------------------------------------- *) (* Condense subdivision by removing points with no relevant zeros. *) (* ------------------------------------------------------------------------- *) let rec condense ps = match ps with int::pt::other -> let rest = condense other in if mem Zero pt then int::pt::rest else rest | _ -> ps;; (* ------------------------------------------------------------------------- *) (* Infer sign on intervals (use with infinities at end) and split if needed *) (* ------------------------------------------------------------------------- *) let rec inferisign ps = match ps with ((l::ls) as x)::(_::ints)::((r::rs)::xs as pts) -> (match (l,r) with (Zero,Zero) -> failwith "inferisign: inconsistent" | (Nonzero,_) | (_,Nonzero) -> failwith "inferisign: indeterminate" | (Zero,_) -> x::(r::ints)::inferisign pts | (_,Zero) -> x::(l::ints)::inferisign pts | (Negative,Negative) | (Positive,Positive) -> x::(l::ints)::inferisign pts | _ -> x::(l::ints)::(Zero::ints)::(r::ints)::inferisign pts) | _ -> ps;; (* ------------------------------------------------------------------------- *) (* Deduce matrix for p,p1,...,pn from matrix for p',p1,...,pn,q0,...,qn *) (* where qi = rem(p,pi) with p0 = p' *) (* ------------------------------------------------------------------------- *) let dedmatrix cont mat = let l = length (hd mat) / 2 in let mat1 = condense(map (inferpsign ** chop_list l) mat) in let mat2 = [swap true (el 1 (hd mat1))]::mat1@[[el 1 (last mat1)]] in let mat3 = butlast(tl(inferisign mat2)) in cont(condense(map (fun l -> hd l :: tl(tl l)) mat3));; (* ------------------------------------------------------------------------- *) (* Pseudo-division making sure the remainder has the same sign. *) (* ------------------------------------------------------------------------- *) let pdivide_pos vars sgns s p = let a = head vars p and (k,r) = pdivide vars s p in let sgn = findsign sgns a in if sgn = Zero then failwith "pdivide_pos: zero head coefficient" else if sgn = Positive or k mod 2 = 0 then r else if sgn = Negative then poly_neg r else poly_mul vars a r;; (* ------------------------------------------------------------------------- *) (* Case splitting for positive/negative (assumed nonzero). *) (* ------------------------------------------------------------------------- *) let split_sign sgns pol cont = match findsign sgns pol with Nonzero -> let fm = Atom(R(">",[pol; zero])) in Or(And(fm,cont(assertsign sgns (pol,Positive))), And(Not fm,cont(assertsign sgns (pol,Negative)))) | _ -> cont sgns;; let split_trichotomy sgns pol cont_z cont_pn = split_zero sgns pol cont_z (fun s' -> split_sign s' pol cont_pn);; (* ------------------------------------------------------------------------- *) (* Main recursive evaluation of sign matrices. *) (* ------------------------------------------------------------------------- *) let rec casesplit vars dun pols cont sgns = match pols with [] -> matrix vars dun cont sgns | p::ops -> split_trichotomy sgns (head vars p) (if is_constant vars p then delconst vars dun p ops cont else casesplit vars dun (behead vars p :: ops) cont) (if is_constant vars p then delconst vars dun p ops cont else casesplit vars (dun@[p]) ops cont) and delconst vars dun p ops cont sgns = let cont' m = cont(map (insertat (length dun) (findsign sgns p)) m) in casesplit vars dun ops cont' sgns and matrix vars pols cont sgns = if pols = [] then try cont [[]] with Failure _ -> False else let p = hd(sort(decreasing (degree vars)) pols) in let p' = poly_diff vars p and i = index p pols in let qs = let p1,p2 = chop_list i pols in p'::p1 @ tl p2 in let gs = map (pdivide_pos vars sgns p) qs in let cont' m = cont(map (fun l -> insertat i (hd l) (tl l)) m) in casesplit vars [] (qs@gs) (dedmatrix cont') sgns;; (* ------------------------------------------------------------------------- *) (* Now the actual quantifier elimination code. *) (* ------------------------------------------------------------------------- *) let basic_real_qelim vars (Exists(x,p)) = let pols = atom_union (function (R(a,[t;Fn("0",[])])) -> [t] | _ -> []) p in let cont mat = if exists (fun m -> testform (zip pols m) p) mat then True else False in casesplit (x::vars) [] pols cont init_sgns;; let real_qelim = simplify ** evalc ** lift_qelim polyatom (simplify ** evalc) basic_real_qelim;; (* ------------------------------------------------------------------------- *) (* First examples. *) (* ------------------------------------------------------------------------- *) START_INTERACTIVE;; real_qelim <>;; real_qelim <>;; real_qelim <>;; #trace testform;; real_qelim <>;; #untrace testform;; real_qelim < f < a * e) ==> f <= a * k>>;; real_qelim <>;; real_qelim < b^2 >= 4 * a * c>>;; real_qelim < a = 0 /\ (b = 0 ==> c = 0) \/ ~(a = 0) /\ b^2 >= 4 * a * c>>;; (* ------------------------------------------------------------------------- *) (* Termination ordering for group theory completion. *) (* ------------------------------------------------------------------------- *) real_qelim <<1 < 2 /\ (forall x. 1 < x ==> 1 < x^2) /\ (forall x y. 1 < x /\ 1 < y ==> 1 < x * (1 + 2 * y))>>;; END_INTERACTIVE;; let rec grpterm tm = match tm with Fn("*",[s;t]) -> let t2 = Fn("*",[Fn("2",[]); grpterm t]) in Fn("*",[grpterm s; Fn("+",[Fn("1",[]); t2])]) | Fn("i",[t]) -> Fn("^",[grpterm t; Fn("2",[])]) | Fn("1",[]) -> Fn("2",[]) | Var x -> tm;; let grpform (Atom(R("=",[s;t]))) = let fm = generalize(Atom(R(">",[grpterm s; grpterm t]))) in relativize(fun x -> Atom(R(">",[Var x;Fn("1",[])]))) fm;; START_INTERACTIVE;; let eqs = complete_and_simplify ["1"; "*"; "i"] [<<1 * x = x>>; <>; <<(x * y) * z = x * y * z>>];; let fm = list_conj (map grpform eqs);; real_qelim fm;; END_INTERACTIVE;; (* ------------------------------------------------------------------------- *) (* A case where using DNF is an improvement. *) (* ------------------------------------------------------------------------- *) let real_qelim' = simplify ** evalc ** lift_qelim polyatom (dnf ** cnnf (fun x -> x) ** evalc) basic_real_qelim;; real_qelim' < a^2 = b) <=> d^4 = 1>>;; (* ------------------------------------------------------------------------- *) (* Didn't seem worth it in the book, but monicization can help a lot. *) (* Now this is just set as an exercise. *) (* ------------------------------------------------------------------------- *) START_INTERACTIVE;; let rec casesplit vars dun pols cont sgns = match pols with [] -> monicize vars dun cont sgns | p::ops -> split_trichotomy sgns (head vars p) (if is_constant vars p then delconst vars dun p ops cont else casesplit vars dun (behead vars p :: ops) cont) (if is_constant vars p then delconst vars dun p ops cont else casesplit vars (dun@[p]) ops cont) and delconst vars dun p ops cont sgns = let cont' m = cont(map (insertat (length dun) (findsign sgns p)) m) in casesplit vars dun ops cont' sgns and matrix vars pols cont sgns = if pols = [] then try cont [[]] with Failure _ -> False else let p = hd(sort(decreasing (degree vars)) pols) in let p' = poly_diff vars p and i = index p pols in let qs = let p1,p2 = chop_list i pols in p'::p1 @ tl p2 in let gs = map (pdivide_pos vars sgns p) qs in let cont' m = cont(map (fun l -> insertat i (hd l) (tl l)) m) in casesplit vars [] (qs@gs) (dedmatrix cont') sgns and monicize vars pols cont sgns = let mols,swaps = unzip(map monic pols) in let sols = setify mols in let indices = map (fun p -> index p sols) mols in let transform m = map2 (fun sw i -> swap sw (el i m)) swaps indices in let cont' mat = cont(map transform mat) in matrix vars sols cont' sgns;; let basic_real_qelim vars (Exists(x,p)) = let pols = atom_union (function (R(a,[t;Fn("0",[])])) -> [t] | _ -> []) p in let cont mat = if exists (fun m -> testform (zip pols m) p) mat then True else False in casesplit (x::vars) [] pols cont init_sgns;; let real_qelim = simplify ** evalc ** lift_qelim polyatom (simplify ** evalc) basic_real_qelim;; let real_qelim' = simplify ** evalc ** lift_qelim polyatom (dnf ** cnnf (fun x -> x) ** evalc) basic_real_qelim;; END_INTERACTIVE;;