Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php on line 456
(* ========================================================================= *) (* Complex quantifier elimination (by simple divisibility a la Tarski). *) (* *) (* Copyright (c) 2003-2007, John Harrison. (See "LICENSE.txt" for details.) *) (* ========================================================================= *) (* ------------------------------------------------------------------------- *) (* Basic arithmetic operations on canonical polynomials. *) (* ------------------------------------------------------------------------- *) let rec poly_add vars pol1 pol2 = match (pol1,pol2) with (Fn("+",[c; Fn("*",[Var x; p])]),Fn("+",[d; Fn("*",[Var y; q])])) -> if earlier vars x y then poly_ladd vars pol2 pol1 else if earlier vars y x then poly_ladd vars pol1 pol2 else let e = poly_add vars c d and r = poly_add vars p q in if r = zero then e else Fn("+",[e; Fn("*",[Var x; r])]) | (_,Fn("+",_)) -> poly_ladd vars pol1 pol2 | (Fn("+",_),pol2) -> poly_ladd vars pol2 pol1 | _ -> numeral2 (+/) pol1 pol2 and poly_ladd vars = fun pol1 (Fn("+",[d; Fn("*",[Var y; q])])) -> Fn("+",[poly_add vars pol1 d; Fn("*",[Var y; q])]);; let rec poly_neg = function (Fn("+",[c; Fn("*",[Var x; p])])) -> Fn("+",[poly_neg c; Fn("*",[Var x; poly_neg p])]) | n -> numeral1 minus_num n;; let poly_sub vars p q = poly_add vars p (poly_neg q);; let rec poly_mul vars pol1 pol2 = match (pol1,pol2) with (Fn("+",[c; Fn("*",[Var x; p])]),Fn("+",[d; Fn("*",[Var y; q])])) -> if earlier vars x y then poly_lmul vars pol2 pol1 else poly_lmul vars pol1 pol2 | (Fn("0",[]),_) | (_,Fn("0",[])) -> zero | (_,Fn("+",_)) -> poly_lmul vars pol1 pol2 | (Fn("+",_),_) -> poly_lmul vars pol2 pol1 | _ -> numeral2 ( */ ) pol1 pol2 and poly_lmul vars = fun pol1 (Fn("+",[d; Fn("*",[Var y; q])])) -> poly_add vars (poly_mul vars pol1 d) (Fn("+",[zero; Fn("*",[Var y; poly_mul vars pol1 q])]));; let poly_pow vars p n = funpow n (poly_mul vars p) (Fn("1",[]));; let poly_div vars p q = poly_mul vars p (numeral1((//) (Int 1)) q);; let poly_var x = Fn("+",[zero; Fn("*",[Var x; Fn("1",[])])]);; (* ------------------------------------------------------------------------- *) (* Convert term into canonical polynomial representative. *) (* ------------------------------------------------------------------------- *) let rec polynate vars tm = match tm with Var x -> poly_var x | Fn("-",[t]) -> poly_neg (polynate vars t) | Fn("+",[s;t]) -> poly_add vars (polynate vars s) (polynate vars t) | Fn("-",[s;t]) -> poly_sub vars (polynate vars s) (polynate vars t) | Fn("*",[s;t]) -> poly_mul vars (polynate vars s) (polynate vars t) | Fn("/",[s;t]) -> poly_div vars (polynate vars s) (polynate vars t) | Fn("^",[p;Fn(n,[])]) -> poly_pow vars (polynate vars p) (int_of_string n) | _ -> if is_numeral tm then tm else failwith "lint: unknown term";; (* ------------------------------------------------------------------------- *) (* Do likewise for atom so the RHS is zero. *) (* ------------------------------------------------------------------------- *) let polyatom vars fm = match fm with Atom(R(a,[s;t])) -> Atom(R(a,[polynate vars (Fn("-",[s;t]));zero])) | _ -> failwith "polyatom: not an atom";; (* ------------------------------------------------------------------------- *) (* Sanity check. *) (* ------------------------------------------------------------------------- *) START_INTERACTIVE;; polyatom ["w"; "x"; "y"; "z"] <<((w + x)^4 + (w + y)^4 + (w + z)^4 + (x + y)^4 + (x + z)^4 + (y + z)^4 + (w - x)^4 + (w - y)^4 + (w - z)^4 + (x - y)^4 + (x - z)^4 + (y - z)^4) / 6 = (w^2 + x^2 + y^2 + z^2)^2>>;; END_INTERACTIVE;; (* ------------------------------------------------------------------------- *) (* Useful utility functions for polynomial terms. *) (* ------------------------------------------------------------------------- *) let rec coefficients vars = function Fn("+",[c; Fn("*",[Var x; q])]) when x = hd vars -> c::(coefficients vars q) | p -> [p];; let degree vars p = length(coefficients vars p) - 1;; let is_constant vars p = degree vars p = 0;; let head vars p = last(coefficients vars p);; let rec behead vars = function Fn("+",[c; Fn("*",[Var x; p])]) when x = hd vars -> let p' = behead vars p in if p' = zero then c else Fn("+",[c; Fn("*",[Var x; p'])]) | _ -> zero;; (* ------------------------------------------------------------------------- *) (* Get the constant multiple of the "maximal" monomial (implicit lex order) *) (* ------------------------------------------------------------------------- *) let rec poly_cmul k p = match p with Fn("+",[c; Fn("*",[Var x; q])]) -> Fn("+",[poly_cmul k c; Fn("*",[Var x; poly_cmul k q])]) | _ -> numeral1 (fun m -> k */ m) p;; let rec headconst p = match p with Fn("+",[c; Fn("*",[Var x; q])]) -> headconst q | Fn(n,[]) -> dest_numeral p;; (* ------------------------------------------------------------------------- *) (* Make a polynomial monic and return negativity flag for head constant *) (* ------------------------------------------------------------------------- *) let monic p = let h = headconst p in if h =/ Int 0 then p,false else poly_cmul (Int 1 // h) p,h pdivide_aux vars (head vars p) (degree vars p) p 0 s;; (* ------------------------------------------------------------------------- *) (* Datatype of signs. *) (* ------------------------------------------------------------------------- *) type sign = Zero | Nonzero | Positive | Negative;; let swap swf s = if not swf then s else match s with Positive -> Negative | Negative -> Positive | _ -> s;; (* ------------------------------------------------------------------------- *) (* Lookup and asserting of polynomial sign, modulo constant multiples. *) (* Note that we are building in a characteristic-zero assumption here. *) (* ------------------------------------------------------------------------- *) let findsign sgns p = try let p',swf = monic p in swap swf (assoc p' sgns) with Failure _ -> failwith "findsign";; let assertsign sgns (p,s) = if p = zero then if s = Zero then sgns else failwith "assertsign" else let p',swf = monic p in let s' = swap swf s in let s0 = try assoc p' sgns with Failure _ -> s' in if s' = s0 or s0 = Nonzero & (s' = Positive or s' = Negative) then (p',s')::(subtract sgns [p',s0]) else failwith "assertsign";; (* ------------------------------------------------------------------------- *) (* Deduce or case-split over zero status of polynomial. *) (* ------------------------------------------------------------------------- *) let split_zero sgns pol cont_z cont_n = try let z = findsign sgns pol in (if z = Zero then cont_z else cont_n) sgns with Failure "findsign" -> let eq = Atom(R("=",[pol; zero])) in Or(And(eq,cont_z (assertsign sgns (pol,Zero))), And(Not eq,cont_n (assertsign sgns (pol,Nonzero))));; (* ------------------------------------------------------------------------- *) (* Whether a polynomial is nonzero in a context. *) (* ------------------------------------------------------------------------- *) let poly_nonzero vars sgns pol = let cs = coefficients vars pol in let dcs,ucs = partition (can (findsign sgns)) cs in if exists (fun p -> findsign sgns p <> Zero) dcs then True else if ucs = [] then False else end_itlist mk_or (map (fun p -> Not(mk_eq p zero)) ucs);; (* ------------------------------------------------------------------------- *) (* Non-divisibility of q by p. *) (* ------------------------------------------------------------------------- *) let rec poly_nondiv vars sgns p s = let _,r = pdivide vars s p in poly_nonzero vars sgns r;; (* ------------------------------------------------------------------------- *) (* Main reduction for exists x. all eqs = 0 and all neqs =/= 0, in context. *) (* ------------------------------------------------------------------------- *) let rec cqelim vars (eqs,neqs) sgns = try let c = find (is_constant vars) eqs in (try let sgns' = assertsign sgns (c,Zero) and eqs' = subtract eqs [c] in And(mk_eq c zero,cqelim vars (eqs',neqs) sgns') with Failure "assertsign" -> False) with Failure _ -> if eqs = [] then list_conj(map (poly_nonzero vars sgns) neqs) else let n = end_itlist min (map (degree vars) eqs) in let p = find (fun p -> degree vars p = n) eqs in let oeqs = subtract eqs [p] in split_zero sgns (head vars p) (cqelim vars (behead vars p::oeqs,neqs)) (fun sgns' -> let cfn s = snd(pdivide vars s p) in if oeqs <> [] then cqelim vars (p::(map cfn oeqs),neqs) sgns' else if neqs = [] then True else let q = end_itlist (poly_mul vars) neqs in poly_nondiv vars sgns' p (poly_pow vars q (degree vars p)));; (* ------------------------------------------------------------------------- *) (* Basic complex quantifier elimination on actual existential formula. *) (* ------------------------------------------------------------------------- *) let init_sgns = [Fn("1",[]),Positive; Fn("0",[]),Zero];; let basic_complex_qelim vars (Exists(x,p)) = let eqs,neqs = partition (non negative) (conjuncts p) in cqelim (x::vars) (map lhs eqs,map (lhs ** negate) neqs) init_sgns;; (* ------------------------------------------------------------------------- *) (* Full quantifier elimination. *) (* ------------------------------------------------------------------------- *) let complex_qelim = simplify ** evalc ** lift_qelim polyatom (dnf ** cnnf (fun x -> x) ** evalc) basic_complex_qelim;; (* ------------------------------------------------------------------------- *) (* Examples. *) (* ------------------------------------------------------------------------- *) START_INTERACTIVE;; complex_qelim < x^4 + 1 = 0>>;; complex_qelim < x^4 + c = 0>>;; complex_qelim < x^4 + c = 0) <=> c = 1>>;; complex_qelim < a * x * y = c /\ a * (x + y) + b = 0>>;; (* ------------------------------------------------------------------------- *) (* More tests, not in the main text. *) (* ------------------------------------------------------------------------- *) let polytest tm = time (polynate (fvt tm)) tm;; let lagrange_4 = polytest <<|(((x1^2) + (x2^2) + (x3^2) + (x4^2)) * ((y1^2) + (y2^2) + (y3^2) + (y4^2))) - ((((((x1*y1) - (x2*y2)) - (x3*y3)) - (x4*y4))^2) + (((((x1*y2) + (x2*y1)) + (x3*y4)) - (x4*y3))^2) + (((((x1*y3) - (x2*y4)) + (x3*y1)) + (x4*y2))^2) + (((((x1*y4) + (x2*y3)) - (x3*y2)) + (x4*y1))^2))|>>;; let lagrange_8 = polytest <<|((p1^2 + q1^2 + r1^2 + s1^2 + t1^2 + u1^2 + v1^2 + w1^2) * (p2^2 + q2^2 + r2^2 + s2^2 + t2^2 + u2^2 + v2^2 + w2^2)) - ((p1 * p2 - q1 * q2 - r1 * r2 - s1 * s2 - t1 * t2 - u1 * u2 - v1 * v2 - w1* w2)^2 + (p1 * q2 + q1 * p2 + r1 * s2 - s1 * r2 + t1 * u2 - u1 * t2 - v1 * w2 + w1* v2)^2 + (p1 * r2 - q1 * s2 + r1 * p2 + s1 * q2 + t1 * v2 + u1 * w2 - v1 * t2 - w1* u2)^2 + (p1 * s2 + q1 * r2 - r1 * q2 + s1 * p2 + t1 * w2 - u1 * v2 + v1 * u2 - w1* t2)^2 + (p1 * t2 - q1 * u2 - r1 * v2 - s1 * w2 + t1 * p2 + u1 * q2 + v1 * r2 + w1* s2)^2 + (p1 * u2 + q1 * t2 - r1 * w2 + s1 * v2 - t1 * q2 + u1 * p2 - v1 * s2 + w1* r2)^2 + (p1 * v2 + q1 * w2 + r1 * t2 - s1 * u2 - t1 * r2 + u1 * s2 + v1 * p2 - w1* q2)^2 + (p1 * w2 - q1 * v2 + r1 * u2 + s1 * t2 - t1 * s2 - u1 * r2 + v1 * q2 + w1* p2)^2)|>>;; let liouville = polytest <<|6 * (x1^2 + x2^2 + x3^2 + x4^2)^2 - (((x1 + x2)^4 + (x1 + x3)^4 + (x1 + x4)^4 + (x2 + x3)^4 + (x2 + x4)^4 + (x3 + x4)^4) + ((x1 - x2)^4 + (x1 - x3)^4 + (x1 - x4)^4 + (x2 - x3)^4 + (x2 - x4)^4 + (x3 - x4)^4))|>>;; let fleck = polytest <<|60 * (x1^2 + x2^2 + x3^2 + x4^2)^3 - (((x1 + x2 + x3)^6 + (x1 + x2 - x3)^6 + (x1 - x2 + x3)^6 + (x1 - x2 - x3)^6 + (x1 + x2 + x4)^6 + (x1 + x2 - x4)^6 + (x1 - x2 + x4)^6 + (x1 - x2 - x4)^6 + (x1 + x3 + x4)^6 + (x1 + x3 - x4)^6 + (x1 - x3 + x4)^6 + (x1 - x3 - x4)^6 + (x2 + x3 + x4)^6 + (x2 + x3 - x4)^6 + (x2 - x3 + x4)^6 + (x2 - x3 - x4)^6) + 2 * ((x1 + x2)^6 + (x1 - x2)^6 + (x1 + x3)^6 + (x1 - x3)^6 + (x1 + x4)^6 + (x1 - x4)^6 + (x2 + x3)^6 + (x2 - x3)^6 + (x2 + x4)^6 + (x2 - x4)^6 + (x3 + x4)^6 + (x3 - x4)^6) + 36 * (x1^6 + x2^6 + x3^6 + x4^6))|>>;; let hurwitz = polytest <<|5040 * (x1^2 + x2^2 + x3^2 + x4^2)^4 - (6 * ((x1 + x2 + x3 + x4)^8 + (x1 + x2 + x3 - x4)^8 + (x1 + x2 - x3 + x4)^8 + (x1 + x2 - x3 - x4)^8 + (x1 - x2 + x3 + x4)^8 + (x1 - x2 + x3 - x4)^8 + (x1 - x2 - x3 + x4)^8 + (x1 - x2 - x3 - x4)^8) + ((2 * x1 + x2 + x3)^8 + (2 * x1 + x2 - x3)^8 + (2 * x1 - x2 + x3)^8 + (2 * x1 - x2 - x3)^8 + (2 * x1 + x2 + x4)^8 + (2 * x1 + x2 - x4)^8 + (2 * x1 - x2 + x4)^8 + (2 * x1 - x2 - x4)^8 + (2 * x1 + x3 + x4)^8 + (2 * x1 + x3 - x4)^8 + (2 * x1 - x3 + x4)^8 + (2 * x1 - x3 - x4)^8 + (2 * x2 + x3 + x4)^8 + (2 * x2 + x3 - x4)^8 + (2 * x2 - x3 + x4)^8 + (2 * x2 - x3 - x4)^8 + (x1 + 2 * x2 + x3)^8 + (x1 + 2 * x2 - x3)^8 + (x1 - 2 * x2 + x3)^8 + (x1 - 2 * x2 - x3)^8 + (x1 + 2 * x2 + x4)^8 + (x1 + 2 * x2 - x4)^8 + (x1 - 2 * x2 + x4)^8 + (x1 - 2 * x2 - x4)^8 + (x1 + 2 * x3 + x4)^8 + (x1 + 2 * x3 - x4)^8 + (x1 - 2 * x3 + x4)^8 + (x1 - 2 * x3 - x4)^8 + (x2 + 2 * x3 + x4)^8 + (x2 + 2 * x3 - x4)^8 + (x2 - 2 * x3 + x4)^8 + (x2 - 2 * x3 - x4)^8 + (x1 + x2 + 2 * x3)^8 + (x1 + x2 - 2 * x3)^8 + (x1 - x2 + 2 * x3)^8 + (x1 - x2 - 2 * x3)^8 + (x1 + x2 + 2 * x4)^8 + (x1 + x2 - 2 * x4)^8 + (x1 - x2 + 2 * x4)^8 + (x1 - x2 - 2 * x4)^8 + (x1 + x3 + 2 * x4)^8 + (x1 + x3 - 2 * x4)^8 + (x1 - x3 + 2 * x4)^8 + (x1 - x3 - 2 * x4)^8 + (x2 + x3 + 2 * x4)^8 + (x2 + x3 - 2 * x4)^8 + (x2 - x3 + 2 * x4)^8 + (x2 - x3 - 2 * x4)^8) + 60 * ((x1 + x2)^8 + (x1 - x2)^8 + (x1 + x3)^8 + (x1 - x3)^8 + (x1 + x4)^8 + (x1 - x4)^8 + (x2 + x3)^8 + (x2 - x3)^8 + (x2 + x4)^8 + (x2 - x4)^8 + (x3 + x4)^8 + (x3 - x4)^8) + 6 * ((2 * x1)^8 + (2 * x2)^8 + (2 * x3)^8 + (2 * x4)^8))|>>;; let schur = polytest <<|22680 * (x1^2 + x2^2 + x3^2 + x4^2)^5 - (9 * ((2 * x1)^10 + (2 * x2)^10 + (2 * x3)^10 + (2 * x4)^10) + 180 * ((x1 + x2)^10 + (x1 - x2)^10 + (x1 + x3)^10 + (x1 - x3)^10 + (x1 + x4)^10 + (x1 - x4)^10 + (x2 + x3)^10 + (x2 - x3)^10 + (x2 + x4)^10 + (x2 - x4)^10 + (x3 + x4)^10 + (x3 - x4)^10) + ((2 * x1 + x2 + x3)^10 + (2 * x1 + x2 - x3)^10 + (2 * x1 - x2 + x3)^10 + (2 * x1 - x2 - x3)^10 + (2 * x1 + x2 + x4)^10 + (2 * x1 + x2 - x4)^10 + (2 * x1 - x2 + x4)^10 + (2 * x1 - x2 - x4)^10 + (2 * x1 + x3 + x4)^10 + (2 * x1 + x3 - x4)^10 + (2 * x1 - x3 + x4)^10 + (2 * x1 - x3 - x4)^10 + (2 * x2 + x3 + x4)^10 + (2 * x2 + x3 - x4)^10 + (2 * x2 - x3 + x4)^10 + (2 * x2 - x3 - x4)^10 + (x1 + 2 * x2 + x3)^10 + (x1 + 2 * x2 - x3)^10 + (x1 - 2 * x2 + x3)^10 + (x1 - 2 * x2 - x3)^10 + (x1 + 2 * x2 + x4)^10 + (x1 + 2 * x2 - x4)^10 + (x1 - 2 * x2 + x4)^10 + (x1 - 2 * x2 - x4)^10 + (x1 + 2 * x3 + x4)^10 + (x1 + 2 * x3 - x4)^10 + (x1 - 2 * x3 + x4)^10 + (x1 - 2 * x3 - x4)^10 + (x2 + 2 * x3 + x4)^10 + (x2 + 2 * x3 - x4)^10 + (x2 - 2 * x3 + x4)^10 + (x2 - 2 * x3 - x4)^10 + (x1 + x2 + 2 * x3)^10 + (x1 + x2 - 2 * x3)^10 + (x1 - x2 + 2 * x3)^10 + (x1 - x2 - 2 * x3)^10 + (x1 + x2 + 2 * x4)^10 + (x1 + x2 - 2 * x4)^10 + (x1 - x2 + 2 * x4)^10 + (x1 - x2 - 2 * x4)^10 + (x1 + x3 + 2 * x4)^10 + (x1 + x3 - 2 * x4)^10 + (x1 - x3 + 2 * x4)^10 + (x1 - x3 - 2 * x4)^10 + (x2 + x3 + 2 * x4)^10 + (x2 + x3 - 2 * x4)^10 + (x2 - x3 + 2 * x4)^10 + (x2 - x3 - 2 * x4)^10) + 9 * ((x1 + x2 + x3 + x4)^10 + (x1 + x2 + x3 - x4)^10 + (x1 + x2 - x3 + x4)^10 + (x1 + x2 - x3 - x4)^10 + (x1 - x2 + x3 + x4)^10 + (x1 - x2 + x3 - x4)^10 + (x1 - x2 - x3 + x4)^10 + (x1 - x2 - x3 - x4)^10))|>>;; let complex_qelim_all = time complex_qelim ** generalize;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim < x^4 + 1 = 0>>;; time complex_qelim < x^4 + 2 = 0>>;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim < x^4 + 2 = 0>>;; time complex_qelim < x^4 + 2 = 0>>;; time complex_qelim <>;; (*** This works by a combination with grobner_decide but seems slow like this: complex_qelim < z = x) /\ x = y ==> a * x * y = c /\ a * (x + y) + b = 0>>;; *** and w/o the condition, it's false I think complex_qelim < z = x \/ z = y) ==> a * x * y = c /\ a * (x + y) + b = 0>>;; *** because the following is! ***) complex_qelim < z = x) ==> a * x * x = c /\ a * (x + x) + b = 0>>;; (*** In fact we have this, tho' I don't know if that's interesting ***) complex_qelim < (a * x * y = c) /\ (a * (x + y) + b = 0)) <=> ~(x = y)>>;; time complex_qelim < (y_1^2 = y_2^2)>>;; time complex_qelim < (x * y)^2 = 6>>;; time complex_qelim < (x^4 + 1 = 0)>>;; time complex_qelim < (x^4 + 1 = 0)>>;; time complex_qelim <<~(exists a x y. (a^2 = 2) /\ (x^2 + a * x + 1 = 0) /\ (y * (x^4 + 1) + 1 = 0))>>;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim < exists x y. (y * x^2 = a) /\ (y * x^2 + x = b)>>;; time complex_qelim < (a * x * y = c) /\ (a * (x + y) + b = 0)>>;; time complex_qelim <<~(forall a b c x y. (a * x^2 + b * x + c = 0) /\ (a * y^2 + b * y + c = 0) ==> (a * x * y = c) /\ (a * (x + y) + b = 0))>>;; time complex_qelim < (y_1^2 = y_2^2)>>;; time complex_qelim < exists x y. (a1 * x + b1 * y = c1) /\ (a2 * x + b2 * y = c2)>>;; (* ------------------------------------------------------------------------- *) (* This seems harder, so see how many quantifiers are feasible. *) (* ------------------------------------------------------------------------- *) time complex_qelim <<(a * x^2 + b * x + c = 0) /\ (a * y^2 + b * y + c = 0) /\ (forall z. (a * z^2 + b * z + c = 0) ==> (z = x) \/ (z = y)) ==> (a * x * y = c) /\ (a * (x + y) + b = 0)>>;; time complex_qelim < (z = x) \/ (z = y)) ==> (a * x * y = c) /\ (a * (x + y) + b = 0)>>;; (**** feasible but lengthy? time complex_qelim < (z = x) \/ (z = y)) ==> (a * x * y = c) /\ (a * (x + y) + b = 0)>>;; time complex_qelim < (z = x) \/ (z = y)) ==> (a * x * y = c) /\ (a * (x + y) + b = 0)>>;; ************) (********* This seems too hard time complex_qelim < (z = x) \/ (z = y)) ==> (a * x * y = c) /\ (a * (x + y) + b = 0)>>;; **************) time complex_qelim <<~(forall x1 y1 x2 y2 x3 y3. exists x0 y0. (x1 - x0)^2 + (y1 - y0)^2 = (x2 - x0)^2 + (y2 - y0)^2 /\ (x2 - x0)^2 + (y2 - y0)^2 = (x3 - x0)^2 + (y3 - y0)^2)>>;; time complex_qelim < (a = 0) /\ (b = 0) /\ (c = 0) \/ ~(a = 0) /\ ~(b^2 = 4 * a * c)>>;; time complex_qelim <<~(forall x1 y1 x2 y2 x3 y3 x0 y0 x0' y0'. (x1 - x0)^2 + (y1 - y0)^2 = (x2 - x0)^2 + (y2 - y0)^2 /\ (x2 - x0)^2 + (y2 - y0)^2 = (x3 - x0)^2 + (y3 - y0)^2 /\ (x1 - x0')^2 + (y1 - y0')^2 = (x2 - x0')^2 + (y2 - y0')^2 /\ (x2 - x0')^2 + (y2 - y0')^2 = (x3 - x0')^2 + (y3 - y0')^2 ==> x0 = x0' /\ y0 = y0')>>;; time complex_qelim < a * (x + y) + b = 0>>;; time complex_qelim < (a * (x + y) + b = 0)>>;; complex_qelim_all <<~(y_1 = 2 * y_3 /\ y_2 = 2 * y_4 /\ y_1 * y_3 = y_2 * y_4 /\ (y_1^2 - y_2^2) * z = 1)>>;; time complex_qelim < (y_1^2 = y_2^2)>>;; (************ complex_qelim_all <<~((c^2 = a^2 + b^2) /\ (c^2 = x0^2 + (y0 - b)^2) /\ (y0 * t1 = a + x0) /\ (y0 * t2 = a - x0) /\ ((1 - t1 * t2) * t = t1 + t2) /\ (u * (b * t - a) = 1) /\ (v1 * a + v2 * x0 + v3 * y0 = 1))>>;; complex_qelim_all <<(c^2 = a^2 + b^2) /\ (c^2 = x0^2 + (y0 - b)^2) /\ (y0 * t1 = a + x0) /\ (y0 * t2 = a - x0) /\ ((1 - t1 * t2) * t = t1 + t2) /\ (~(a = 0) \/ ~(x0 = 0) \/ ~(y0 = 0)) ==> (b * t = a)>>;; *********) complex_qelim_all <<(x1 = u3) /\ (x1 * (u2 - u1) = x2 * u3) /\ (x4 * (x2 - u1) = x1 * (x3 - u1)) /\ (x3 * u3 = x4 * u2) /\ ~(u1 = 0) /\ ~(u3 = 0) ==> (x3^2 + x4^2 = (u2 - x3)^2 + (u3 - x4)^2)>>;; complex_qelim_all <<(u1 * x1 - u1 * u3 = 0) /\ (u3 * x2 - (u2 - u1) * x1 = 0) /\ (x1 * x4 - (x2 - u1) * x3 - u1 * x1 = 0) /\ (u3 * x4 - u2 * x3 = 0) /\ ~(u1 = 0) /\ ~(u3 = 0) ==> (2 * u2 * x4 + 2 * u3 * x3 - u3^2 - u2^2 = 0)>>;; complex_qelim_all <<(y1 * y3 + x1 * x3 = 0) /\ (y3 * (y2 - y3) + (x2 - x3) * x3 = 0) /\ ~(x3 = 0) /\ ~(y3 = 0) ==> (y1 * (x2 - x3) = x1 * (y2 - y3))>>;; (********** complex_qelim_all <<(2 * u2 * x2 + 2 * u3 * x1 - u3^2 - u2^2 = 0) /\ (2 * u1 * x2 - u1^2 = 0) /\ (-(x3^2) + 2 * x2 * x3 + 2 * u4 * x1 - u4^2 = 0) /\ (u3 * x5 + (-(u2) + u1) * x4 - u1 * u3 = 0) /\ ((u2 - u1) * x5 + u3 * x4 + (-(u2) + u1) * x3 - u3 * u4 = 0) /\ (u3 * x7 - u2 * x6 = 0) /\ (u2 * x7 + u3 * x6 - u2 * x3 - u3 * u4 = 0) /\ ~(4 * u1 * u3 = 0) /\ ~(2 * u1 = 0) /\ ~(-(u3^2) - u2^2 + 2 * u1 * u2 - u1^2 = 0) /\ ~(u3 = 0) /\ ~(-(u3^2) - u2^2 = 0) /\ ~(u2 = 0) ==> (x4 * x7 + (-(x5) + x3) * x6 - x3 * x4 = 0)>>;; time complex_qelim <>;; time complex_qelim <>;; *********) time complex_qelim < a * x * y = c /\ a * (x + y) + b = 0>>;; complex_qelim_all < a * x * y = c /\ a * (x + y) + b = 0>>;; (* ------------------------------------------------------------------------- *) (* The Colmerauer example. *) (* ------------------------------------------------------------------------- *) (********* This works, but is quite slow. And it's false! Presumably we actually need to use ordering properties associated with absolute values let colmerauer = complex_qelim_all <<(x_1 + x_3 = (x_2) \/ x_1 + x_3 = -(x_2)) /\ (x_2 + x_4 = (x_3) \/ x_2 + x_4 = -(x_3)) /\ (x_3 + x_5 = (x_4) \/ x_3 + x_5 = -(x_4)) /\ (x_4 + x_6 = (x_5) \/ x_4 + x_6 = -(x_5)) /\ (x_5 + x_7 = (x_6) \/ x_5 + x_7 = -(x_6)) /\ (x_6 + x_8 = (x_7) \/ x_6 + x_8 = -(x_7)) /\ (x_7 + x_9 = (x_8) \/ x_7 + x_9 = -(x_8)) /\ (x_8 + x_10 = (x_9) \/ x_8 + x_10 = -(x_9)) /\ (x_9 + x_11 = (x_10) \/ x_9 + x_11 = -(x_10)) ==> x_1 = x_10 /\ x_2 = x_11>>;; ***********) (* ------------------------------------------------------------------------- *) (* Checking resultants from Maple. *) (* ------------------------------------------------------------------------- *) (*** interface(prettyprint=0); resultant(a * x^2 + b * x + c, 2 * a * x + b,x); ***) time complex_qelim < (4*a^2*c-b^2*a = 0)>>;; time complex_qelim < d^2*c-e*d*b+a*e^2 = 0>>;; time complex_qelim < d^2*c^2-2*d*c*a*f+a^2*f^2-e*d*b*c-e*b*a*f+a*e^2*c+f*d*b^2 = 0>>;; (**** No hope for this one I think time complex_qelim < e^3*d^2+3*e*d*g*a*f-2*e^2*d*g*b-g^2*a*f*b+g^2*e*b^2-f*e^2*c*d+f^2*c*g*a-f*e*c* g*b+f^2*e*b*d-f^3*a*d+g*e^2*c^2-2*e*c*a*g^2+a^2*g^3 = 0>>;; ****) (* ------------------------------------------------------------------------- *) (* Some trigonometric addition formulas (checking stuff from Maple). *) (* ------------------------------------------------------------------------- *) time complex_qelim < (2 * y^2 - 1)^2 + (2 * x * y)^2 = 1>>;; (* ------------------------------------------------------------------------- *) (* The examples from my thesis. *) (* ------------------------------------------------------------------------- *) time complex_qelim < 2 * s - (2 * s * c * c - s^3) = 3 * s^3>>;; time complex_qelim <>;; time complex_qelim < (((9 * u^8) * v) * v - (u * u^9)) * 128 + (((7 * u^6) * v) * v - (u * u^7)) * 144 + (((5 * u^4) * v) * v - (u * u^5)) * 168 + (((3 * u^2) * v) * v - (u * u^3)) * 210 + (v * v - (u * u)) * 315 + 1280 * u^10 = 315>>;; (* ------------------------------------------------------------------------- *) (* Deliberately silly examples from Poizat's model theory book (6.6). *) (* ------------------------------------------------------------------------- *) time complex_qelim <>;; time complex_qelim <>;; (* ------------------------------------------------------------------------- *) (* Actually prove simple equivalences. *) (* ------------------------------------------------------------------------- *) time complex_qelim < ~(x = 0) \/ ~(y = 0)>>;; time complex_qelim < ~(x = 0) \/ ~(y = 0) \/ z = 0>>;; (* ------------------------------------------------------------------------- *) (* Invertibility of 2x2 matrix in terms of nonzero determinant. *) (* ------------------------------------------------------------------------- *) time complex_qelim <>;; time complex_qelim < ~(a * d = b * c)>>;; (* ------------------------------------------------------------------------- *) (* Inspired by Cardano's formula for a cubic. Not all complex cbrts work. *) (* ------------------------------------------------------------------------- *) time complex_qelim < x^3 + m * x = n>>;; time complex_qelim < exists ct cu. ct^3 = t /\ cu^3 = u /\ (x = ct - cu ==> x^3 + m * x = n)>>;; (* ------------------------------------------------------------------------- *) (* SOS in rational functions for Motzkin polynomial (dehomogenized). *) (* Of course these are just trivial normalization, nothing deep. *) (* ------------------------------------------------------------------------- *) time complex_qelim <>;; time complex_qelim <>;; time complex_qelim <>;; time complex_qelim <>;; (* ------------------------------------------------------------------------- *) (* A cute bilinear identity -- see ch14 of Rajwade's "Squares" for more. *) (* ------------------------------------------------------------------------- *) polytest <<|(x_1^2 + x_2^2 + x_3^2 + x_4^2 + x_5^2 + x_6^2 + x_7^2 + x_8^2 + x_9^2) * (y_1^2 + y_2^2 + y_3^2 + y_4^2 + y_5^2 + y_6^2 + y_7^2 + y_8^2 + y_9^2 + y_10^2 + y_11^2 + y_12^2 + y_13^2 + y_14^2 + y_15^2 + y_16^2) - ((0 + x_1 * y_1 + x_2 * y_2 + x_3 * y_3 + x_4 * y_4 + x_5 * y_5 + x_6 * y_6 + x_7 * y_7 + x_8 * y_8 + x_9 * y_9)^2 + (0 - x_2 * y_1 + x_1 * y_2 + x_4 * y_3 - x_3 * y_4 + x_6 * y_5 - x_5 * y_6 - x_8 * y_7 + x_7 * y_8 + x_9 * y_10)^2 + (0 - x_3 * y_1 - x_4 * y_2 + x_1 * y_3 + x_2 * y_4 + x_7 * y_5 + x_8 * y_6 - x_5 * y_7 - x_6 * y_8 + x_9 * y_11)^2 + (0 - x_4 * y_1 + x_3 * y_2 - x_2 * y_3 + x_1 * y_4 + x_8 * y_5 - x_7 * y_6 + x_6 * y_7 - x_5 * y_8 + x_9 * y_12)^2 + (0 - x_5 * y_1 - x_6 * y_2 - x_7 * y_3 - x_8 * y_4 + x_1 * y_5 + x_2 * y_6 + x_3 * y_7 + x_4 * y_8 + x_9 * y_13)^2 + (0 - x_6 * y_1 + x_5 * y_2 - x_8 * y_3 + x_7 * y_4 - x_2 * y_5 + x_1 * y_6 - x_4 * y_7 + x_3 * y_8 + x_9 * y_14)^2 + (0 - x_7 * y_1 + x_8 * y_2 + x_5 * y_3 - x_6 * y_4 - x_3 * y_5 + x_4 * y_6 + x_1 * y_7 - x_2 * y_8 + x_9 * y_15)^2 + (0 - x_8 * y_1 - x_7 * y_2 + x_6 * y_3 + x_5 * y_4 - x_4 * y_5 - x_3 * y_6 + x_2 * y_7 + x_1 * y_8 + x_9 * y_16)^2 + (0 - x_9 * y_1 + x_1 * y_9 - x_2 * y_10 - x_3 * y_11 - x_4 * y_12 - x_5 * y_13 - x_6 * y_14 - x_7 * y_15 - x_8 * y_16)^2 + (0 - x_9 * y_2 + x_2 * y_9 + x_1 * y_10 - x_4 * y_11 + x_3 * y_12 - x_6 * y_13 + x_5 * y_14 + x_8 * y_15 - x_7 * y_16)^2 + (0 - x_9 * y_3 + x_3 * y_9 + x_4 * y_10 + x_1 * y_11 - x_2 * y_12 - x_7 * y_13 - x_8 * y_14 + x_5 * y_15 + x_6 * y_16)^2 + (0 - x_9 * y_4 + x_4 * y_9 - x_3 * y_10 + x_2 * y_11 + x_1 * y_12 - x_8 * y_13 + x_7 * y_14 - x_6 * y_15 + x_5 * y_16)^2 + (0 - x_9 * y_5 + x_5 * y_9 + x_6 * y_10 + x_7 * y_11 + x_8 * y_12 + x_1 * y_13 - x_2 * y_14 - x_3 * y_15 - x_4 * y_16)^2 + (0 - x_9 * y_6 + x_6 * y_9 - x_5 * y_10 + x_8 * y_11 - x_7 * y_12 + x_2 * y_13 + x_1 * y_14 + x_4 * y_15 - x_3 * y_16)^2 + (0 - x_9 * y_7 + x_7 * y_9 - x_8 * y_10 - x_5 * y_11 + x_6 * y_12 + x_3 * y_13 - x_4 * y_14 + x_1 * y_15 + x_2 * y_16)^2 + (0 - x_9 * y_8 + x_8 * y_9 + x_7 * y_10 - x_6 * y_11 - x_5 * y_12 + x_4 * y_13 + x_3 * y_14 - x_2 * y_15 + x_1 * y_16)^2)|>>;; (* ------------------------------------------------------------------------- *) (* This is essentially the Cauchy-Riemann conditions for a differential. *) (* ------------------------------------------------------------------------- *) time complex_qelim < (a = d)>>;; time complex_qelim < (a = d) /\ (b = -(c))>>;; time complex_qelim < (a = d) /\ (b = -(c))>>;; (* ------------------------------------------------------------------------- *) (* Finding non-trivial perpendiculars to lines. *) (* ------------------------------------------------------------------------- *) complex_qelim <>;; END_INTERACTIVE;;