From: Cezary Kaliszyk Date: Tue, 27 Aug 2013 11:24:06 +0000 (+0200) Subject: Update from HH X-Git-Url: http://colo12-c703.uibk.ac.at/git/?a=commitdiff_plain;p=Complex%20Analysis%2F.git Update from HH --- cb6a49e45a2591f628639c525c88e9f844be134b diff --git a/Complex/complex_grobner.ml b/Complex/complex_grobner.ml new file mode 100644 index 0000000..30d8b17 --- /dev/null +++ b/Complex/complex_grobner.ml @@ -0,0 +1,501 @@ +(* ========================================================================= *) +(* Grobner basis algorithm. *) +(* ========================================================================= *) + +needs "Complex/complexnumbers.ml";; +needs "Complex/quelim.ml";; + +prioritize_complex();; + +(* ------------------------------------------------------------------------- *) +(* Utility functions. *) +(* ------------------------------------------------------------------------- *) + +let allpairs f l1 l2 = + itlist ((@) o C map l2 o f) l1 [];; + +let rec merge ord l1 l2 = + match l1 with + [] -> l2 + | h1::t1 -> match l2 with + [] -> l1 + | h2::t2 -> if ord h1 h2 then h1::(merge ord t1 l2) + else h2::(merge ord l1 t2);; + +let sort ord = + let rec mergepairs l1 l2 = + match (l1,l2) with + ([s],[]) -> s + | (l,[]) -> mergepairs [] l + | (l,[s1]) -> mergepairs (s1::l) [] + | (l,(s1::s2::ss)) -> mergepairs ((merge ord s1 s2)::l) ss in + fun l -> if l = [] then [] else mergepairs [] (map (fun x -> [x]) l);; + +(* ------------------------------------------------------------------------- *) +(* Type for recording history, i.e. how a polynomial was obtained. *) +(* ------------------------------------------------------------------------- *) + +type history = + Start of int + | Mmul of (num * (int list)) * history + | Add of history * history;; + +(* ------------------------------------------------------------------------- *) +(* Conversion of leaves, i.e. variables and constant rational constants. *) +(* ------------------------------------------------------------------------- *) + +let grob_var vars tm = + let res = map (fun i -> if i = tm then 1 else 0) vars in + if exists (fun x -> x <> 0) res then [Int 1,res] else failwith "grob_var";; + +let grob_const = + let cx_tm = `Cx` in + fun vars tm -> + try let l,r = dest_comb tm in + if l = cx_tm then + let x = rat_of_term r in + if x =/ Int 0 then [] else [x,map (fun v -> 0) vars] + else failwith "" + with Failure _ -> failwith "grob_const";; + +(* ------------------------------------------------------------------------- *) +(* Monomial ordering. *) +(* ------------------------------------------------------------------------- *) + +let morder_lt = + let rec lexorder l1 l2 = + match (l1,l2) with + [],[] -> false + | (x1::o1,x2::o2) -> x1 > x2 or x1 = x2 & lexorder o1 o2 + | _ -> failwith "morder: inconsistent monomial lengths" in + fun m1 m2 -> let n1 = itlist (+) m1 0 + and n2 = itlist (+) m2 0 in + n1 < n2 or n1 = n2 & lexorder m1 m2;; + +let morder_le m1 m2 = morder_lt m1 m2 or (m1 = m2);; + +let morder_gt m1 m2 = morder_lt m2 m1;; + +(* ------------------------------------------------------------------------- *) +(* Arithmetic on canonical polynomials. *) +(* ------------------------------------------------------------------------- *) + +let grob_neg = map (fun (c,m) -> (minus_num c,m));; + +let rec grob_add l1 l2 = + match (l1,l2) with + ([],l2) -> l2 + | (l1,[]) -> l1 + | ((c1,m1)::o1,(c2,m2)::o2) -> + if m1 = m2 then + let c = c1+/c2 and rest = grob_add o1 o2 in + if c =/ Int 0 then rest else (c,m1)::rest + else if morder_lt m2 m1 then (c1,m1)::(grob_add o1 l2) + else (c2,m2)::(grob_add l1 o2);; + +let grob_sub l1 l2 = grob_add l1 (grob_neg l2);; + +let grob_mmul (c1,m1) (c2,m2) = (c1*/c2,map2 (+) m1 m2);; + +let rec grob_cmul cm pol = map (grob_mmul cm) pol;; + +let rec grob_mul l1 l2 = + match l1 with + [] -> [] + | (h1::t1) -> grob_add (grob_cmul h1 l2) (grob_mul t1 l2);; + +let rec grob_pow vars l n = + if n < 0 then failwith "grob_pow: negative power" + else if n = 0 then [Int 1,map (fun v -> 0) vars] + else grob_mul l (grob_pow vars l (n - 1));; + +(* ------------------------------------------------------------------------- *) +(* Monomial division operation. *) +(* ------------------------------------------------------------------------- *) + +let mdiv (c1,m1) (c2,m2) = + (c1//c2, + map2 (fun n1 n2 -> if n1 < n2 then failwith "mdiv" else n1-n2) m1 m2);; + +(* ------------------------------------------------------------------------- *) +(* Lowest common multiple of two monomials. *) +(* ------------------------------------------------------------------------- *) + +let mlcm (c1,m1) (c2,m2) = (Int 1,map2 max m1 m2);; + +(* ------------------------------------------------------------------------- *) +(* Reduce monomial cm by polynomial pol, returning replacement for cm. *) +(* ------------------------------------------------------------------------- *) + +let reduce1 cm (pol,hpol) = + match pol with + [] -> failwith "reduce1" + | cm1::cms -> try let (c,m) = mdiv cm cm1 in + (grob_cmul (minus_num c,m) cms, + Mmul((minus_num c,m),hpol)) + with Failure _ -> failwith "reduce1";; + +(* ------------------------------------------------------------------------- *) +(* Try this for all polynomials in a basis. *) +(* ------------------------------------------------------------------------- *) + +let reduceb cm basis = tryfind (fun p -> reduce1 cm p) basis;; + +(* ------------------------------------------------------------------------- *) +(* Reduction of a polynomial (always picking largest monomial possible). *) +(* ------------------------------------------------------------------------- *) + +let rec reduce basis (pol,hist) = + match pol with + [] -> (pol,hist) + | cm::ptl -> try let q,hnew = reduceb cm basis in + reduce basis (grob_add q ptl,Add(hnew,hist)) + with Failure _ -> + let q,hist' = reduce basis (ptl,hist) in + cm::q,hist';; + +(* ------------------------------------------------------------------------- *) +(* Check for orthogonality w.r.t. LCM. *) +(* ------------------------------------------------------------------------- *) + +let orthogonal l p1 p2 = + snd l = snd(grob_mmul (hd p1) (hd p2));; + +(* ------------------------------------------------------------------------- *) +(* Compute S-polynomial of two polynomials. *) +(* ------------------------------------------------------------------------- *) + +let spoly cm ph1 ph2 = + match (ph1,ph2) with + ([],h),p -> ([],h) + | p,([],h) -> ([],h) + | (cm1::ptl1,his1),(cm2::ptl2,his2) -> + (grob_sub (grob_cmul (mdiv cm cm1) ptl1) + (grob_cmul (mdiv cm cm2) ptl2), + Add(Mmul(mdiv cm cm1,his1), + Mmul(mdiv (minus_num(fst cm),snd cm) cm2,his2)));; + +(* ------------------------------------------------------------------------- *) +(* Make a polynomial monic. *) +(* ------------------------------------------------------------------------- *) + +let monic (pol,hist) = + if pol = [] then (pol,hist) else + let c',m' = hd pol in + (map (fun (c,m) -> (c//c',m)) pol, + Mmul((Int 1 // c',map (K 0) m'),hist));; + +(* ------------------------------------------------------------------------- *) +(* The most popular heuristic is to order critical pairs by LCM monomial. *) +(* ------------------------------------------------------------------------- *) + +let forder ((c1,m1),_) ((c2,m2),_) = morder_lt m1 m2;; + +(* ------------------------------------------------------------------------- *) +(* Stupid stuff forced on us by lack of equality test on num type. *) +(* ------------------------------------------------------------------------- *) + +let rec poly_lt p q = + match (p,q) with + p,[] -> false + | [],q -> true + | (c1,m1)::o1,(c2,m2)::o2 -> + c1 c1 =/ c2 & m1 = m2) p1 p2;; + +let memx ((p1,h1),(p2,h2)) ppairs = + not (exists (fun ((q1,_),(q2,_)) -> poly_eq p1 q1 & poly_eq p2 q2) ppairs);; + +(* ------------------------------------------------------------------------- *) +(* Buchberger's second criterion. *) +(* ------------------------------------------------------------------------- *) + +let criterion2 basis (lcm,((p1,h1),(p2,h2))) opairs = + exists (fun g -> not(poly_eq (fst g) p1) & not(poly_eq (fst g) p2) & + can (mdiv lcm) (hd(fst g)) & + not(memx (align(g,(p1,h1))) (map snd opairs)) & + not(memx (align(g,(p2,h2))) (map snd opairs))) basis;; + +(* ------------------------------------------------------------------------- *) +(* Test for hitting constant polynomial. *) +(* ------------------------------------------------------------------------- *) + +let constant_poly p = + length p = 1 & forall ((=) 0) (snd(hd p));; + +(* ------------------------------------------------------------------------- *) +(* Grobner basis algorithm. *) +(* ------------------------------------------------------------------------- *) + +let rec grobner histories basis pairs = + print_string(string_of_int(length basis)^" basis elements and "^ + string_of_int(length pairs)^" critical pairs"); + print_newline(); + match pairs with + [] -> rev histories,basis + | (l,(p1,p2))::opairs -> + let (sp,hist) = monic (reduce basis (spoly l p1 p2)) in + if sp = [] or criterion2 basis (l,(p1,p2)) opairs + then grobner histories basis opairs else + let sph = sp,Start(length histories) in + if constant_poly sp + then grobner ((sp,hist)::histories) (sph::basis) [] else + let rawcps = + map (fun p -> mlcm (hd(fst p)) (hd sp),align(p,sph)) basis in + let newcps = filter + (fun (l,(p,q)) -> not(orthogonal l (fst p) (fst q))) rawcps in + grobner ((sp,hist)::histories) (sph::basis) + (merge forder opairs (sort forder newcps));; + +(* ------------------------------------------------------------------------- *) +(* Overall function. *) +(* ------------------------------------------------------------------------- *) + +let groebner pols = + let npols = map2 (fun p n -> p,Start n) pols (0--(length pols - 1)) in + let phists = filter (fun (p,_) -> p <> []) npols in + let bas0 = map monic phists in + let bas = map2 (fun (p,h) n -> p,Start n) bas0 + ((length bas0)--(2 * length bas0 - 1)) in + let prs0 = allpairs (fun x y -> x,y) bas bas in + let prs1 = filter (fun ((x,_),(y,_)) -> poly_lt x y) prs0 in + let prs2 = map (fun (p,q) -> mlcm (hd(fst p)) (hd(fst q)),(p,q)) prs1 in + let prs3 = filter (fun (l,(p,q)) -> not(orthogonal l (fst p) (fst q))) prs2 in + grobner (rev bas0 @ rev phists) bas (mergesort forder prs3);; + +(* ------------------------------------------------------------------------- *) +(* Alternative orthography. *) +(* ------------------------------------------------------------------------- *) + +let gr'o'bner = groebner;; + +(* ------------------------------------------------------------------------- *) +(* Conversion from HOL term. *) +(* ------------------------------------------------------------------------- *) + +let grobify_term = + let neg_tm = `(--):complex->complex` + and add_tm = `(+):complex->complex->complex` + and sub_tm = `(-):complex->complex->complex` + and mul_tm = `(*):complex->complex->complex` + and pow_tm = `(pow):complex->num->complex` in + let rec grobify_term vars tm = + try grob_var vars tm with Failure _ -> + try grob_const vars tm with Failure _ -> + let lop,r = dest_comb tm in + if lop = neg_tm then grob_neg(grobify_term vars r) else + let op,l = dest_comb lop in + if op = pow_tm then + grob_pow vars (grobify_term vars l) (dest_small_numeral r) + else + (if op = add_tm then grob_add else if op = sub_tm then grob_sub + else if op = mul_tm then grob_mul else failwith "unknown term") + (grobify_term vars l) (grobify_term vars r) in + fun vars tm -> + try grobify_term vars tm with Failure _ -> failwith "grobify_term";; + +let grobvars = + let neg_tm = `(--):complex->complex` + and add_tm = `(+):complex->complex->complex` + and sub_tm = `(-):complex->complex->complex` + and mul_tm = `(*):complex->complex->complex` + and pow_tm = `(pow):complex->num->complex` in + let rec grobvars tm acc = + if is_complex_const tm then acc + else if not (is_comb tm) then tm::acc else + let lop,r = dest_comb tm in + if lop = neg_tm then grobvars r acc + else if not (is_comb lop) then tm::acc else + let op,l = dest_comb lop in + if op = pow_tm then grobvars l acc + else if op = mul_tm or op = sub_tm or op = add_tm + then grobvars l (grobvars r acc) + else tm::acc in + fun tm -> setify(grobvars tm []);; + +let grobify_equations = + let zero_tm = `Cx(&0)` + and sub_tm = `(-):complex->complex->complex` + and complex_ty = `:complex` in + let grobify_equation vars tm = + let l,r = dest_eq tm in + if r <> zero_tm then grobify_term vars (mk_comb(mk_comb(sub_tm,l),r)) + else grobify_term vars l in + fun tm -> + let cjs = conjuncts tm in + let rawvars = itlist + (fun eq acc -> let l,r = dest_eq eq in + union (union (grobvars l) (grobvars r)) acc) cjs [] in + let vars = sort (fun x y -> x < y) rawvars in + vars,map(grobify_equation vars) cjs;; + +(* ------------------------------------------------------------------------- *) +(* Printer. *) +(* ------------------------------------------------------------------------- *) + +let string_of_monomial vars (c,m) = + let xns = filter (fun (x,y) -> y <> 0) (map2 (fun x y -> x,y) vars m) in + let xnstrs = map + (fun (x,n) -> x^(if n = 1 then "" else "^"^(string_of_int n))) xns in + if xns = [] then Num.string_of_num c else + let basstr = if c =/ Int 1 then "" else (Num.string_of_num c)^" * " in + basstr ^ end_itlist (fun s t -> s^" * "^t) xnstrs;; + +let string_of_polynomial vars l = + if l = [] then "0" else + end_itlist (fun s t -> s^" + "^t) (map (string_of_monomial vars) l);; + +(* ------------------------------------------------------------------------- *) +(* Resolve a proof. *) +(* ------------------------------------------------------------------------- *) + +let rec resolve_proof vars prf = + match prf with + Start n -> + [n,[Int 1,map (K 0) vars]] + | Mmul(pol,lin) -> + let lis = resolve_proof vars lin in + map (fun (n,p) -> n,grob_cmul pol p) lis + | Add(lin1,lin2) -> + let lis1 = resolve_proof vars lin1 + and lis2 = resolve_proof vars lin2 in + let dom = setify(union (map fst lis1) (map fst lis2)) in + map (fun n -> let a = try assoc n lis1 with Failure _ -> [] + and b = try assoc n lis2 with Failure _ -> [] in + n,grob_add a b) dom;; + +(* ------------------------------------------------------------------------- *) +(* Convert a polynomial back to HOL. *) +(* ------------------------------------------------------------------------- *) + +let holify_polynomial = + let complex_ty = `:complex` + and pow_tm = `(pow):complex->num->complex` + and mk_mul = mk_binop `(*):complex->complex->complex` + and mk_add = mk_binop `(+):complex->complex->complex` + and zero_tm = `Cx(&0)` + and add_tm = `(+):complex->complex->complex` + and mul_tm = `(*):complex->complex->complex` + and complex_term_of_rat = curry mk_comb `Cx` o term_of_rat in + let holify_varpow (v,n) = + if n = 1 then v + else list_mk_comb(pow_tm,[v;mk_small_numeral n]) in + let holify_monomial vars (c,m) = + let xps = map holify_varpow (filter (fun (_,n) -> n <> 0) (zip vars m)) in + end_itlist mk_mul (complex_term_of_rat c :: xps) in + let holify_polynomial vars p = + if p = [] then zero_tm + else end_itlist mk_add (map (holify_monomial vars) p) in + holify_polynomial;; + +(* ------------------------------------------------------------------------- *) +(* Recursively find the set of basis elements involved. *) +(* ------------------------------------------------------------------------- *) + +let dependencies = + let rec dependencies prf acc = + match prf with + Start n -> n::acc + | Mmul(pol,lin) -> dependencies lin acc + | Add(lin1,lin2) -> dependencies lin1 (dependencies lin2 acc) in + fun prf -> setify(dependencies prf []);; + +let rec involved deps sofar todo = + match todo with + [] -> sort (<) sofar + | (h::hs) -> + if mem h sofar then involved deps sofar hs + else involved deps (h::sofar) (el h deps @ hs);; + +(* ------------------------------------------------------------------------- *) +(* Refute a conjunction of equations in HOL. *) +(* ------------------------------------------------------------------------- *) + +let GROBNER_REFUTE = + let add_tm = `(+):complex->complex->complex` + and mul_tm = `(*):complex->complex->complex` in + let APPLY_pth = MATCH_MP(SIMPLE_COMPLEX_ARITH + `(x = y) ==> (x + Cx(--(&1)) * (y + Cx(&1)) = Cx(&0)) ==> F`) + and MK_ADD th1 th2 = MK_COMB(AP_TERM add_tm th1,th2) in + let rec holify_lincombs vars cjs prfs = + match prfs with + [] -> cjs + | (p::ps) -> + if p = [] then holify_lincombs vars (cjs @ [TRUTH]) ps else + let holis = map (fun (n,q) -> n,holify_polynomial vars q) p in + let ths = + map (fun (n,m) -> AP_TERM (mk_comb(mul_tm,m)) (el n cjs)) holis in + let th = CONV_RULE(BINOP_CONV COMPLEX_POLY_CONV) + (end_itlist MK_ADD ths) in + holify_lincombs vars (cjs @ [th]) ps in + fun tm -> + let vars,pols = grobify_equations tm in + let (prfs,gb) = groebner pols in + let (_,prf) = + find (fun (p,h) -> length p = 1 & forall ((=)0) (snd(hd p))) gb in + let deps = map (dependencies o snd) prfs + and depl = dependencies prf in + let need = involved deps [] depl in + let pprfs = + map2 (fun p n -> if mem n need then resolve_proof vars (snd p) else []) + prfs (0--(length prfs - 1)) + and ppr = resolve_proof vars prf in + let ths = CONJUNCTS(ASSUME tm) in + let th = last + (holify_lincombs vars ths (snd(chop_list(length ths) (pprfs @ [ppr])))) in + CONV_RULE COMPLEX_RAT_EQ_CONV th;; + +(* ------------------------------------------------------------------------- *) +(* Overall conversion. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_ARITH = + let pth0 = SIMPLE_COMPLEX_ARITH `(x = y) <=> (x - y = Cx(&0))` + and pth1 = prove + (`!x. ~(x = Cx(&0)) <=> ?z. z * x + Cx(&1) = Cx(&0)`, + CONV_TAC(time FULL_COMPLEX_QUELIM_CONV)) + and pth2a = prove + (`!x y u v. ~(x = y) \/ ~(u = v) <=> + ?w z. w * (x - y) + z * (u - v) + Cx(&1) = Cx(&0)`, + CONV_TAC(time FULL_COMPLEX_QUELIM_CONV)) + and pth2b = prove + (`!x y. ~(x = y) <=> ?z. z * (x - y) + Cx(&1) = Cx(&0)`, + CONV_TAC(time FULL_COMPLEX_QUELIM_CONV)) + and pth3 = TAUT `(p ==> F) ==> (~q <=> p) ==> q` in + let GEN_PRENEX_CONV = + GEN_REWRITE_CONV REDEPTH_CONV + [AND_FORALL_THM; + LEFT_AND_FORALL_THM; + RIGHT_AND_FORALL_THM; + LEFT_OR_FORALL_THM; + RIGHT_OR_FORALL_THM; + OR_EXISTS_THM; + LEFT_OR_EXISTS_THM; + RIGHT_OR_EXISTS_THM; + LEFT_AND_EXISTS_THM; + RIGHT_AND_EXISTS_THM] in + let INITIAL_CONV = + NNF_CONV THENC + GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth1] THENC + GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth2a] THENC + GEN_REWRITE_CONV ONCE_DEPTH_CONV [pth2b] THENC + ONCE_DEPTH_CONV(GEN_REWRITE_CONV I [pth0] o + check ((<>) `Cx(&0)` o rand)) THENC + GEN_PRENEX_CONV THENC + DNF_CONV in + fun tm -> + let avs = frees tm in + let tm' = list_mk_forall(avs,tm) in + let th1 = INITIAL_CONV(mk_neg tm') in + let evs,bod = strip_exists(rand(concl th1)) in + if is_forall bod then failwith "COMPLEX_ARITH: non-universal formula" else + let djs = disjuncts bod in + let th2 = end_itlist SIMPLE_DISJ_CASES(map GROBNER_REFUTE djs) in + let th3 = itlist SIMPLE_CHOOSE evs th2 in + SPECL avs (MATCH_MP (MATCH_MP pth3 (DISCH_ALL th3)) th1);; diff --git a/Complex/complex_real.ml b/Complex/complex_real.ml new file mode 100644 index 0000000..8098b25 --- /dev/null +++ b/Complex/complex_real.ml @@ -0,0 +1,9 @@ +(* ========================================================================= *) +(* Trivial restriction of complex Groebner bases to reals. *) +(* ========================================================================= *) + +let GROBNER_REAL_ARITH = + let trans_conv = GEN_REWRITE_CONV TOP_SWEEP_CONV + [GSYM CX_INJ; CX_POW; CX_MUL; CX_ADD; CX_NEG; CX_SUB] in + fun tm -> let th = trans_conv tm in + EQ_MP (SYM th) (COMPLEX_ARITH(rand(concl th)));; diff --git a/Complex/complex_transc.ml b/Complex/complex_transc.ml new file mode 100644 index 0000000..71dc894 --- /dev/null +++ b/Complex/complex_transc.ml @@ -0,0 +1,321 @@ +(* ========================================================================= *) +(* Complex transcendental functions. *) +(* ========================================================================= *) + +needs "Library/transc.ml";; +needs "Library/floor.ml";; +needs "Complex/complexnumbers.ml";; + +unparse_as_infix "exp";; +remove_interface "exp";; + +(* ------------------------------------------------------------------------- *) +(* Complex square roots. *) +(* ------------------------------------------------------------------------- *) + +let csqrt = new_definition + `csqrt(z) = if Im(z) = &0 then + if &0 <= Re(z) then complex(sqrt(Re(z)),&0) + else complex(&0,sqrt(--Re(z))) + else complex(sqrt((norm(z) + Re(z)) / &2), + (Im(z) / abs(Im(z))) * + sqrt((norm(z) - Re(z)) / &2))`;; + +let COMPLEX_NORM_GE_RE_IM = prove + (`!z. abs(Re(z)) <= norm(z) /\ abs(Im(z)) <= norm(z)`, + GEN_TAC THEN ONCE_REWRITE_TAC[GSYM POW_2_SQRT_ABS] THEN + REWRITE_TAC[complex_norm] THEN + CONJ_TAC THEN + MATCH_MP_TAC SQRT_MONO_LE THEN + ASM_SIMP_TAC[REAL_LE_ADDR; REAL_LE_ADDL; REAL_POW_2; REAL_LE_SQUARE]);; + +let CSQRT = prove + (`!z. csqrt(z) pow 2 = z`, + GEN_TAC THEN REWRITE_TAC[COMPLEX_POW_2; csqrt] THEN COND_CASES_TAC THENL + [COND_CASES_TAC THEN + ASM_REWRITE_TAC[CX_DEF; complex_mul; RE; IM; REAL_MUL_RZERO; REAL_MUL_LZERO; + REAL_SUB_LZERO; REAL_SUB_RZERO; REAL_ADD_LID; COMPLEX_EQ] THEN + REWRITE_TAC[REAL_NEG_EQ; GSYM REAL_POW_2] THEN + ASM_SIMP_TAC[SQRT_POW_2; REAL_ARITH `~(&0 <= x) ==> &0 <= --x`]; + ALL_TAC] THEN + REWRITE_TAC[complex_mul; RE; IM] THEN + ONCE_REWRITE_TAC[REAL_ARITH + `(s * s - (i * s') * (i * s') = s * s - (i * i) * (s' * s')) /\ + (s * i * s' + (i * s')* s = &2 * i * s * s')`] THEN + REWRITE_TAC[GSYM REAL_POW_2] THEN + SUBGOAL_THEN `&0 <= norm(z) + Re(z) /\ &0 <= norm(z) - Re(z)` + STRIP_ASSUME_TAC THENL + [MP_TAC(SPEC `z:complex` COMPLEX_NORM_GE_RE_IM) THEN REAL_ARITH_TAC; + ALL_TAC] THEN + ASM_SIMP_TAC[REAL_LE_DIV; REAL_POS; GSYM SQRT_MUL; SQRT_POW_2] THEN + REWRITE_TAC[COMPLEX_EQ; RE; IM] THEN CONJ_TAC THENL + [ASM_SIMP_TAC[REAL_POW_DIV; REAL_POW2_ABS; + REAL_POW_EQ_0; REAL_DIV_REFL] THEN + REWRITE_TAC[real_div; REAL_MUL_LID; GSYM REAL_SUB_RDISTRIB] THEN + REWRITE_TAC[REAL_ARITH `(m + r) - (m - r) = r * &2`] THEN + REWRITE_TAC[GSYM REAL_MUL_ASSOC] THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN + REWRITE_TAC[REAL_MUL_RID]; ALL_TAC] THEN + REWRITE_TAC[real_div] THEN + ONCE_REWRITE_TAC[AC REAL_MUL_AC + `(a * b) * a' * b = (a * a') * (b * b:real)`] THEN + REWRITE_TAC[REAL_DIFFSQ] THEN + REWRITE_TAC[complex_norm; GSYM REAL_POW_2] THEN + SIMP_TAC[SQRT_POW_2; REAL_LE_ADD; + REWRITE_RULE[GSYM REAL_POW_2] REAL_LE_SQUARE] THEN + REWRITE_TAC[REAL_ADD_SUB; GSYM REAL_POW_MUL] THEN + REWRITE_TAC[POW_2_SQRT_ABS] THEN + REWRITE_TAC[REAL_ABS_MUL; REAL_ABS_INV; REAL_ABS_NUM] THEN + ONCE_REWRITE_TAC[AC REAL_MUL_AC + `&2 * (i * a') * a * h = i * (&2 * h) * a * a'`] THEN + CONV_TAC REAL_RAT_REDUCE_CONV THEN + REWRITE_TAC[REAL_MUL_LID; GSYM real_div] THEN + ASM_SIMP_TAC[REAL_DIV_REFL; REAL_ABS_ZERO; REAL_MUL_RID]);; + +(* ------------------------------------------------------------------------- *) +(* Complex exponential. *) +(* ------------------------------------------------------------------------- *) + +let cexp = new_definition + `cexp z = Cx(exp(Re z)) * complex(cos(Im z),sin(Im z))`;; + +let CX_CEXP = prove + (`!x. Cx(exp x) = cexp(Cx x)`, + REWRITE_TAC[cexp; CX_DEF; RE; IM; SIN_0; COS_0] THEN + REWRITE_TAC[GSYM CX_DEF; GSYM CX_MUL; REAL_MUL_RID]);; + +let CEXP_0 = prove + (`cexp(Cx(&0)) = Cx(&1)`, + REWRITE_TAC[GSYM CX_CEXP; REAL_EXP_0]);; + +let CEXP_ADD = prove + (`!w z. cexp(w + z) = cexp(w) * cexp(z)`, + REWRITE_TAC[COMPLEX_EQ; cexp; complex_mul; complex_add; RE; IM; CX_DEF] THEN + REWRITE_TAC[REAL_EXP_ADD; SIN_ADD; COS_ADD] THEN CONV_TAC REAL_RING);; + +let CEXP_MUL = prove + (`!n z. cexp(Cx(&n) * z) = cexp(z) pow n`, + INDUCT_TAC THEN REWRITE_TAC[complex_pow] THEN + REWRITE_TAC[COMPLEX_MUL_LZERO; CEXP_0] THEN + REWRITE_TAC[GSYM REAL_OF_NUM_SUC; COMPLEX_ADD_RDISTRIB; CX_ADD] THEN + ASM_REWRITE_TAC[CEXP_ADD; COMPLEX_MUL_LID] THEN + REWRITE_TAC[COMPLEX_MUL_AC]);; + +let CEXP_NONZERO = prove + (`!z. ~(cexp z = Cx(&0))`, + GEN_TAC THEN REWRITE_TAC[cexp; COMPLEX_ENTIRE; CX_INJ; REAL_EXP_NZ] THEN + REWRITE_TAC[CX_DEF; RE; IM; COMPLEX_EQ] THEN + MP_TAC(SPEC `Im z` SIN_CIRCLE) THEN CONV_TAC REAL_RING);; + +let CEXP_NEG_LMUL = prove + (`!z. cexp(--z) * cexp(z) = Cx(&1)`, + REWRITE_TAC[GSYM CEXP_ADD; COMPLEX_ADD_LINV; CEXP_0]);; + +let CEXP_NEG_RMUL = prove + (`!z. cexp(z) * cexp(--z) = Cx(&1)`, + REWRITE_TAC[GSYM CEXP_ADD; COMPLEX_ADD_RINV; CEXP_0]);; + +let CEXP_NEG = prove + (`!z. cexp(--z) = inv(cexp z)`, + MESON_TAC[CEXP_NEG_LMUL; COMPLEX_MUL_LINV_UNIQ]);; + +let CEXP_SUB = prove + (`!w z. cexp(w - z) = cexp(w) / cexp(z)`, + REWRITE_TAC[complex_sub; complex_div; CEXP_NEG; CEXP_ADD]);; + +(* ------------------------------------------------------------------------- *) +(* Complex trig functions. *) +(* ------------------------------------------------------------------------- *) + +let ccos = new_definition + `ccos z = (cexp(ii * z) + cexp(--ii * z)) / Cx(&2)`;; + +let csin = new_definition + `csin z = (cexp(ii * z) - cexp(--ii * z)) / (Cx(&2) * ii)`;; + +let CX_CSIN,CX_CCOS = (CONJ_PAIR o prove) + (`(!x. Cx(sin x) = csin(Cx x)) /\ (!x. Cx(cos x) = ccos(Cx x))`, + REWRITE_TAC[csin; ccos; cexp; CX_DEF; ii; RE; IM; complex_mul; complex_add; + REAL_MUL_RZERO; REAL_MUL_LZERO; REAL_SUB_RZERO; + REAL_MUL_LID; complex_neg; REAL_EXP_0; REAL_ADD_LID; + REAL_MUL_LNEG; REAL_NEG_0; REAL_ADD_RID; complex_sub; + SIN_NEG; COS_NEG; GSYM REAL_MUL_2; GSYM real_sub; + complex_div; REAL_SUB_REFL; complex_inv; REAL_SUB_RNEG] THEN + CONJ_TAC THEN GEN_TAC THEN CONV_TAC REAL_RAT_REDUCE_CONV THEN + REWRITE_TAC[REAL_MUL_RZERO] THEN + AP_TERM_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN CONV_TAC REAL_RING);; + +let CSIN_0 = prove + (`csin(Cx(&0)) = Cx(&0)`, + REWRITE_TAC[GSYM CX_CSIN; SIN_0]);; + +let CCOS_0 = prove + (`ccos(Cx(&0)) = Cx(&1)`, + REWRITE_TAC[GSYM CX_CCOS; COS_0]);; + +let CSIN_CIRCLE = prove + (`!z. csin(z) pow 2 + ccos(z) pow 2 = Cx(&1)`, + GEN_TAC THEN REWRITE_TAC[csin; ccos] THEN + MP_TAC(SPEC `ii * z` CEXP_NEG_LMUL) THEN + MP_TAC COMPLEX_POW_II_2 THEN REWRITE_TAC[COMPLEX_MUL_LNEG] THEN + CONV_TAC COMPLEX_FIELD);; + +let CSIN_ADD = prove + (`!w z. csin(w + z) = csin(w) * ccos(z) + ccos(w) * csin(z)`, + REPEAT GEN_TAC THEN + REWRITE_TAC[csin; ccos; COMPLEX_ADD_LDISTRIB; CEXP_ADD] THEN + MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);; + +let CCOS_ADD = prove + (`!w z. ccos(w + z) = ccos(w) * ccos(z) - csin(w) * csin(z)`, + REPEAT GEN_TAC THEN + REWRITE_TAC[csin; ccos; COMPLEX_ADD_LDISTRIB; CEXP_ADD] THEN + MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_FIELD);; + +let CSIN_NEG = prove + (`!z. csin(--z) = --(csin(z))`, + REWRITE_TAC[csin; COMPLEX_MUL_LNEG; COMPLEX_MUL_RNEG; COMPLEX_NEG_NEG] THEN + GEN_TAC THEN MP_TAC COMPLEX_POW_II_2 THEN + CONV_TAC COMPLEX_FIELD);; + +let CCOS_NEG = prove + (`!z. ccos(--z) = ccos(z)`, + REWRITE_TAC[ccos; COMPLEX_MUL_LNEG; COMPLEX_MUL_RNEG; COMPLEX_NEG_NEG] THEN + GEN_TAC THEN MP_TAC COMPLEX_POW_II_2 THEN + CONV_TAC COMPLEX_FIELD);; + +let CSIN_DOUBLE = prove + (`!z. csin(Cx(&2) * z) = Cx(&2) * csin(z) * ccos(z)`, + REWRITE_TAC[COMPLEX_RING `Cx(&2) * x = x + x`; CSIN_ADD] THEN + CONV_TAC COMPLEX_RING);; + +let CCOS_DOUBLE = prove + (`!z. ccos(Cx(&2) * z) = (ccos(z) pow 2) - (csin(z) pow 2)`, + REWRITE_TAC[COMPLEX_RING `Cx(&2) * x = x + x`; CCOS_ADD] THEN + CONV_TAC COMPLEX_RING);; + +(* ------------------------------------------------------------------------- *) +(* Euler and de Moivre formulas. *) +(* ------------------------------------------------------------------------- *) + +let CEXP_EULER = prove + (`!z. cexp(ii * z) = ccos(z) + ii * csin(z)`, + REWRITE_TAC[ccos; csin] THEN MP_TAC COMPLEX_POW_II_2 THEN + CONV_TAC COMPLEX_FIELD);; + +let DEMOIVRE = prove + (`!z n. (ccos z + ii * csin z) pow n = + ccos(Cx(&n) * z) + ii * csin(Cx(&n) * z)`, + REWRITE_TAC[GSYM CEXP_EULER; GSYM CEXP_MUL] THEN + REWRITE_TAC[COMPLEX_MUL_AC]);; + +(* ------------------------------------------------------------------------- *) +(* Some lemmas. *) +(* ------------------------------------------------------------------------- *) + +let EXISTS_COMPLEX = prove + (`!P. (?z. P (Re z) (Im z)) <=> ?x y. P x y`, + MESON_TAC[RE; IM; COMPLEX]);; + +let COMPLEX_UNIMODULAR_POLAR = prove + (`!z. (norm z = &1) ==> ?x. z = complex(cos(x),sin(x))`, + GEN_TAC THEN + DISCH_THEN(MP_TAC o C AP_THM `2` o AP_TERM `(pow):real->num->real`) THEN + REWRITE_TAC[complex_norm] THEN + SIMP_TAC[REAL_POW_2; REWRITE_RULE[REAL_POW_2] SQRT_POW_2; + REAL_LE_SQUARE; REAL_LE_ADD] THEN + REWRITE_TAC[GSYM REAL_POW_2; REAL_MUL_LID] THEN + DISCH_THEN(X_CHOOSE_TAC `t:real` o MATCH_MP CIRCLE_SINCOS) THEN + EXISTS_TAC `t:real` THEN ASM_REWRITE_TAC[COMPLEX_EQ; RE; IM]);; + +let SIN_INTEGER_2PI = prove + (`!n. integer n ==> sin((&2 * pi) * n) = &0`, + REWRITE_TAC[integer; REAL_ARITH `abs(x) = &n <=> x = &n \/ x = -- &n`] THEN + REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[REAL_MUL_RNEG; SIN_NEG] THEN + REWRITE_TAC[GSYM REAL_MUL_ASSOC; SIN_DOUBLE] THEN + REWRITE_TAC[REAL_ARITH `pi * &n = &n * pi`; SIN_NPI] THEN + REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO; REAL_NEG_0]);; + +let COS_INTEGER_2PI = prove + (`!n. integer n ==> cos((&2 * pi) * n) = &1`, + REWRITE_TAC[integer; REAL_ARITH `abs(x) = &n <=> x = &n \/ x = -- &n`] THEN + REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[REAL_MUL_RNEG; COS_NEG] THEN + REWRITE_TAC[GSYM REAL_MUL_ASSOC; COS_DOUBLE] THEN + REWRITE_TAC[REAL_ARITH `pi * &n = &n * pi`; SIN_NPI; COS_NPI] THEN + REWRITE_TAC[REAL_POW_POW] THEN ONCE_REWRITE_TAC[MULT_SYM] THEN + REWRITE_TAC[GSYM REAL_POW_POW; REAL_POW_2] THEN + CONV_TAC REAL_RAT_REDUCE_CONV THEN + REWRITE_TAC[REAL_POW_ONE; REAL_SUB_RZERO]);; + +let SINCOS_PRINCIPAL_VALUE = prove + (`!x. ?y. (--pi < y /\ y <= pi) /\ (sin(y) = sin(x) /\ cos(y) = cos(x))`, + GEN_TAC THEN EXISTS_TAC `pi - (&2 * pi) * frac((pi - x) / (&2 * pi))` THEN + CONJ_TAC THENL + [SIMP_TAC[REAL_ARITH `--p < p - x <=> x < (&2 * p) * &1`; + REAL_ARITH `p - x <= p <=> (&2 * p) * &0 <= x`; + REAL_LT_LMUL_EQ; REAL_LE_LMUL_EQ; REAL_LT_MUL; + PI_POS; REAL_OF_NUM_LT; ARITH; FLOOR_FRAC]; + REWRITE_TAC[FRAC_FLOOR; REAL_SUB_LDISTRIB] THEN + SIMP_TAC[REAL_DIV_LMUL; REAL_ENTIRE; REAL_OF_NUM_EQ; ARITH; REAL_LT_IMP_NZ; + PI_POS; REAL_ARITH `a - (a - b - c):real = b + c`; SIN_ADD; COS_ADD] THEN + SIMP_TAC[FLOOR_FRAC; SIN_INTEGER_2PI; COS_INTEGER_2PI] THEN + CONV_TAC REAL_RING]);; + +(* ------------------------------------------------------------------------- *) +(* Complex logarithms (the conventional principal value). *) +(* ------------------------------------------------------------------------- *) + +let clog = new_definition + `clog z = @w. cexp(w) = z /\ --pi < Im(w) /\ Im(w) <= pi`;; + +let CLOG_WORKS = prove + (`!z. ~(z = Cx(&0)) + ==> cexp(clog z) = z /\ --pi < Im(clog z) /\ Im(clog z) <= pi`, + GEN_TAC THEN DISCH_TAC THEN REWRITE_TAC[clog] THEN CONV_TAC SELECT_CONV THEN + REWRITE_TAC[cexp; EXISTS_COMPLEX] THEN + EXISTS_TAC `ln(norm(z:complex))` THEN + SUBGOAL_THEN `exp(ln(norm(z:complex))) = norm(z)` SUBST1_TAC THENL + [ASM_MESON_TAC[REAL_EXP_LN; COMPLEX_NORM_NZ]; ALL_TAC] THEN + MP_TAC(SPEC `z / Cx(norm z)` COMPLEX_UNIMODULAR_POLAR) THEN ANTS_TAC THENL + [ASM_SIMP_TAC[COMPLEX_NORM_DIV; COMPLEX_NORM_CX] THEN + ASM_SIMP_TAC[COMPLEX_ABS_NORM; REAL_DIV_REFL; COMPLEX_NORM_ZERO]; + ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `x:real` STRIP_ASSUME_TAC) THEN + MP_TAC(SPEC `x:real` SINCOS_PRINCIPAL_VALUE) THEN + MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `y:real` THEN + STRIP_TAC THEN ASM_REWRITE_TAC[] THEN + ASM_MESON_TAC[CX_INJ; COMPLEX_DIV_LMUL; COMPLEX_NORM_ZERO]);; + +let CEXP_CLOG = prove + (`!z. ~(z = Cx(&0)) ==> cexp(clog z) = z`, + SIMP_TAC[CLOG_WORKS]);; + +(* ------------------------------------------------------------------------- *) +(* Unwinding number. *) +(* ------------------------------------------------------------------------- *) + +let unwinding = new_definition + `unwinding(z) = (z - clog(cexp z)) / (Cx(&2 * pi) * ii)`;; + +let COMPLEX_II_NZ = prove + (`~(ii = Cx(&0))`, + MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_RING);; + +let UNWINDING_2PI = prove + (`Cx(&2 * pi) * ii * unwinding(z) = z - clog(cexp z)`, + REWRITE_TAC[unwinding; COMPLEX_MUL_ASSOC] THEN + MATCH_MP_TAC COMPLEX_DIV_LMUL THEN + REWRITE_TAC[COMPLEX_ENTIRE; CX_INJ; COMPLEX_II_NZ] THEN + MP_TAC PI_POS THEN REAL_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* An example of how to get nice identities with unwinding number. *) +(* ------------------------------------------------------------------------- *) + +let CLOG_MUL = prove + (`!w z. ~(w = Cx(&0)) /\ ~(z = Cx(&0)) + ==> clog(w * z) = + clog(w) + clog(z) - + Cx(&2 * pi) * ii * unwinding(clog w + clog z)`, + REWRITE_TAC[UNWINDING_2PI; + COMPLEX_RING `w + z - ((w + z) - c) = c:complex`] THEN + ASM_SIMP_TAC[CEXP_ADD; CEXP_CLOG]);; diff --git a/Complex/complexnumbers.ml b/Complex/complexnumbers.ml new file mode 100644 index 0000000..ef2f172 --- /dev/null +++ b/Complex/complexnumbers.ml @@ -0,0 +1,912 @@ +(* ========================================================================= *) +(* Basic definitions and properties of complex numbers. *) +(* ========================================================================= *) + +needs "Library/transc.ml";; + +prioritize_real();; + +(* ------------------------------------------------------------------------- *) +(* Definition of complex number type. *) +(* ------------------------------------------------------------------------- *) + +let complex_tybij_raw = + new_type_definition "complex" ("complex","coords") + (prove(`?x:real#real. T`,REWRITE_TAC[]));; + +let complex_tybij = REWRITE_RULE [] complex_tybij_raw;; + +(* ------------------------------------------------------------------------- *) +(* Real and imaginary parts of a number. *) +(* ------------------------------------------------------------------------- *) + +let RE_DEF = new_definition + `Re(z) = FST(coords(z))`;; + +let IM_DEF = new_definition + `Im(z) = SND(coords(z))`;; + +(* ------------------------------------------------------------------------- *) +(* Set up overloading. *) +(* ------------------------------------------------------------------------- *) + +do_list overload_interface + ["+",`complex_add:complex->complex->complex`; + "-",`complex_sub:complex->complex->complex`; + "*",`complex_mul:complex->complex->complex`; + "/",`complex_div:complex->complex->complex`; + "--",`complex_neg:complex->complex`; + "pow",`complex_pow:complex->num->complex`; + "inv",`complex_inv:complex->complex`];; + +let prioritize_complex() = prioritize_overload(mk_type("complex",[]));; + +(* ------------------------------------------------------------------------- *) +(* Complex absolute value (modulus). *) +(* ------------------------------------------------------------------------- *) + +make_overloadable "norm" `:A->real`;; +overload_interface("norm",`complex_norm:complex->real`);; + +let complex_norm = new_definition + `norm(z) = sqrt(Re(z) pow 2 + Im(z) pow 2)`;; + +(* ------------------------------------------------------------------------- *) +(* Imaginary unit (too inconvenient to use "i"!) *) +(* ------------------------------------------------------------------------- *) + +let ii = new_definition + `ii = complex(&0,&1)`;; + +(* ------------------------------------------------------------------------- *) +(* Injection from reals. *) +(* ------------------------------------------------------------------------- *) + +let CX_DEF = new_definition + `Cx(a) = complex(a,&0)`;; + +(* ------------------------------------------------------------------------- *) +(* Arithmetic operations. *) +(* ------------------------------------------------------------------------- *) + +let complex_neg = new_definition + `--z = complex(--(Re(z)),--(Im(z)))`;; + +let complex_add = new_definition + `w + z = complex(Re(w) + Re(z),Im(w) + Im(z))`;; + +let complex_sub = new_definition + `w - z = w + --z`;; + +let complex_mul = new_definition + `w * z = complex(Re(w) * Re(z) - Im(w) * Im(z), + Re(w) * Im(z) + Im(w) * Re(z))`;; + +let complex_inv = new_definition + `inv(z) = complex(Re(z) / (Re(z) pow 2 + Im(z) pow 2), + --(Im(z)) / (Re(z) pow 2 + Im(z) pow 2))`;; + +let complex_div = new_definition + `w / z = w * inv(z)`;; + +let complex_pow = new_recursive_definition num_RECURSION + `(x pow 0 = Cx(&1)) /\ + (!n. x pow (SUC n) = x * x pow n)`;; + +(* ------------------------------------------------------------------------- *) +(* Various handy rewrites. *) +(* ------------------------------------------------------------------------- *) + +let RE = prove + (`(Re(complex(x,y)) = x)`, + REWRITE_TAC[RE_DEF; complex_tybij]);; + +let IM = prove + (`Im(complex(x,y)) = y`, + REWRITE_TAC[IM_DEF; complex_tybij]);; + +let COMPLEX = prove + (`complex(Re(z),Im(z)) = z`, + REWRITE_TAC[IM_DEF; RE_DEF; complex_tybij]);; + +let COMPLEX_EQ = prove + (`!w z. (w = z) <=> (Re(w) = Re(z)) /\ (Im(w) = Im(z))`, + REWRITE_TAC[RE_DEF; IM_DEF; GSYM PAIR_EQ] THEN MESON_TAC[complex_tybij]);; + +(* ------------------------------------------------------------------------- *) +(* Crude tactic to automate very simple algebraic equivalences. *) +(* ------------------------------------------------------------------------- *) + +let SIMPLE_COMPLEX_ARITH_TAC = + REWRITE_TAC[COMPLEX_EQ; RE; IM; CX_DEF; + complex_add; complex_neg; complex_sub; complex_mul] THEN + REAL_ARITH_TAC;; + +let SIMPLE_COMPLEX_ARITH tm = prove(tm,SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Basic algebraic properties that can be proved automatically by this. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_ADD_SYM = prove + (`!x y. x + y = y + x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_ASSOC = prove + (`!x y z. x + y + z = (x + y) + z`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_LID = prove + (`!x. Cx(&0) + x = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_LINV = prove + (`!x. --x + x = Cx(&0)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_SYM = prove + (`!x y. x * y = y * x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_ASSOC = prove + (`!x y z. x * y * z = (x * y) * z`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_LID = prove + (`!x. Cx(&1) * x = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_LDISTRIB = prove + (`!x y z. x * (y + z) = x * y + x * z`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_AC = prove + (`(m + n = n + m) /\ ((m + n) + p = m + n + p) /\ (m + n + p = n + m + p)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_AC = prove + (`(m * n = n * m) /\ ((m * n) * p = m * n * p) /\ (m * n * p = n * m * p)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_RID = prove + (`!x. x + Cx(&0) = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_RID = prove + (`!x. x * Cx(&1) = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_RINV = prove + (`!x. x + --x = Cx(&0)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_RDISTRIB = prove + (`!x y z. (x + y) * z = x * z + y * z`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_ADD_LCANCEL = prove + (`!x y z. (x + y = x + z) <=> (y = z)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_ADD_RCANCEL = prove + (`!x y z. (x + z = y + z) <=> (x = y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_RZERO = prove + (`!x. x * Cx(&0) = Cx(&0)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_LZERO = prove + (`!x. Cx(&0) * x = Cx(&0)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_NEG = prove + (`!x. --(--x) = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_RNEG = prove + (`!x y. x * --y = --(x * y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_LNEG = prove + (`!x y. --x * y = --(x * y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_ADD = prove + (`!x y. --(x + y) = --x + --y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_0 = prove + (`--Cx(&0) = Cx(&0)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_ADD_LCANCEL_0 = prove + (`!x y. (x + y = x) <=> (y = Cx(&0))`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_ADD_RCANCEL_0 = prove + (`!x y. (x + y = y) <=> (x = Cx(&0))`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_LNEG_UNIQ = prove + (`!x y. (x + y = Cx(&0)) <=> (x = --y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_RNEG_UNIQ = prove + (`!x y. (x + y = Cx(&0)) <=> (y = --x)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_LMUL = prove + (`!x y. --(x * y) = --x * y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_RMUL = prove + (`!x y. --(x * y) = x * --y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_MUL2 = prove + (`!x y. --x * --y = x * y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_ADD = prove + (`!x y. x - y + y = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_ADD2 = prove + (`!x y. y + x - y = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_REFL = prove + (`!x. x - x = Cx(&0)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_0 = prove + (`!x y. (x - y = Cx(&0)) <=> (x = y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_EQ_0 = prove + (`!x. (--x = Cx(&0)) <=> (x = Cx(&0))`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_SUB = prove + (`!x y. --(x - y) = y - x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_SUB = prove + (`!x y. (x + y) - x = y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_EQ = prove + (`!x y. (--x = y) <=> (x = --y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_NEG_MINUS1 = prove + (`!x. --x = --Cx(&1) * x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_SUB = prove + (`!x y. x - y - x = --y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD2_SUB2 = prove + (`!a b c d. (a + b) - (c + d) = a - c + b - d`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_LZERO = prove + (`!x. Cx(&0) - x = --x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_RZERO = prove + (`!x. x - Cx(&0) = x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_LNEG = prove + (`!x y. --x - y = --(x + y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_RNEG = prove + (`!x y. x - --y = x + y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_NEG2 = prove + (`!x y. --x - --y = y - x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_TRIANGLE = prove + (`!a b c. a - b + b - c = a - c`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_SUB_LADD = prove + (`!x y z. (x = y - z) <=> (x + z = y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_SUB_RADD = prove + (`!x y z. (x - y = z) <=> (x = z + y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_SUB2 = prove + (`!x y. x - (x - y) = y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_ADD_SUB2 = prove + (`!x y. x - (x + y) = --y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_DIFFSQ = prove + (`!x y. (x + y) * (x - y) = x * x - y * y`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_EQ_NEG2 = prove + (`!x y. (--x = --y) <=> (x = y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_LDISTRIB = prove + (`!x y z. x * (y - z) = x * y - x * z`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_SUB_RDISTRIB = prove + (`!x y z. (x - y) * z = x * z - y * z`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_MUL_2 = prove + (`!x. &2 * x = x + x`, + SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Homomorphic embedding properties for Cx mapping. *) +(* ------------------------------------------------------------------------- *) + +let CX_INJ = prove + (`!x y. (Cx(x) = Cx(y)) <=> (x = y)`, + REWRITE_TAC[CX_DEF; COMPLEX_EQ; RE; IM]);; + +let CX_NEG = prove + (`!x. Cx(--x) = --(Cx(x))`, + REWRITE_TAC[CX_DEF; complex_neg; RE; IM; REAL_NEG_0]);; + +let CX_INV = prove + (`!x. Cx(inv x) = inv(Cx x)`, + GEN_TAC THEN + REWRITE_TAC[CX_DEF; complex_inv; RE; IM] THEN + REWRITE_TAC[real_div; REAL_NEG_0; REAL_MUL_LZERO] THEN + REWRITE_TAC[COMPLEX_EQ; REAL_POW_2; REAL_MUL_RZERO; RE; IM] THEN + REWRITE_TAC[REAL_ADD_RID; REAL_INV_MUL] THEN + ASM_CASES_TAC `x = &0` THEN + ASM_REWRITE_TAC[REAL_INV_0; REAL_MUL_LZERO] THEN + REWRITE_TAC[REAL_MUL_ASSOC] THEN + GEN_REWRITE_TAC LAND_CONV [GSYM REAL_MUL_LID] THEN + AP_THM_TAC THEN AP_TERM_TAC THEN ASM_MESON_TAC[REAL_MUL_RINV]);; + +let CX_ADD = prove + (`!x y. Cx(x + y) = Cx(x) + Cx(y)`, + REWRITE_TAC[CX_DEF; complex_add; RE; IM; REAL_ADD_LID]);; + +let CX_SUB = prove + (`!x y. Cx(x - y) = Cx(x) - Cx(y)`, + REWRITE_TAC[complex_sub; real_sub; CX_ADD; CX_NEG]);; + +let CX_MUL = prove + (`!x y. Cx(x * y) = Cx(x) * Cx(y)`, + REWRITE_TAC[CX_DEF; complex_mul; RE; IM; REAL_MUL_LZERO; REAL_MUL_RZERO] THEN + REWRITE_TAC[REAL_SUB_RZERO; REAL_ADD_RID]);; + +let CX_DIV = prove + (`!x y. Cx(x / y) = Cx(x) / Cx(y)`, + REWRITE_TAC[complex_div; real_div; CX_MUL; CX_INV]);; + +let CX_POW = prove + (`!x n. Cx(x pow n) = Cx(x) pow n`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; real_pow; CX_MUL]);; + +let CX_ABS = prove + (`!x. Cx(abs x) = Cx(norm(Cx(x)))`, + REWRITE_TAC[CX_DEF; complex_norm; COMPLEX_EQ; RE; IM] THEN + REWRITE_TAC[REAL_POW_2; REAL_MUL_LZERO; REAL_ADD_RID] THEN + REWRITE_TAC[GSYM REAL_POW_2; POW_2_SQRT_ABS]);; + +let COMPLEX_NORM_CX = prove + (`!x. norm(Cx(x)) = abs(x)`, + REWRITE_TAC[GSYM CX_INJ; CX_ABS]);; + +(* ------------------------------------------------------------------------- *) +(* A convenient lemma that we need a few times below. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_ENTIRE = prove + (`!x y. (x * y = Cx(&0)) <=> (x = Cx(&0)) \/ (y = Cx(&0))`, + REWRITE_TAC[COMPLEX_EQ; complex_mul; RE; IM; CX_DEF; GSYM REAL_SOS_EQ_0] THEN + CONV_TAC REAL_RING);; + +(* ------------------------------------------------------------------------- *) +(* Powers. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_POW_ADD = prove + (`!x m n. x pow (m + n) = x pow m * x pow n`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[ADD_CLAUSES; complex_pow; + COMPLEX_MUL_LID; COMPLEX_MUL_ASSOC]);; + +let COMPLEX_POW_POW = prove + (`!x m n. (x pow m) pow n = x pow (m * n)`, + GEN_TAC THEN GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; MULT_CLAUSES; COMPLEX_POW_ADD]);; + +let COMPLEX_POW_1 = prove + (`!x. x pow 1 = x`, + REWRITE_TAC[num_CONV `1`] THEN REWRITE_TAC[complex_pow; COMPLEX_MUL_RID]);; + +let COMPLEX_POW_2 = prove + (`!x. x pow 2 = x * x`, + REWRITE_TAC[num_CONV `2`] THEN REWRITE_TAC[complex_pow; COMPLEX_POW_1]);; + +let COMPLEX_POW_NEG = prove + (`!x n. (--x) pow n = if EVEN n then x pow n else --(x pow n)`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; EVEN] THEN + ASM_CASES_TAC `EVEN n` THEN + ASM_REWRITE_TAC[COMPLEX_MUL_RNEG; COMPLEX_MUL_LNEG; COMPLEX_NEG_NEG]);; + +let COMPLEX_POW_ONE = prove + (`!n. Cx(&1) pow n = Cx(&1)`, + INDUCT_TAC THEN ASM_REWRITE_TAC[complex_pow; COMPLEX_MUL_LID]);; + +let COMPLEX_POW_MUL = prove + (`!x y n. (x * y) pow n = (x pow n) * (y pow n)`, + GEN_TAC THEN GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; COMPLEX_MUL_LID; COMPLEX_MUL_AC]);; + +let COMPLEX_POW_II_2 = prove + (`ii pow 2 = --Cx(&1)`, + REWRITE_TAC[ii; COMPLEX_POW_2; complex_mul; CX_DEF; RE; IM; complex_neg] THEN + CONV_TAC REAL_RAT_REDUCE_CONV);; + +let COMPLEX_POW_EQ_0 = prove + (`!x n. (x pow n = Cx(&0)) <=> (x = Cx(&0)) /\ ~(n = 0)`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[NOT_SUC; complex_pow; COMPLEX_ENTIRE] THENL + [SIMPLE_COMPLEX_ARITH_TAC; CONV_TAC TAUT]);; + +(* ------------------------------------------------------------------------- *) +(* Norms (aka "moduli"). *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_NORM_CX = prove + (`!x. norm(Cx x) = abs(x)`, + GEN_TAC THEN REWRITE_TAC[complex_norm; CX_DEF; RE; IM] THEN + REWRITE_TAC[REAL_POW_2; REAL_MUL_LZERO; REAL_ADD_RID] THEN + REWRITE_TAC[GSYM REAL_POW_2; POW_2_SQRT_ABS]);; + +let COMPLEX_NORM_POS = prove + (`!z. &0 <= norm(z)`, + SIMP_TAC[complex_norm; SQRT_POS_LE; REAL_POW_2; + REAL_LE_SQUARE; REAL_LE_ADD]);; + +let COMPLEX_ABS_NORM = prove + (`!z. abs(norm z) = norm z`, + REWRITE_TAC[real_abs; COMPLEX_NORM_POS]);; + +let COMPLEX_NORM_ZERO = prove + (`!z. (norm z = &0) <=> (z = Cx(&0))`, + GEN_TAC THEN REWRITE_TAC[complex_norm] THEN + GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM SQRT_0] THEN + SIMP_TAC[REAL_POW_2; REAL_LE_SQUARE; REAL_LE_ADD; REAL_POS; SQRT_INJ] THEN + REWRITE_TAC[COMPLEX_EQ; RE; IM; CX_DEF] THEN + SIMP_TAC[REAL_LE_SQUARE; REAL_ARITH + `&0 <= x /\ &0 <= y ==> ((x + y = &0) <=> (x = &0) /\ (y = &0))`] THEN + REWRITE_TAC[REAL_ENTIRE]);; + +let COMPLEX_NORM_NUM = prove + (`norm(Cx(&n)) = &n`, + REWRITE_TAC[COMPLEX_NORM_CX; REAL_ABS_NUM]);; + +let COMPLEX_NORM_0 = prove + (`norm(Cx(&0)) = &0`, + MESON_TAC[COMPLEX_NORM_ZERO]);; + +let COMPLEX_NORM_NZ = prove + (`!z. &0 < norm(z) <=> ~(z = Cx(&0))`, + MESON_TAC[COMPLEX_NORM_ZERO; COMPLEX_ABS_NORM; REAL_ABS_NZ]);; + +let COMPLEX_NORM_NEG = prove + (`!z. norm(--z) = norm(z)`, + REWRITE_TAC[complex_neg; complex_norm; REAL_POW_2; RE; IM] THEN + GEN_TAC THEN AP_TERM_TAC THEN REAL_ARITH_TAC);; + +let COMPLEX_NORM_MUL = prove + (`!w z. norm(w * z) = norm(w) * norm(z)`, + REPEAT GEN_TAC THEN + REWRITE_TAC[complex_norm; complex_mul; RE; IM] THEN + SIMP_TAC[GSYM SQRT_MUL; REAL_POW_2; REAL_LE_ADD; REAL_LE_SQUARE] THEN + AP_TERM_TAC THEN REAL_ARITH_TAC);; + +let COMPLEX_NORM_POW = prove + (`!z n. norm(z pow n) = norm(z) pow n`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; real_pow; COMPLEX_NORM_NUM; COMPLEX_NORM_MUL]);; + +let COMPLEX_NORM_INV = prove + (`!z. norm(inv z) = inv(norm z)`, + GEN_TAC THEN REWRITE_TAC[complex_norm; complex_inv; RE; IM] THEN + REWRITE_TAC[REAL_POW_2; real_div] THEN + REWRITE_TAC[REAL_ARITH `(r * d) * r * d + (--i * d) * --i * d = + (r * r + i * i) * d * d:real`] THEN + ASM_CASES_TAC `Re z * Re z + Im z * Im z = &0` THENL + [ASM_REWRITE_TAC[REAL_INV_0; SQRT_0; REAL_MUL_LZERO]; ALL_TAC] THEN + CONV_TAC SYM_CONV THEN MATCH_MP_TAC REAL_MUL_RINV_UNIQ THEN + SIMP_TAC[GSYM SQRT_MUL; REAL_LE_MUL; REAL_LE_INV_EQ; REAL_LE_ADD; + REAL_LE_SQUARE] THEN + ONCE_REWRITE_TAC[AC REAL_MUL_AC + `a * a * b * b:real = (a * b) * (a * b)`] THEN + ASM_SIMP_TAC[REAL_MUL_RINV; REAL_MUL_LID; SQRT_1]);; + +let COMPLEX_NORM_DIV = prove + (`!w z. norm(w / z) = norm(w) / norm(z)`, + REWRITE_TAC[complex_div; real_div; COMPLEX_NORM_INV; COMPLEX_NORM_MUL]);; + +let COMPLEX_NORM_TRIANGLE = prove + (`!w z. norm(w + z) <= norm(w) + norm(z)`, + REPEAT GEN_TAC THEN REWRITE_TAC[complex_norm; complex_add; RE; IM] THEN + MATCH_MP_TAC(REAL_ARITH `&0 <= y /\ abs(x) <= abs(y) ==> x <= y`) THEN + SIMP_TAC[SQRT_POS_LE; REAL_POW_2; REAL_LE_ADD; REAL_LE_SQUARE; + REAL_LE_SQUARE_ABS; SQRT_POW_2] THEN + GEN_REWRITE_TAC RAND_CONV[REAL_ARITH + `(a + b) * (a + b) = a * a + b * b + &2 * a * b`] THEN + REWRITE_TAC[GSYM REAL_POW_2] THEN + SIMP_TAC[SQRT_POW_2; REAL_POW_2; REAL_LE_ADD; REAL_LE_SQUARE] THEN + REWRITE_TAC[REAL_ARITH + `(rw + rz) * (rw + rz) + (iw + iz) * (iw + iz) <= + (rw * rw + iw * iw) + (rz * rz + iz * iz) + &2 * other <=> + rw * rz + iw * iz <= other`] THEN + SIMP_TAC[GSYM SQRT_MUL; REAL_POW_2; REAL_LE_ADD; REAL_LE_SQUARE] THEN + MATCH_MP_TAC(REAL_ARITH `&0 <= y /\ abs(x) <= abs(y) ==> x <= y`) THEN + SIMP_TAC[SQRT_POS_LE; REAL_POW_2; REAL_LE_ADD; REAL_LE_SQUARE; + REAL_LE_SQUARE_ABS; SQRT_POW_2; REAL_LE_MUL] THEN + REWRITE_TAC[REAL_ARITH + `(rw * rz + iw * iz) * (rw * rz + iw * iz) <= + (rw * rw + iw * iw) * (rz * rz + iz * iz) <=> + &0 <= (rw * iz - rz * iw) * (rw * iz - rz * iw)`] THEN + REWRITE_TAC[REAL_LE_SQUARE]);; + +let COMPLEX_NORM_TRIANGLE_SUB = prove + (`!w z. norm(w) <= norm(w + z) + norm(z)`, + MESON_TAC[COMPLEX_NORM_TRIANGLE; COMPLEX_NORM_NEG; COMPLEX_ADD_ASSOC; + COMPLEX_ADD_RINV; COMPLEX_ADD_RID]);; + +let COMPLEX_NORM_ABS_NORM = prove + (`!w z. abs(norm w - norm z) <= norm(w - z)`, + REPEAT GEN_TAC THEN + MATCH_MP_TAC(REAL_ARITH + `a - b <= x /\ b - a <= x ==> abs(a - b) <= x:real`) THEN + MESON_TAC[COMPLEX_NEG_SUB; COMPLEX_NORM_NEG; REAL_LE_SUB_RADD; complex_sub; + COMPLEX_NORM_TRIANGLE_SUB]);; + +(* ------------------------------------------------------------------------- *) +(* Complex conjugate. *) +(* ------------------------------------------------------------------------- *) + +let cnj = new_definition + `cnj(z) = complex(Re(z),--(Im(z)))`;; + +(* ------------------------------------------------------------------------- *) +(* Conjugation is an automorphism. *) +(* ------------------------------------------------------------------------- *) + +let CNJ_INJ = prove + (`!w z. (cnj(w) = cnj(z)) <=> (w = z)`, + REWRITE_TAC[cnj; COMPLEX_EQ; RE; IM; REAL_EQ_NEG2]);; + +let CNJ_CNJ = prove + (`!z. cnj(cnj z) = z`, + REWRITE_TAC[cnj; COMPLEX_EQ; RE; IM; REAL_NEG_NEG]);; + +let CNJ_CX = prove + (`!x. cnj(Cx x) = Cx x`, + REWRITE_TAC[cnj; COMPLEX_EQ; CX_DEF; REAL_NEG_0; RE; IM]);; + +let COMPLEX_NORM_CNJ = prove + (`!z. norm(cnj z) = norm(z)`, + REWRITE_TAC[complex_norm; cnj; REAL_POW_2] THEN + REWRITE_TAC[REAL_MUL_LNEG; REAL_MUL_RNEG; RE; IM; REAL_NEG_NEG]);; + +let CNJ_NEG = prove + (`!z. cnj(--z) = --(cnj z)`, + REWRITE_TAC[cnj; complex_neg; COMPLEX_EQ; RE; IM]);; + +let CNJ_INV = prove + (`!z. cnj(inv z) = inv(cnj z)`, + REWRITE_TAC[cnj; complex_inv; COMPLEX_EQ; RE; IM] THEN + REWRITE_TAC[real_div; REAL_NEG_NEG; REAL_POW_2; + REAL_MUL_LNEG; REAL_MUL_RNEG]);; + +let CNJ_ADD = prove + (`!w z. cnj(w + z) = cnj(w) + cnj(z)`, + REWRITE_TAC[cnj; complex_add; COMPLEX_EQ; RE; IM] THEN + REWRITE_TAC[REAL_NEG_ADD; REAL_MUL_LNEG; REAL_MUL_RNEG; REAL_NEG_NEG]);; + +let CNJ_SUB = prove + (`!w z. cnj(w - z) = cnj(w) - cnj(z)`, + REWRITE_TAC[complex_sub; CNJ_ADD; CNJ_NEG]);; + +let CNJ_MUL = prove + (`!w z. cnj(w * z) = cnj(w) * cnj(z)`, + REWRITE_TAC[cnj; complex_mul; COMPLEX_EQ; RE; IM] THEN + REWRITE_TAC[REAL_NEG_ADD; REAL_MUL_LNEG; REAL_MUL_RNEG; REAL_NEG_NEG]);; + +let CNJ_DIV = prove + (`!w z. cnj(w / z) = cnj(w) / cnj(z)`, + REWRITE_TAC[complex_div; CNJ_MUL; CNJ_INV]);; + +let CNJ_POW = prove + (`!z n. cnj(z pow n) = cnj(z) pow n`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; CNJ_MUL; CNJ_CX]);; + +(* ------------------------------------------------------------------------- *) +(* Conversion of (complex-type) rational constant to ML rational number. *) +(* ------------------------------------------------------------------------- *) + +let is_complex_const = + let cx_tm = `Cx` in + fun tm -> + is_comb tm & + let l,r = dest_comb tm in l = cx_tm & is_ratconst r;; + +let dest_complex_const = + let cx_tm = `Cx` in + fun tm -> + let l,r = dest_comb tm in + if l = cx_tm then rat_of_term r + else failwith "dest_complex_const";; + +let mk_complex_const = + let cx_tm = `Cx` in + fun r -> + mk_comb(cx_tm,term_of_rat r);; + +(* ------------------------------------------------------------------------- *) +(* Conversions to perform operations if coefficients are rational constants. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_RAT_MUL_CONV = + GEN_REWRITE_CONV I [GSYM CX_MUL] THENC RAND_CONV REAL_RAT_MUL_CONV;; + +let COMPLEX_RAT_ADD_CONV = + GEN_REWRITE_CONV I [GSYM CX_ADD] THENC RAND_CONV REAL_RAT_ADD_CONV;; + +let COMPLEX_RAT_EQ_CONV = + GEN_REWRITE_CONV I [CX_INJ] THENC REAL_RAT_EQ_CONV;; + +let COMPLEX_RAT_POW_CONV = + let x_tm = `x:real` + and n_tm = `n:num` in + let pth = SYM(SPECL [x_tm; n_tm] CX_POW) in + fun tm -> + let lop,r = dest_comb tm in + let op,bod = dest_comb lop in + let th1 = INST [rand bod,x_tm; r,n_tm] pth in + let tm1,tm2 = dest_comb(concl th1) in + if rand tm1 <> tm then failwith "COMPLEX_RAT_POW_CONV" else + let tm3,tm4 = dest_comb tm2 in + TRANS th1 (AP_TERM tm3 (REAL_RAT_REDUCE_CONV tm4));; + +(* ------------------------------------------------------------------------- *) +(* Instantiate polynomial normalizer. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_POLY_CLAUSES = prove + (`(!x y z. x + (y + z) = (x + y) + z) /\ + (!x y. x + y = y + x) /\ + (!x. Cx(&0) + x = x) /\ + (!x y z. x * (y * z) = (x * y) * z) /\ + (!x y. x * y = y * x) /\ + (!x. Cx(&1) * x = x) /\ + (!x. Cx(&0) * x = Cx(&0)) /\ + (!x y z. x * (y + z) = x * y + x * z) /\ + (!x. x pow 0 = Cx(&1)) /\ + (!x n. x pow (SUC n) = x * x pow n)`, + REWRITE_TAC[complex_pow] THEN SIMPLE_COMPLEX_ARITH_TAC) +and COMPLEX_POLY_NEG_CLAUSES = prove + (`(!x. --x = Cx(-- &1) * x) /\ + (!x y. x - y = x + Cx(-- &1) * y)`, + SIMPLE_COMPLEX_ARITH_TAC);; + +let COMPLEX_POLY_NEG_CONV,COMPLEX_POLY_ADD_CONV,COMPLEX_POLY_SUB_CONV, + COMPLEX_POLY_MUL_CONV,COMPLEX_POLY_POW_CONV,COMPLEX_POLY_CONV = + SEMIRING_NORMALIZERS_CONV COMPLEX_POLY_CLAUSES COMPLEX_POLY_NEG_CLAUSES + (is_complex_const, + COMPLEX_RAT_ADD_CONV,COMPLEX_RAT_MUL_CONV,COMPLEX_RAT_POW_CONV) + (<);; + +let COMPLEX_RAT_INV_CONV = + GEN_REWRITE_CONV I [GSYM CX_INV] THENC RAND_CONV REAL_RAT_INV_CONV;; + +let COMPLEX_POLY_CONV = + let neg_tm = `(--):complex->complex` + and inv_tm = `inv:complex->complex` + and add_tm = `(+):complex->complex->complex` + and sub_tm = `(-):complex->complex->complex` + and mul_tm = `(*):complex->complex->complex` + and div_tm = `(/):complex->complex->complex` + and pow_tm = `(pow):complex->num->complex` + and div_conv = REWR_CONV complex_div in + let rec COMPLEX_POLY_CONV tm = + if not(is_comb tm) or is_complex_const tm then REFL tm else + let lop,r = dest_comb tm in + if lop = neg_tm then + let th1 = AP_TERM lop (COMPLEX_POLY_CONV r) in + TRANS th1 (COMPLEX_POLY_NEG_CONV (rand(concl th1))) + else if lop = inv_tm then + let th1 = AP_TERM lop (COMPLEX_POLY_CONV r) in + TRANS th1 (TRY_CONV COMPLEX_RAT_INV_CONV (rand(concl th1))) + else if not(is_comb lop) then REFL tm else + let op,l = dest_comb lop in + if op = pow_tm then + let th1 = AP_THM (AP_TERM op (COMPLEX_POLY_CONV l)) r in + TRANS th1 (TRY_CONV COMPLEX_POLY_POW_CONV (rand(concl th1))) + else if op = add_tm or op = mul_tm or op = sub_tm then + let th1 = MK_COMB(AP_TERM op (COMPLEX_POLY_CONV l), + COMPLEX_POLY_CONV r) in + let fn = if op = add_tm then COMPLEX_POLY_ADD_CONV + else if op = mul_tm then COMPLEX_POLY_MUL_CONV + else COMPLEX_POLY_SUB_CONV in + TRANS th1 (fn (rand(concl th1))) + else if op = div_tm then + let th1 = div_conv tm in + TRANS th1 (COMPLEX_POLY_CONV (rand(concl th1))) + else REFL tm in + COMPLEX_POLY_CONV;; + +(* ------------------------------------------------------------------------- *) +(* Complex number version of usual ring procedure. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_MUL_LINV = prove + (`!z. ~(z = Cx(&0)) ==> (inv(z) * z = Cx(&1))`, + REWRITE_TAC[complex_mul; complex_inv; RE; IM; COMPLEX_EQ; CX_DEF] THEN + REWRITE_TAC[GSYM REAL_SOS_EQ_0] THEN CONV_TAC REAL_FIELD);; + +let COMPLEX_MUL_RINV = prove + (`!z. ~(z = Cx(&0)) ==> (z * inv(z) = Cx(&1))`, + ONCE_REWRITE_TAC[COMPLEX_MUL_SYM] THEN REWRITE_TAC[COMPLEX_MUL_LINV]);; + +let COMPLEX_RING,complex_ideal_cofactors = + let ring_pow_tm = `(pow):complex->num->complex` + and COMPLEX_INTEGRAL = prove + (`(!x. Cx(&0) * x = Cx(&0)) /\ + (!x y z. (x + y = x + z) <=> (y = z)) /\ + (!w x y z. (w * y + x * z = w * z + x * y) <=> (w = x) \/ (y = z))`, + REWRITE_TAC[COMPLEX_ENTIRE; SIMPLE_COMPLEX_ARITH + `(w * y + x * z = w * z + x * y) <=> + (w - x) * (y - z) = Cx(&0)`] THEN + SIMPLE_COMPLEX_ARITH_TAC) + and COMPLEX_RABINOWITSCH = prove + (`!x y:complex. ~(x = y) <=> ?z. (x - y) * z = Cx(&1)`, + REPEAT GEN_TAC THEN + GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM COMPLEX_SUB_0] THEN + MESON_TAC[COMPLEX_MUL_RINV; COMPLEX_MUL_LZERO; + SIMPLE_COMPLEX_ARITH `~(Cx(&1) = Cx(&0))`]) + and init = ALL_CONV in + let pure,ideal = + RING_AND_IDEAL_CONV + (dest_complex_const,mk_complex_const,COMPLEX_RAT_EQ_CONV, + `(--):complex->complex`,`(+):complex->complex->complex`, + `(-):complex->complex->complex`,`(inv):complex->complex`, + `(*):complex->complex->complex`,`(/):complex->complex->complex`, + `(pow):complex->num->complex`, + COMPLEX_INTEGRAL,COMPLEX_RABINOWITSCH,COMPLEX_POLY_CONV) in + (fun tm -> let th = init tm in EQ_MP (SYM th) (pure(rand(concl th)))), + ideal;; + +(* ------------------------------------------------------------------------- *) +(* Most basic properties of inverses. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_INV_0 = prove + (`inv(Cx(&0)) = Cx(&0)`, + REWRITE_TAC[complex_inv; CX_DEF; RE; IM; real_div; REAL_MUL_LZERO; + REAL_NEG_0]);; + +let COMPLEX_INV_MUL = prove + (`!w z. inv(w * z) = inv(w) * inv(z)`, + REPEAT GEN_TAC THEN + MAP_EVERY ASM_CASES_TAC [`w = Cx(&0)`; `z = Cx(&0)`] THEN + ASM_REWRITE_TAC[COMPLEX_INV_0; COMPLEX_MUL_LZERO; COMPLEX_MUL_RZERO] THEN + REPEAT(POP_ASSUM MP_TAC) THEN + REWRITE_TAC[complex_mul; complex_inv; RE; IM; COMPLEX_EQ; CX_DEF] THEN + REWRITE_TAC[GSYM REAL_SOS_EQ_0] THEN CONV_TAC REAL_FIELD);; + +let COMPLEX_INV_1 = prove + (`inv(Cx(&1)) = Cx(&1)`, + REWRITE_TAC[complex_inv; CX_DEF; RE; IM] THEN + CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[REAL_DIV_1]);; + +let COMPLEX_POW_INV = prove + (`!x n. (inv x) pow n = inv(x pow n)`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[complex_pow; COMPLEX_INV_1; COMPLEX_INV_MUL]);; + +let COMPLEX_INV_INV = prove + (`!x:complex. inv(inv x) = x`, + GEN_TAC THEN ASM_CASES_TAC `x = Cx(&0)` THEN + ASM_REWRITE_TAC[COMPLEX_INV_0] THEN + POP_ASSUM MP_TAC THEN + MAP_EVERY (fun t -> MP_TAC(SPEC t COMPLEX_MUL_RINV)) + [`x:complex`; `inv(x):complex`] THEN + CONV_TAC COMPLEX_RING);; + +(* ------------------------------------------------------------------------- *) +(* And also field procedure. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_FIELD = + let prenex_conv = + TOP_DEPTH_CONV BETA_CONV THENC + PURE_REWRITE_CONV[FORALL_SIMP; EXISTS_SIMP; complex_div; + COMPLEX_INV_INV; COMPLEX_INV_MUL; GSYM REAL_POW_INV] THENC + NNFC_CONV THENC DEPTH_BINOP_CONV `(/\)` CONDS_CELIM_CONV THENC + PRENEX_CONV + and setup_conv = NNF_CONV THENC WEAK_CNF_CONV THENC CONJ_CANON_CONV + and is_inv = + let inv_tm = `inv:complex->complex` + and is_div = is_binop `(/):complex->complex->complex` in + fun tm -> (is_div tm or (is_comb tm & rator tm = inv_tm)) & + not(is_complex_const(rand tm)) + and lemma_inv = MESON[COMPLEX_MUL_RINV] + `!x. x = Cx(&0) \/ x * inv(x) = Cx(&1)` + and dcases = MATCH_MP(TAUT `(p \/ q) /\ (r \/ s) ==> (p \/ r) \/ q /\ s`) in + let cases_rule th1 th2 = dcases (CONJ th1 th2) in + let BASIC_COMPLEX_FIELD tm = + let is_freeinv t = is_inv t & free_in t tm in + let itms = setify(map rand (find_terms is_freeinv tm)) in + let dth = if itms = [] then TRUTH + else end_itlist cases_rule (map (C SPEC lemma_inv) itms) in + let tm' = mk_imp(concl dth,tm) in + let th1 = setup_conv tm' in + let ths = map COMPLEX_RING (conjuncts(rand(concl th1))) in + let th2 = EQ_MP (SYM th1) (end_itlist CONJ ths) in + MP (EQ_MP (SYM th1) (end_itlist CONJ ths)) dth in + fun tm -> + let th0 = prenex_conv tm in + let tm0 = rand(concl th0) in + let avs,bod = strip_forall tm0 in + let th1 = setup_conv bod in + let ths = map BASIC_COMPLEX_FIELD (conjuncts(rand(concl th1))) in + EQ_MP (SYM th0) (GENL avs (EQ_MP (SYM th1) (end_itlist CONJ ths)));; + +(* ------------------------------------------------------------------------- *) +(* Properties of inverses, divisions are now mostly automatic. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_POW_DIV = prove + (`!x y n. (x / y) pow n = (x pow n) / (y pow n)`, + REWRITE_TAC[complex_div; COMPLEX_POW_MUL; COMPLEX_POW_INV]);; + +let COMPLEX_DIV_REFL = prove + (`!x. ~(x = Cx(&0)) ==> (x / x = Cx(&1))`, + CONV_TAC COMPLEX_FIELD);; + +let COMPLEX_EQ_MUL_LCANCEL = prove + (`!x y z. (x * y = x * z) <=> (x = Cx(&0)) \/ (y = z)`, + CONV_TAC COMPLEX_FIELD);; + +let COMPLEX_EQ_MUL_RCANCEL = prove + (`!x y z. (x * z = y * z) <=> (x = y) \/ (z = Cx(&0))`, + CONV_TAC COMPLEX_FIELD);; + +let COMPLEX_MUL_RINV_UNIQ = prove + (`!w z. w * z = Cx(&1) ==> inv w = z`, + CONV_TAC COMPLEX_FIELD);; + +let COMPLEX_MUL_LINV_UNIQ = prove + (`!w z. w * z = Cx(&1) ==> inv z = w`, + CONV_TAC COMPLEX_FIELD);; + +let COMPLEX_DIV_LMUL = prove + (`!w z. ~(z = Cx(&0)) ==> z * w / z = w`, + CONV_TAC COMPLEX_FIELD);; + +let COMPLEX_DIV_RMUL = prove + (`!w z. ~(z = Cx(&0)) ==> w / z * z = w`, + CONV_TAC COMPLEX_FIELD);; diff --git a/Complex/cpoly.ml b/Complex/cpoly.ml new file mode 100644 index 0000000..5b84dc0 --- /dev/null +++ b/Complex/cpoly.ml @@ -0,0 +1,977 @@ +(* ========================================================================= *) +(* Properties of complex polynomials (not canonically represented). *) +(* ========================================================================= *) + +needs "Complex/complexnumbers.ml";; + +prioritize_complex();; + +parse_as_infix("++",(16,"right"));; +parse_as_infix("**",(20,"right"));; +parse_as_infix("##",(20,"right"));; +parse_as_infix("divides",(14,"right"));; +parse_as_infix("exp",(22,"right"));; + +do_list override_interface + ["++",`poly_add:complex list->complex list->complex list`; + "**",`poly_mul:complex list->complex list->complex list`; + "##",`poly_cmul:complex->complex list->complex list`; + "neg",`poly_neg:complex list->complex list`; + "divides",`poly_divides:complex list->complex list->bool`; + "exp",`poly_exp:complex list -> num -> complex list`; + "diff",`poly_diff:complex list->complex list`];; + +let SIMPLE_COMPLEX_ARITH tm = prove(tm,SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Polynomials. *) +(* ------------------------------------------------------------------------- *) + +let poly = new_recursive_definition list_RECURSION + `(poly [] x = Cx(&0)) /\ + (poly (CONS h t) x = h + x * poly t x)`;; + +(* ------------------------------------------------------------------------- *) +(* Arithmetic operations on polynomials. *) +(* ------------------------------------------------------------------------- *) + +let poly_add = new_recursive_definition list_RECURSION + `([] ++ l2 = l2) /\ + ((CONS h t) ++ l2 = + (if l2 = [] then CONS h t + else CONS (h + HD l2) (t ++ TL l2)))`;; + +let poly_cmul = new_recursive_definition list_RECURSION + `(c ## [] = []) /\ + (c ## (CONS h t) = CONS (c * h) (c ## t))`;; + +let poly_neg = new_definition + `neg = (##) (--(Cx(&1)))`;; + +let poly_mul = new_recursive_definition list_RECURSION + `([] ** l2 = []) /\ + ((CONS h t) ** l2 = + if t = [] then h ## l2 + else (h ## l2) ++ CONS (Cx(&0)) (t ** l2))`;; + +let poly_exp = new_recursive_definition num_RECURSION + `(p exp 0 = [Cx(&1)]) /\ + (p exp (SUC n) = p ** p exp n)`;; + +(* ------------------------------------------------------------------------- *) +(* Useful clausifications. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ADD_CLAUSES = prove + (`([] ++ p2 = p2) /\ + (p1 ++ [] = p1) /\ + ((CONS h1 t1) ++ (CONS h2 t2) = CONS (h1 + h2) (t1 ++ t2))`, + REWRITE_TAC[poly_add; NOT_CONS_NIL; HD; TL] THEN + SPEC_TAC(`p1:complex list`,`p1:complex list`) THEN + LIST_INDUCT_TAC THEN ASM_REWRITE_TAC[poly_add]);; + +let POLY_CMUL_CLAUSES = prove + (`(c ## [] = []) /\ + (c ## (CONS h t) = CONS (c * h) (c ## t))`, + REWRITE_TAC[poly_cmul]);; + +let POLY_NEG_CLAUSES = prove + (`(neg [] = []) /\ + (neg (CONS h t) = CONS (--h) (neg t))`, + REWRITE_TAC[poly_neg; POLY_CMUL_CLAUSES; + COMPLEX_MUL_LNEG; COMPLEX_MUL_LID]);; + +let POLY_MUL_CLAUSES = prove + (`([] ** p2 = []) /\ + ([h1] ** p2 = h1 ## p2) /\ + ((CONS h1 (CONS k1 t1)) ** p2 = + h1 ## p2 ++ CONS (Cx(&0)) (CONS k1 t1 ** p2))`, + REWRITE_TAC[poly_mul; NOT_CONS_NIL]);; + +(* ------------------------------------------------------------------------- *) +(* Various natural consequences of syntactic definitions. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ADD = prove + (`!p1 p2 x. poly (p1 ++ p2) x = poly p1 x + poly p2 x`, + LIST_INDUCT_TAC THEN REWRITE_TAC[poly_add; poly; COMPLEX_ADD_LID] THEN + LIST_INDUCT_TAC THEN + ASM_REWRITE_TAC[NOT_CONS_NIL; HD; TL; poly; COMPLEX_ADD_RID] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_CMUL = prove + (`!p c x. poly (c ## p) x = c * poly p x`, + LIST_INDUCT_TAC THEN ASM_REWRITE_TAC[poly; poly_cmul] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_NEG = prove + (`!p x. poly (neg p) x = --(poly p x)`, + REWRITE_TAC[poly_neg; POLY_CMUL] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_MUL = prove + (`!x p1 p2. poly (p1 ** p2) x = poly p1 x * poly p2 x`, + GEN_TAC THEN LIST_INDUCT_TAC THEN + REWRITE_TAC[poly_mul; poly; COMPLEX_MUL_LZERO; POLY_CMUL; POLY_ADD] THEN + SPEC_TAC(`h:complex`,`h:complex`) THEN + SPEC_TAC(`t:complex list`,`t:complex list`) THEN + LIST_INDUCT_TAC THEN + REWRITE_TAC[poly_mul; POLY_CMUL; POLY_ADD; poly; POLY_CMUL; + COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; NOT_CONS_NIL] THEN + ASM_REWRITE_TAC[POLY_ADD; POLY_CMUL; poly] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_EXP = prove + (`!p n x. poly (p exp n) x = (poly p x) pow n`, + GEN_TAC THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[poly_exp; complex_pow; POLY_MUL] THEN + REWRITE_TAC[poly] THEN SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Lemmas. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ADD_RZERO = prove + (`!p. poly (p ++ []) = poly p`, + REWRITE_TAC[FUN_EQ_THM; POLY_ADD; poly; COMPLEX_ADD_RID]);; + +let POLY_MUL_ASSOC = prove + (`!p q r. poly (p ** (q ** r)) = poly ((p ** q) ** r)`, + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; COMPLEX_MUL_ASSOC]);; + +let POLY_EXP_ADD = prove + (`!d n p. poly(p exp (n + d)) = poly(p exp n ** p exp d)`, + REWRITE_TAC[FUN_EQ_THM; POLY_MUL] THEN + INDUCT_TAC THEN ASM_REWRITE_TAC[POLY_MUL; ADD_CLAUSES; poly_exp; poly] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Key property that f(a) = 0 ==> (x - a) divides p(x). Very delicate! *) +(* ------------------------------------------------------------------------- *) + +let POLY_LINEAR_REM = prove + (`!t h. ?q r. CONS h t = [r] ++ [--a; Cx(&1)] ** q`, + LIST_INDUCT_TAC THEN REWRITE_TAC[] THENL + [GEN_TAC THEN EXISTS_TAC `[]:complex list` THEN + EXISTS_TAC `h:complex` THEN + REWRITE_TAC[poly_add; poly_mul; poly_cmul; NOT_CONS_NIL] THEN + REWRITE_TAC[HD; TL; COMPLEX_ADD_RID]; + X_GEN_TAC `k:complex` THEN + POP_ASSUM(STRIP_ASSUME_TAC o SPEC `h:complex`) THEN + EXISTS_TAC `CONS (r:complex) q` THEN EXISTS_TAC `r * a + k` THEN + ASM_REWRITE_TAC[POLY_ADD_CLAUSES; POLY_MUL_CLAUSES; poly_cmul] THEN + REWRITE_TAC[CONS_11] THEN CONJ_TAC THENL + [SIMPLE_COMPLEX_ARITH_TAC; ALL_TAC] THEN + SPEC_TAC(`q:complex list`,`q:complex list`) THEN + LIST_INDUCT_TAC THEN + REWRITE_TAC[POLY_ADD_CLAUSES; POLY_MUL_CLAUSES; poly_cmul] THEN + REWRITE_TAC[COMPLEX_ADD_RID; COMPLEX_MUL_LID] THEN + REWRITE_TAC[COMPLEX_ADD_AC]]);; + +let POLY_LINEAR_DIVIDES = prove + (`!a p. (poly p a = Cx(&0)) <=> (p = []) \/ ?q. p = [--a; Cx(&1)] ** q`, + GEN_TAC THEN LIST_INDUCT_TAC THENL + [REWRITE_TAC[poly]; ALL_TAC] THEN + EQ_TAC THEN STRIP_TAC THENL + [DISJ2_TAC THEN STRIP_ASSUME_TAC(SPEC_ALL POLY_LINEAR_REM) THEN + EXISTS_TAC `q:complex list` THEN ASM_REWRITE_TAC[] THEN + SUBGOAL_THEN `r = Cx(&0)` SUBST_ALL_TAC THENL + [UNDISCH_TAC `poly (CONS h t) a = Cx(&0)` THEN + ASM_REWRITE_TAC[] THEN REWRITE_TAC[POLY_ADD; POLY_MUL] THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; + COMPLEX_MUL_RID] THEN + REWRITE_TAC[COMPLEX_ADD_LINV] THEN SIMPLE_COMPLEX_ARITH_TAC; + REWRITE_TAC[poly_mul] THEN REWRITE_TAC[NOT_CONS_NIL] THEN + SPEC_TAC(`q:complex list`,`q:complex list`) THEN LIST_INDUCT_TAC THENL + [REWRITE_TAC[poly_cmul; poly_add; NOT_CONS_NIL; + HD; TL; COMPLEX_ADD_LID]; + REWRITE_TAC[poly_cmul; poly_add; NOT_CONS_NIL; + HD; TL; COMPLEX_ADD_LID]]]; + ASM_REWRITE_TAC[] THEN REWRITE_TAC[poly]; + ASM_REWRITE_TAC[] THEN REWRITE_TAC[poly] THEN + REWRITE_TAC[POLY_MUL] THEN REWRITE_TAC[poly] THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; COMPLEX_MUL_RID] THEN + REWRITE_TAC[COMPLEX_ADD_LINV] THEN SIMPLE_COMPLEX_ARITH_TAC]);; + +(* ------------------------------------------------------------------------- *) +(* Thanks to the finesse of the above, we can use length rather than degree. *) +(* ------------------------------------------------------------------------- *) + +let POLY_LENGTH_MUL = prove + (`!q. LENGTH([--a; Cx(&1)] ** q) = SUC(LENGTH q)`, + let lemma = prove + (`!p h k a. LENGTH (k ## p ++ CONS h (a ## p)) = SUC(LENGTH p)`, + LIST_INDUCT_TAC THEN + ASM_REWRITE_TAC[poly_cmul; POLY_ADD_CLAUSES; LENGTH]) in + REWRITE_TAC[poly_mul; NOT_CONS_NIL; lemma]);; + +(* ------------------------------------------------------------------------- *) +(* Thus a nontrivial polynomial of degree n has no more than n roots. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ROOTS_INDEX_LEMMA = prove + (`!n. !p. ~(poly p = poly []) /\ (LENGTH p = n) + ==> ?i. !x. (poly p x = Cx(&0)) ==> ?m. m <= n /\ (x = i m)`, + INDUCT_TAC THENL + [REWRITE_TAC[LENGTH_EQ_NIL] THEN MESON_TAC[]; + REPEAT STRIP_TAC THEN ASM_CASES_TAC `?a. poly p a = Cx(&0)` THENL + [UNDISCH_TAC `?a. poly p a = Cx(&0)` THEN + DISCH_THEN(CHOOSE_THEN MP_TAC) THEN + GEN_REWRITE_TAC LAND_CONV [POLY_LINEAR_DIVIDES] THEN + DISCH_THEN(DISJ_CASES_THEN MP_TAC) THENL [ASM_MESON_TAC[]; ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `q:complex list` SUBST_ALL_TAC) THEN + FIRST_ASSUM(UNDISCH_TAC o check is_forall o concl) THEN + UNDISCH_TAC `~(poly ([-- a; Cx(&1)] ** q) = poly [])` THEN + POP_ASSUM MP_TAC THEN REWRITE_TAC[POLY_LENGTH_MUL; SUC_INJ] THEN + DISCH_TAC THEN ASM_CASES_TAC `poly q = poly []` THENL + [ASM_REWRITE_TAC[POLY_MUL; poly; COMPLEX_MUL_RZERO; FUN_EQ_THM]; + DISCH_THEN(K ALL_TAC)] THEN + DISCH_THEN(MP_TAC o SPEC `q:complex list`) THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_TAC `i:num->complex`) THEN + EXISTS_TAC `\m. if m = SUC n then a:complex else i m` THEN + REWRITE_TAC[POLY_MUL; LE; COMPLEX_ENTIRE] THEN + X_GEN_TAC `x:complex` THEN DISCH_THEN(DISJ_CASES_THEN MP_TAC) THENL + [DISCH_THEN(fun th -> EXISTS_TAC `SUC n` THEN MP_TAC th) THEN + REWRITE_TAC[poly] THEN SIMPLE_COMPLEX_ARITH_TAC; + DISCH_THEN(ANTE_RES_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `m:num` STRIP_ASSUME_TAC) THEN + EXISTS_TAC `m:num` THEN ASM_REWRITE_TAC[] THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[] THEN + UNDISCH_TAC `m:num <= n` THEN ASM_REWRITE_TAC[] THEN ARITH_TAC]; + UNDISCH_TAC `~(?a. poly p a = Cx(&0))` THEN + REWRITE_TAC[NOT_EXISTS_THM] THEN DISCH_TAC THEN ASM_REWRITE_TAC[]]]);; + +let POLY_ROOTS_INDEX_LENGTH = prove + (`!p. ~(poly p = poly []) + ==> ?i. !x. (poly p(x) = Cx(&0)) ==> ?n. n <= LENGTH p /\ (x = i n)`, + MESON_TAC[POLY_ROOTS_INDEX_LEMMA]);; + +let POLY_ROOTS_FINITE_LEMMA = prove + (`!p. ~(poly p = poly []) + ==> ?N i. !x. (poly p(x) = Cx(&0)) ==> ?n:num. n < N /\ (x = i n)`, + MESON_TAC[POLY_ROOTS_INDEX_LENGTH; LT_SUC_LE]);; + +let FINITE_LEMMA = prove + (`!i N P. (!x. P x ==> ?n:num. n < N /\ (x = i n)) + ==> ?a. !x. P x ==> norm(x) < a`, + GEN_TAC THEN ONCE_REWRITE_TAC[RIGHT_IMP_EXISTS_THM] THEN INDUCT_TAC THENL + [REWRITE_TAC[LT] THEN MESON_TAC[]; ALL_TAC] THEN + X_GEN_TAC `P:complex->bool` THEN + POP_ASSUM(MP_TAC o SPEC `\z. P z /\ ~(z = (i:num->complex) N)`) THEN + DISCH_THEN(X_CHOOSE_TAC `a:real`) THEN + EXISTS_TAC `abs(a) + norm(i(N:num)) + &1` THEN + POP_ASSUM MP_TAC THEN REWRITE_TAC[LT] THEN + SUBGOAL_THEN `(!x. norm(x) < abs(a) + norm(x) + &1) /\ + (!x y. norm(x) < a ==> norm(x) < abs(a) + norm(y) + &1)` + (fun th -> MP_TAC th THEN MESON_TAC[]) THEN + CONJ_TAC THENL [REAL_ARITH_TAC; ALL_TAC] THEN + REPEAT GEN_TAC THEN MP_TAC(SPEC `y:complex` COMPLEX_NORM_POS) THEN + REAL_ARITH_TAC);; + +let POLY_ROOTS_FINITE = prove + (`!p. ~(poly p = poly []) <=> + ?N i. !x. (poly p(x) = Cx(&0)) ==> ?n:num. n < N /\ (x = i n)`, + GEN_TAC THEN EQ_TAC THEN REWRITE_TAC[POLY_ROOTS_FINITE_LEMMA] THEN + REWRITE_TAC[FUN_EQ_THM; LEFT_IMP_EXISTS_THM; NOT_FORALL_THM; poly] THEN + MP_TAC(GENL [`i:num->complex`; `N:num`] + (SPECL [`i:num->complex`; `N:num`; `\x. poly p x = Cx(&0)`] + FINITE_LEMMA)) THEN + REWRITE_TAC[] THEN MESON_TAC[REAL_ARITH `~(abs(x) < x)`; COMPLEX_NORM_CX]);; + +(* ------------------------------------------------------------------------- *) +(* Hence get entirety and cancellation for polynomials. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ENTIRE_LEMMA = prove + (`!p q. ~(poly p = poly []) /\ ~(poly q = poly []) + ==> ~(poly (p ** q) = poly [])`, + REPEAT GEN_TAC THEN REWRITE_TAC[POLY_ROOTS_FINITE] THEN + DISCH_THEN(CONJUNCTS_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `N2:num` (X_CHOOSE_TAC `i2:num->complex`)) THEN + DISCH_THEN(X_CHOOSE_THEN `N1:num` (X_CHOOSE_TAC `i1:num->complex`)) THEN + EXISTS_TAC `N1 + N2:num` THEN + EXISTS_TAC `\n:num. if n < N1 then i1(n):complex else i2(n - N1)` THEN + X_GEN_TAC `x:complex` THEN REWRITE_TAC[COMPLEX_ENTIRE; POLY_MUL] THEN + DISCH_THEN(DISJ_CASES_THEN (ANTE_RES_THEN (X_CHOOSE_TAC `n:num`))) THENL + [EXISTS_TAC `n:num` THEN ASM_REWRITE_TAC[] THEN + FIRST_ASSUM(MP_TAC o CONJUNCT1) THEN ARITH_TAC; + EXISTS_TAC `N1 + n:num` THEN ASM_REWRITE_TAC[LT_ADD_LCANCEL] THEN + REWRITE_TAC[ARITH_RULE `~(m + n < m:num)`] THEN + AP_TERM_TAC THEN ARITH_TAC]);; + +let POLY_ENTIRE = prove + (`!p q. (poly (p ** q) = poly []) <=> + (poly p = poly []) \/ (poly q = poly [])`, + REPEAT GEN_TAC THEN EQ_TAC THENL + [MESON_TAC[POLY_ENTIRE_LEMMA]; + REWRITE_TAC[FUN_EQ_THM; POLY_MUL] THEN + STRIP_TAC THEN + ASM_REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_MUL_LZERO; poly]]);; + +let POLY_MUL_LCANCEL = prove + (`!p q r. (poly (p ** q) = poly (p ** r)) <=> + (poly p = poly []) \/ (poly q = poly r)`, + let lemma1 = prove + (`!p q. (poly (p ++ neg q) = poly []) <=> (poly p = poly q)`, + REWRITE_TAC[FUN_EQ_THM; POLY_ADD; POLY_NEG; poly] THEN + REWRITE_TAC[SIMPLE_COMPLEX_ARITH `(p + --q = Cx(&0)) <=> (p = q)`]) in + let lemma2 = prove + (`!p q r. poly (p ** q ++ neg(p ** r)) = poly (p ** (q ++ neg(r)))`, + REWRITE_TAC[FUN_EQ_THM; POLY_ADD; POLY_NEG; POLY_MUL] THEN + SIMPLE_COMPLEX_ARITH_TAC) in + ONCE_REWRITE_TAC[GSYM lemma1] THEN + REWRITE_TAC[lemma2; POLY_ENTIRE] THEN + REWRITE_TAC[lemma1]);; + +let POLY_EXP_EQ_0 = prove + (`!p n. (poly (p exp n) = poly []) <=> (poly p = poly []) /\ ~(n = 0)`, + REPEAT GEN_TAC THEN REWRITE_TAC[FUN_EQ_THM; poly] THEN + REWRITE_TAC[LEFT_AND_FORALL_THM] THEN AP_TERM_TAC THEN ABS_TAC THEN + SPEC_TAC(`n:num`,`n:num`) THEN INDUCT_TAC THEN + REWRITE_TAC[poly_exp; poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; + CX_INJ; REAL_OF_NUM_EQ; ARITH; NOT_SUC] THEN + ASM_REWRITE_TAC[POLY_MUL; poly; COMPLEX_ENTIRE] THEN + CONV_TAC TAUT);; + +let POLY_PRIME_EQ_0 = prove + (`!a. ~(poly [a ; Cx(&1)] = poly [])`, + GEN_TAC THEN REWRITE_TAC[FUN_EQ_THM; poly] THEN + DISCH_THEN(MP_TAC o SPEC `Cx(&1) - a`) THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_EXP_PRIME_EQ_0 = prove + (`!a n. ~(poly ([a ; Cx(&1)] exp n) = poly [])`, + MESON_TAC[POLY_EXP_EQ_0; POLY_PRIME_EQ_0]);; + +(* ------------------------------------------------------------------------- *) +(* Can also prove a more "constructive" notion of polynomial being trivial. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ZERO_LEMMA = prove + (`!h t. (poly (CONS h t) = poly []) ==> (h = Cx(&0)) /\ (poly t = poly [])`, + let lemma = REWRITE_RULE[FUN_EQ_THM; poly] POLY_ROOTS_FINITE in + REPEAT GEN_TAC THEN REWRITE_TAC[FUN_EQ_THM; poly] THEN + ASM_CASES_TAC `h = Cx(&0)` THEN ASM_REWRITE_TAC[] THENL + [REWRITE_TAC[COMPLEX_ADD_LID]; + DISCH_THEN(MP_TAC o SPEC `Cx(&0)`) THEN + POP_ASSUM MP_TAC THEN SIMPLE_COMPLEX_ARITH_TAC] THEN + CONV_TAC CONTRAPOS_CONV THEN + DISCH_THEN(MP_TAC o REWRITE_RULE[lemma]) THEN + DISCH_THEN(X_CHOOSE_THEN `N:num` (X_CHOOSE_TAC `i:num->complex`)) THEN + MP_TAC(SPECL + [`i:num->complex`; `N:num`; `\x. poly t x = Cx(&0)`] FINITE_LEMMA) THEN + ASM_REWRITE_TAC[] THEN DISCH_THEN(X_CHOOSE_TAC `a:real`) THEN + DISCH_THEN(MP_TAC o SPEC `Cx(abs(a) + &1)`) THEN + REWRITE_TAC[COMPLEX_ENTIRE; DE_MORGAN_THM] THEN CONJ_TAC THENL + [REWRITE_TAC[CX_INJ] THEN REAL_ARITH_TAC; + DISCH_THEN(MP_TAC o MATCH_MP + (ASSUME `!x. (poly t x = Cx(&0)) ==> norm(x) < a`)) THEN + REWRITE_TAC[COMPLEX_NORM_CX] THEN REAL_ARITH_TAC]);; + +let POLY_ZERO = prove + (`!p. (poly p = poly []) <=> ALL (\c. c = Cx(&0)) p`, + LIST_INDUCT_TAC THEN ASM_REWRITE_TAC[ALL] THEN EQ_TAC THENL + [DISCH_THEN(MP_TAC o MATCH_MP POLY_ZERO_LEMMA) THEN ASM_REWRITE_TAC[]; + POP_ASSUM(SUBST1_TAC o SYM) THEN STRIP_TAC THEN + ASM_REWRITE_TAC[FUN_EQ_THM; poly] THEN SIMPLE_COMPLEX_ARITH_TAC]);; + +(* ------------------------------------------------------------------------- *) +(* Basics of divisibility. *) +(* ------------------------------------------------------------------------- *) + +let divides = new_definition + `p1 divides p2 <=> ?q. poly p2 = poly (p1 ** q)`;; + +let POLY_PRIMES = prove + (`!a p q. [a; Cx(&1)] divides (p ** q) <=> + [a; Cx(&1)] divides p \/ [a; Cx(&1)] divides q`, + REPEAT GEN_TAC THEN REWRITE_TAC[divides; POLY_MUL; FUN_EQ_THM; poly] THEN + REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; COMPLEX_MUL_RID] THEN + EQ_TAC THENL + [DISCH_THEN(X_CHOOSE_THEN `r:complex list` (MP_TAC o SPEC `--a`)) THEN + REWRITE_TAC[COMPLEX_ENTIRE; GSYM complex_sub; + COMPLEX_SUB_REFL; COMPLEX_MUL_LZERO] THEN + DISCH_THEN DISJ_CASES_TAC THENL [DISJ1_TAC; DISJ2_TAC] THEN + (POP_ASSUM(MP_TAC o REWRITE_RULE[POLY_LINEAR_DIVIDES]) THEN + REWRITE_TAC[COMPLEX_NEG_NEG] THEN + DISCH_THEN(DISJ_CASES_THEN2 SUBST_ALL_TAC + (X_CHOOSE_THEN `s:complex list` SUBST_ALL_TAC)) THENL + [EXISTS_TAC `[]:complex list` THEN REWRITE_TAC[poly; COMPLEX_MUL_RZERO]; + EXISTS_TAC `s:complex list` THEN GEN_TAC THEN + REWRITE_TAC[POLY_MUL; poly] THEN SIMPLE_COMPLEX_ARITH_TAC]); + DISCH_THEN(DISJ_CASES_THEN(X_CHOOSE_TAC `s:complex list`)) THEN + ASM_REWRITE_TAC[] THENL + [EXISTS_TAC `s ** q`; EXISTS_TAC `p ** s`] THEN + GEN_TAC THEN REWRITE_TAC[POLY_MUL] THEN SIMPLE_COMPLEX_ARITH_TAC]);; + +let POLY_DIVIDES_REFL = prove + (`!p. p divides p`, + GEN_TAC THEN REWRITE_TAC[divides] THEN EXISTS_TAC `[Cx(&1)]` THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly] THEN SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_DIVIDES_TRANS = prove + (`!p q r. p divides q /\ q divides r ==> p divides r`, + REPEAT GEN_TAC THEN REWRITE_TAC[divides] THEN + DISCH_THEN(CONJUNCTS_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `s:complex list` ASSUME_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `t:complex list` ASSUME_TAC) THEN + EXISTS_TAC `t ** s` THEN + ASM_REWRITE_TAC[FUN_EQ_THM; POLY_MUL; COMPLEX_MUL_ASSOC]);; + +let POLY_DIVIDES_EXP = prove + (`!p m n. m <= n ==> (p exp m) divides (p exp n)`, + REPEAT GEN_TAC THEN REWRITE_TAC[LE_EXISTS] THEN + DISCH_THEN(X_CHOOSE_THEN `d:num` SUBST1_TAC) THEN + SPEC_TAC(`d:num`,`d:num`) THEN INDUCT_TAC THEN + REWRITE_TAC[ADD_CLAUSES; POLY_DIVIDES_REFL] THEN + MATCH_MP_TAC POLY_DIVIDES_TRANS THEN + EXISTS_TAC `p exp (m + d)` THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[divides] THEN EXISTS_TAC `p:complex list` THEN + REWRITE_TAC[poly_exp; FUN_EQ_THM; POLY_MUL] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_EXP_DIVIDES = prove + (`!p q m n. (p exp n) divides q /\ m <= n ==> (p exp m) divides q`, + MESON_TAC[POLY_DIVIDES_TRANS; POLY_DIVIDES_EXP]);; + +let POLY_DIVIDES_ADD = prove + (`!p q r. p divides q /\ p divides r ==> p divides (q ++ r)`, + REPEAT GEN_TAC THEN REWRITE_TAC[divides] THEN + DISCH_THEN(CONJUNCTS_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `s:complex list` ASSUME_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `t:complex list` ASSUME_TAC) THEN + EXISTS_TAC `t ++ s` THEN + ASM_REWRITE_TAC[FUN_EQ_THM; POLY_ADD; POLY_MUL] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_DIVIDES_SUB = prove + (`!p q r. p divides q /\ p divides (q ++ r) ==> p divides r`, + REPEAT GEN_TAC THEN REWRITE_TAC[divides] THEN + DISCH_THEN(CONJUNCTS_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `s:complex list` ASSUME_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `t:complex list` ASSUME_TAC) THEN + EXISTS_TAC `s ++ neg(t)` THEN + POP_ASSUM_LIST(MP_TAC o end_itlist CONJ) THEN + REWRITE_TAC[FUN_EQ_THM; POLY_ADD; POLY_MUL; POLY_NEG] THEN + DISCH_THEN(STRIP_ASSUME_TAC o GSYM) THEN + REWRITE_TAC[COMPLEX_ADD_LDISTRIB; COMPLEX_MUL_RNEG] THEN + ASM_REWRITE_TAC[] THEN SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_DIVIDES_SUB2 = prove + (`!p q r. p divides r /\ p divides (q ++ r) ==> p divides q`, + REPEAT STRIP_TAC THEN MATCH_MP_TAC POLY_DIVIDES_SUB THEN + EXISTS_TAC `r:complex list` THEN ASM_REWRITE_TAC[] THEN + UNDISCH_TAC `p divides (q ++ r)` THEN + REWRITE_TAC[divides; POLY_ADD; FUN_EQ_THM; POLY_MUL] THEN + DISCH_THEN(X_CHOOSE_TAC `s:complex list`) THEN + EXISTS_TAC `s:complex list` THEN + X_GEN_TAC `x:complex` THEN POP_ASSUM(MP_TAC o SPEC `x:complex`) THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_DIVIDES_ZERO = prove + (`!p q. (poly p = poly []) ==> q divides p`, + REPEAT GEN_TAC THEN DISCH_TAC THEN REWRITE_TAC[divides] THEN + EXISTS_TAC `[]:complex list` THEN + ASM_REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly; COMPLEX_MUL_RZERO]);; + +(* ------------------------------------------------------------------------- *) +(* At last, we can consider the order of a root. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ORDER_EXISTS = prove + (`!a d. !p. (LENGTH p = d) /\ ~(poly p = poly []) + ==> ?n. ([--a; Cx(&1)] exp n) divides p /\ + ~(([--a; Cx(&1)] exp (SUC n)) divides p)`, + GEN_TAC THEN + (STRIP_ASSUME_TAC o prove_recursive_functions_exist num_RECURSION) + `(!p q. mulexp 0 p q = q) /\ + (!p q n. mulexp (SUC n) p q = p ** (mulexp n p q))` THEN + SUBGOAL_THEN + `!d. !p. (LENGTH p = d) /\ ~(poly p = poly []) + ==> ?n q. (p = mulexp (n:num) [--a; Cx(&1)] q) /\ + ~(poly q a = Cx(&0))` + MP_TAC THENL + [INDUCT_TAC THENL + [REWRITE_TAC[LENGTH_EQ_NIL] THEN MESON_TAC[]; ALL_TAC] THEN + X_GEN_TAC `p:complex list` THEN + ASM_CASES_TAC `poly p a = Cx(&0)` THENL + [STRIP_TAC THEN UNDISCH_TAC `poly p a = Cx(&0)` THEN + DISCH_THEN(MP_TAC o REWRITE_RULE[POLY_LINEAR_DIVIDES]) THEN + DISCH_THEN(DISJ_CASES_THEN MP_TAC) THENL [ASM_MESON_TAC[]; ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `q:complex list` SUBST_ALL_TAC) THEN + UNDISCH_TAC + `!p. (LENGTH p = d) /\ ~(poly p = poly []) + ==> ?n q. (p = mulexp (n:num) [--a; Cx(&1)] q) /\ + ~(poly q a = Cx(&0))` THEN + DISCH_THEN(MP_TAC o SPEC `q:complex list`) THEN + RULE_ASSUM_TAC(REWRITE_RULE[POLY_LENGTH_MUL; POLY_ENTIRE; + DE_MORGAN_THM; SUC_INJ]) THEN + ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_THEN `n:num` + (X_CHOOSE_THEN `s:complex list` STRIP_ASSUME_TAC)) THEN + EXISTS_TAC `SUC n` THEN EXISTS_TAC `s:complex list` THEN + ASM_REWRITE_TAC[]; + STRIP_TAC THEN EXISTS_TAC `0` THEN EXISTS_TAC `p:complex list` THEN + ASM_REWRITE_TAC[]]; + DISCH_TAC THEN REPEAT GEN_TAC THEN + DISCH_THEN(ANTE_RES_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `n:num` + (X_CHOOSE_THEN `s:complex list` STRIP_ASSUME_TAC)) THEN + EXISTS_TAC `n:num` THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[divides] THEN CONJ_TAC THENL + [EXISTS_TAC `s:complex list` THEN + SPEC_TAC(`n:num`,`n:num`) THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[poly_exp; FUN_EQ_THM; POLY_MUL; poly] THEN + SIMPLE_COMPLEX_ARITH_TAC; + DISCH_THEN(X_CHOOSE_THEN `r:complex list` MP_TAC) THEN + SPEC_TAC(`n:num`,`n:num`) THEN INDUCT_TAC THEN + ASM_REWRITE_TAC[] THENL + [UNDISCH_TAC `~(poly s a = Cx(&0))` THEN CONV_TAC CONTRAPOS_CONV THEN + REWRITE_TAC[] THEN DISCH_THEN SUBST1_TAC THEN + REWRITE_TAC[poly; poly_exp; POLY_MUL] THEN SIMPLE_COMPLEX_ARITH_TAC; + REWRITE_TAC[] THEN ONCE_ASM_REWRITE_TAC[] THEN + ONCE_REWRITE_TAC[poly_exp] THEN + REWRITE_TAC[GSYM POLY_MUL_ASSOC; POLY_MUL_LCANCEL] THEN + REWRITE_TAC[DE_MORGAN_THM] THEN CONJ_TAC THENL + [REWRITE_TAC[FUN_EQ_THM] THEN + DISCH_THEN(MP_TAC o SPEC `a + Cx(&1)`) THEN + REWRITE_TAC[poly] THEN SIMPLE_COMPLEX_ARITH_TAC; + DISCH_THEN(ANTE_RES_THEN MP_TAC) THEN REWRITE_TAC[]]]]]);; + +let POLY_ORDER = prove + (`!p a. ~(poly p = poly []) + ==> ?!n. ([--a; Cx(&1)] exp n) divides p /\ + ~(([--a; Cx(&1)] exp (SUC n)) divides p)`, + MESON_TAC[POLY_ORDER_EXISTS; POLY_EXP_DIVIDES; LE_SUC_LT; LT_CASES]);; + +(* ------------------------------------------------------------------------- *) +(* Definition of order. *) +(* ------------------------------------------------------------------------- *) + +let order = new_definition + `order a p = @n. ([--a; Cx(&1)] exp n) divides p /\ + ~(([--a; Cx(&1)] exp (SUC n)) divides p)`;; + +let ORDER = prove + (`!p a n. ([--a; Cx(&1)] exp n) divides p /\ + ~(([--a; Cx(&1)] exp (SUC n)) divides p) <=> + (n = order a p) /\ + ~(poly p = poly [])`, + REPEAT GEN_TAC THEN REWRITE_TAC[order] THEN + EQ_TAC THEN STRIP_TAC THENL + [SUBGOAL_THEN `~(poly p = poly [])` ASSUME_TAC THENL + [FIRST_ASSUM(UNDISCH_TAC o check is_neg o concl) THEN + CONV_TAC CONTRAPOS_CONV THEN REWRITE_TAC[divides] THEN + DISCH_THEN SUBST1_TAC THEN EXISTS_TAC `[]:complex list` THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly; COMPLEX_MUL_RZERO]; + ASM_REWRITE_TAC[] THEN CONV_TAC SYM_CONV THEN + MATCH_MP_TAC SELECT_UNIQUE THEN REWRITE_TAC[]]; + ONCE_ASM_REWRITE_TAC[] THEN CONV_TAC SELECT_CONV] THEN + ASM_MESON_TAC[POLY_ORDER]);; + +let ORDER_THM = prove + (`!p a. ~(poly p = poly []) + ==> ([--a; Cx(&1)] exp (order a p)) divides p /\ + ~(([--a; Cx(&1)] exp (SUC(order a p))) divides p)`, + MESON_TAC[ORDER]);; + +let ORDER_UNIQUE = prove + (`!p a n. ~(poly p = poly []) /\ + ([--a; Cx(&1)] exp n) divides p /\ + ~(([--a; Cx(&1)] exp (SUC n)) divides p) + ==> (n = order a p)`, + MESON_TAC[ORDER]);; + +let ORDER_POLY = prove + (`!p q a. (poly p = poly q) ==> (order a p = order a q)`, + REPEAT STRIP_TAC THEN + ASM_REWRITE_TAC[order; divides; FUN_EQ_THM; POLY_MUL]);; + +let ORDER_ROOT = prove + (`!p a. (poly p a = Cx(&0)) <=> (poly p = poly []) \/ ~(order a p = 0)`, + REPEAT GEN_TAC THEN ASM_CASES_TAC `poly p = poly []` THEN + ASM_REWRITE_TAC[poly] THEN EQ_TAC THENL + [DISCH_THEN(MP_TAC o REWRITE_RULE[POLY_LINEAR_DIVIDES]) THEN + ASM_CASES_TAC `p:complex list = []` THENL [ASM_MESON_TAC[]; ALL_TAC] THEN + ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_THEN `q:complex list` SUBST_ALL_TAC) THEN DISCH_TAC THEN + FIRST_ASSUM(MP_TAC o SPEC `a:complex` o MATCH_MP ORDER_THM) THEN + ASM_REWRITE_TAC[poly_exp; DE_MORGAN_THM] THEN DISJ2_TAC THEN + REWRITE_TAC[divides] THEN EXISTS_TAC `q:complex list` THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly] THEN SIMPLE_COMPLEX_ARITH_TAC; + DISCH_TAC THEN + FIRST_ASSUM(MP_TAC o SPEC `a:complex` o MATCH_MP ORDER_THM) THEN + UNDISCH_TAC `~(order a p = 0)` THEN + SPEC_TAC(`order a p`,`n:num`) THEN + INDUCT_TAC THEN ASM_REWRITE_TAC[poly_exp; NOT_SUC] THEN + DISCH_THEN(MP_TAC o CONJUNCT1) THEN REWRITE_TAC[divides] THEN + DISCH_THEN(X_CHOOSE_THEN `s:complex list` SUBST1_TAC) THEN + REWRITE_TAC[POLY_MUL; poly] THEN SIMPLE_COMPLEX_ARITH_TAC]);; + +let ORDER_DIVIDES = prove + (`!p a n. ([--a; Cx(&1)] exp n) divides p <=> + (poly p = poly []) \/ n <= order a p`, + REPEAT GEN_TAC THEN ASM_CASES_TAC `poly p = poly []` THEN + ASM_REWRITE_TAC[] THENL + [ASM_REWRITE_TAC[divides] THEN EXISTS_TAC `[]:complex list` THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly; COMPLEX_MUL_RZERO]; + ASM_MESON_TAC[ORDER_THM; POLY_EXP_DIVIDES; NOT_LE; LE_SUC_LT]]);; + +let ORDER_DECOMP = prove + (`!p a. ~(poly p = poly []) + ==> ?q. (poly p = poly (([--a; Cx(&1)] exp (order a p)) ** q)) /\ + ~([--a; Cx(&1)] divides q)`, + REPEAT STRIP_TAC THEN FIRST_ASSUM(MP_TAC o MATCH_MP ORDER_THM) THEN + DISCH_THEN(CONJUNCTS_THEN2 MP_TAC ASSUME_TAC o SPEC `a:complex`) THEN + DISCH_THEN(X_CHOOSE_TAC `q:complex list` o REWRITE_RULE[divides]) THEN + EXISTS_TAC `q:complex list` THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_TAC `r: complex list` o REWRITE_RULE[divides]) THEN + UNDISCH_TAC `~([-- a; Cx(&1)] exp SUC (order a p) divides p)` THEN + ASM_REWRITE_TAC[] THEN REWRITE_TAC[divides] THEN + EXISTS_TAC `r:complex list` THEN + ASM_REWRITE_TAC[POLY_MUL; FUN_EQ_THM; poly_exp] THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Important composition properties of orders. *) +(* ------------------------------------------------------------------------- *) + +let ORDER_MUL = prove + (`!a p q. ~(poly (p ** q) = poly []) ==> + (order a (p ** q) = order a p + order a q)`, + REPEAT GEN_TAC THEN + DISCH_THEN(fun th -> ASSUME_TAC th THEN MP_TAC th) THEN + REWRITE_TAC[POLY_ENTIRE; DE_MORGAN_THM] THEN STRIP_TAC THEN + SUBGOAL_THEN `(order a p + order a q = order a (p ** q)) /\ + ~(poly (p ** q) = poly [])` + MP_TAC THENL [ALL_TAC; MESON_TAC[]] THEN + REWRITE_TAC[GSYM ORDER] THEN CONJ_TAC THENL + [MP_TAC(CONJUNCT1 (SPEC `a:complex` + (MATCH_MP ORDER_THM (ASSUME `~(poly p = poly [])`)))) THEN + DISCH_THEN(X_CHOOSE_TAC `r: complex list` o REWRITE_RULE[divides]) THEN + MP_TAC(CONJUNCT1 (SPEC `a:complex` + (MATCH_MP ORDER_THM (ASSUME `~(poly q = poly [])`)))) THEN + DISCH_THEN(X_CHOOSE_TAC `s: complex list` o REWRITE_RULE[divides]) THEN + REWRITE_TAC[divides; FUN_EQ_THM] THEN EXISTS_TAC `s ** r` THEN + ASM_REWRITE_TAC[POLY_MUL; POLY_EXP_ADD] THEN SIMPLE_COMPLEX_ARITH_TAC; + X_CHOOSE_THEN `r: complex list` STRIP_ASSUME_TAC + (SPEC `a:complex` (MATCH_MP ORDER_DECOMP + (ASSUME `~(poly p = poly [])`))) THEN + X_CHOOSE_THEN `s: complex list` STRIP_ASSUME_TAC + (SPEC `a:complex` (MATCH_MP ORDER_DECOMP + (ASSUME `~(poly q = poly [])`))) THEN + ASM_REWRITE_TAC[divides; FUN_EQ_THM; POLY_EXP_ADD; POLY_MUL; poly_exp] THEN + DISCH_THEN(X_CHOOSE_THEN `t:complex list` STRIP_ASSUME_TAC) THEN + SUBGOAL_THEN `[--a; Cx(&1)] divides (r ** s)` MP_TAC THENL + [ALL_TAC; ASM_REWRITE_TAC[POLY_PRIMES]] THEN + REWRITE_TAC[divides] THEN EXISTS_TAC `t:complex list` THEN + SUBGOAL_THEN `poly ([-- a; Cx(&1)] exp (order a p) ** r ** s) = + poly ([-- a; Cx(&1)] exp (order a p) ** + ([-- a; Cx(&1)] ** t))` + MP_TAC THENL + [ALL_TAC; MESON_TAC[POLY_MUL_LCANCEL; POLY_EXP_PRIME_EQ_0]] THEN + SUBGOAL_THEN `poly ([-- a; Cx(&1)] exp (order a q) ** + [-- a; Cx(&1)] exp (order a p) ** r ** s) = + poly ([-- a; Cx(&1)] exp (order a q) ** + [-- a; Cx(&1)] exp (order a p) ** + [-- a; Cx(&1)] ** t)` + MP_TAC THENL + [ALL_TAC; MESON_TAC[POLY_MUL_LCANCEL; POLY_EXP_PRIME_EQ_0]] THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; POLY_ADD] THEN + FIRST_ASSUM(UNDISCH_TAC o check is_forall o concl) THEN + REWRITE_TAC[COMPLEX_MUL_AC]]);; + +(* ------------------------------------------------------------------------- *) +(* Normalization of a polynomial. *) +(* ------------------------------------------------------------------------- *) + +let normalize = new_recursive_definition list_RECURSION + `(normalize [] = []) /\ + (normalize (CONS h t) = + if normalize t = [] then if h = Cx(&0) then [] else [h] + else CONS h (normalize t))`;; + +let POLY_NORMALIZE = prove + (`!p. poly (normalize p) = poly p`, + LIST_INDUCT_TAC THEN REWRITE_TAC[normalize; poly] THEN + ASM_CASES_TAC `h = Cx(&0)` THEN ASM_REWRITE_TAC[] THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[poly; FUN_EQ_THM] THEN + UNDISCH_TAC `poly (normalize t) = poly t` THEN + DISCH_THEN(SUBST1_TAC o SYM) THEN ASM_REWRITE_TAC[poly] THEN + REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_ADD_LID]);; + +let LENGTH_NORMALIZE_LE = prove + (`!p. LENGTH(normalize p) <= LENGTH p`, + LIST_INDUCT_TAC THEN REWRITE_TAC[LENGTH; normalize; LE_REFL] THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[LENGTH; LE_SUC] THEN + COND_CASES_TAC THEN REWRITE_TAC[LENGTH] THEN ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* The degree of a polynomial. *) +(* ------------------------------------------------------------------------- *) + +let degree = new_definition + `degree p = PRE(LENGTH(normalize p))`;; + +let DEGREE_ZERO = prove + (`!p. (poly p = poly []) ==> (degree p = 0)`, + REPEAT STRIP_TAC THEN REWRITE_TAC[degree] THEN + SUBGOAL_THEN `normalize p = []` SUBST1_TAC THENL + [POP_ASSUM MP_TAC THEN SPEC_TAC(`p:complex list`,`p:complex list`) THEN + REWRITE_TAC[POLY_ZERO] THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[normalize; ALL] THEN + STRIP_TAC THEN ASM_REWRITE_TAC[] THEN + SUBGOAL_THEN `normalize t = []` (fun th -> REWRITE_TAC[th]) THEN + FIRST_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[]; + REWRITE_TAC[LENGTH; PRE]]);; + +(* ------------------------------------------------------------------------- *) +(* Show that the degree is welldefined. *) +(* ------------------------------------------------------------------------- *) + +let POLY_CONS_EQ = prove + (`(poly (CONS h1 t1) = poly (CONS h2 t2)) <=> + (h1 = h2) /\ (poly t1 = poly t2)`, + REWRITE_TAC[FUN_EQ_THM] THEN EQ_TAC THENL [ALL_TAC; SIMP_TAC[poly]] THEN + ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH `(a = b) <=> (a + --b = Cx(&0))`] THEN + REWRITE_TAC[GSYM POLY_NEG; GSYM POLY_ADD] THEN DISCH_TAC THEN + SUBGOAL_THEN `poly (CONS h1 t1 ++ neg(CONS h2 t2)) = poly []` MP_TAC THENL + [ASM_REWRITE_TAC[poly; FUN_EQ_THM]; ALL_TAC] THEN + REWRITE_TAC[poly_neg; poly_cmul; poly_add; NOT_CONS_NIL; HD; TL] THEN + DISCH_THEN(MP_TAC o MATCH_MP POLY_ZERO_LEMMA) THEN + SIMP_TAC[poly; COMPLEX_MUL_LNEG; COMPLEX_MUL_LID]);; + +let POLY_NORMALIZE_ZERO = prove + (`!p. (poly p = poly []) <=> (normalize p = [])`, + REWRITE_TAC[POLY_ZERO] THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[ALL; normalize] THEN + ASM_CASES_TAC `normalize t = []` THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[NOT_CONS_NIL] THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[NOT_CONS_NIL]);; + +let POLY_NORMALIZE_EQ_LEMMA = prove + (`!p q. (poly p = poly q) ==> (normalize p = normalize q)`, + LIST_INDUCT_TAC THENL + [MESON_TAC[POLY_NORMALIZE_ZERO]; ALL_TAC] THEN + LIST_INDUCT_TAC THENL + [MESON_TAC[POLY_NORMALIZE_ZERO]; ALL_TAC] THEN + REWRITE_TAC[POLY_CONS_EQ] THEN + STRIP_TAC THEN ASM_REWRITE_TAC[normalize] THEN + FIRST_X_ASSUM(MP_TAC o SPEC `t':complex list`) THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN SUBST1_TAC THEN REFL_TAC);; + +let POLY_NORMALIZE_EQ = prove + (`!p q. (poly p = poly q) <=> (normalize p = normalize q)`, + MESON_TAC[POLY_NORMALIZE_EQ_LEMMA; POLY_NORMALIZE]);; + +let DEGREE_WELLDEF = prove + (`!p q. (poly p = poly q) ==> (degree p = degree q)`, + SIMP_TAC[degree; POLY_NORMALIZE_EQ]);; + +(* ------------------------------------------------------------------------- *) +(* Degree of a product with a power of linear terms. *) +(* ------------------------------------------------------------------------- *) + +let NORMALIZE_EQ = prove + (`!p. ~(LAST p = Cx(&0)) ==> (normalize p = p)`, + MATCH_MP_TAC list_INDUCT THEN REWRITE_TAC[NOT_CONS_NIL] THEN + REWRITE_TAC[normalize; LAST] THEN REPEAT GEN_TAC THEN + REPEAT(COND_CASES_TAC THEN ASM_SIMP_TAC[normalize]));; + +let NORMAL_DEGREE = prove + (`!p. ~(LAST p = Cx(&0)) ==> (degree p = LENGTH p - 1)`, + SIMP_TAC[degree; NORMALIZE_EQ] THEN REPEAT STRIP_TAC THEN ARITH_TAC);; + +let LAST_LINEAR_MUL_LEMMA = prove + (`!p a b x. + LAST(a ## p ++ CONS x (b ## p)) = if p = [] then x else b * LAST(p)`, + LIST_INDUCT_TAC THEN + REWRITE_TAC[poly_cmul; poly_add; LAST; HD; TL; NOT_CONS_NIL] THEN + REPEAT GEN_TAC THEN + SUBGOAL_THEN `~(a ## t ++ CONS (b * h) (b ## t) = [])` + ASSUME_TAC THENL + [SPEC_TAC(`t:complex list`,`t:complex list`) THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[poly_cmul; poly_add; NOT_CONS_NIL]; + ALL_TAC] THEN + ASM_REWRITE_TAC[] THEN COND_CASES_TAC THEN ASM_REWRITE_TAC[]);; + +let LAST_LINEAR_MUL = prove + (`!p. ~(p = []) ==> (LAST([a; Cx(&1)] ** p) = LAST p)`, + SIMP_TAC[poly_mul; NOT_CONS_NIL; LAST_LINEAR_MUL_LEMMA; COMPLEX_MUL_LID]);; + +let NORMAL_NORMALIZE = prove + (`!p. ~(normalize p = []) ==> ~(LAST(normalize p) = Cx(&0))`, + LIST_INDUCT_TAC THEN REWRITE_TAC[normalize] THEN + POP_ASSUM MP_TAC THEN ASM_CASES_TAC `normalize t = []` THEN + ASM_REWRITE_TAC[LAST; NOT_CONS_NIL] THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[LAST]);; + +let LINEAR_MUL_DEGREE = prove + (`!p a. ~(poly p = poly []) ==> (degree([a; Cx(&1)] ** p) = degree(p) + 1)`, + REPEAT STRIP_TAC THEN + SUBGOAL_THEN `degree([a; Cx(&1)] ** normalize p) = degree(normalize p) + 1` + MP_TAC THENL + [FIRST_ASSUM(ASSUME_TAC o + GEN_REWRITE_RULE RAND_CONV [POLY_NORMALIZE_ZERO]) THEN + FIRST_ASSUM(MP_TAC o MATCH_MP NORMAL_NORMALIZE) THEN + FIRST_ASSUM(MP_TAC o MATCH_MP LAST_LINEAR_MUL) THEN + SIMP_TAC[NORMAL_DEGREE] THEN REPEAT STRIP_TAC THEN + SUBST1_TAC(SYM(SPEC `a:complex` COMPLEX_NEG_NEG)) THEN + REWRITE_TAC[POLY_LENGTH_MUL] THEN + UNDISCH_TAC `~(normalize p = [])` THEN + SPEC_TAC(`normalize p`,`p:complex list`) THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[NOT_CONS_NIL; LENGTH] THEN ARITH_TAC; + MATCH_MP_TAC EQ_IMP THEN BINOP_TAC THEN + TRY(AP_THM_TAC THEN AP_TERM_TAC) THEN MATCH_MP_TAC DEGREE_WELLDEF THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; POLY_NORMALIZE]]);; + +let LINEAR_POW_MUL_DEGREE = prove + (`!n a p. degree([a; Cx(&1)] exp n ** p) = + if poly p = poly [] then 0 else degree p + n`, + INDUCT_TAC THEN REWRITE_TAC[poly_exp] THENL + [GEN_TAC THEN COND_CASES_TAC THEN ASM_REWRITE_TAC[] THENL + [MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC `degree(p)` THEN CONJ_TAC THENL + [MATCH_MP_TAC DEGREE_WELLDEF THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly; COMPLEX_MUL_RZERO; + COMPLEX_ADD_RID; COMPLEX_MUL_LID]; + MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC `degree []` THEN CONJ_TAC THENL + [MATCH_MP_TAC DEGREE_WELLDEF THEN ASM_REWRITE_TAC[]; + REWRITE_TAC[degree; LENGTH; normalize; ARITH]]]; + REWRITE_TAC[ADD_CLAUSES] THEN MATCH_MP_TAC DEGREE_WELLDEF THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; poly; COMPLEX_MUL_RZERO; + COMPLEX_ADD_RID; COMPLEX_MUL_LID]]; + ALL_TAC] THEN + REPEAT GEN_TAC THEN MATCH_MP_TAC EQ_TRANS THEN + EXISTS_TAC `degree([a; Cx (&1)] exp n ** ([a; Cx (&1)] ** p))` THEN + CONJ_TAC THENL + [MATCH_MP_TAC DEGREE_WELLDEF THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; COMPLEX_MUL_AC]; ALL_TAC] THEN + ASM_REWRITE_TAC[POLY_ENTIRE] THEN + ASM_CASES_TAC `poly p = poly []` THEN ASM_REWRITE_TAC[] THEN + ASM_SIMP_TAC[LINEAR_MUL_DEGREE] THEN + SUBGOAL_THEN `~(poly [a; Cx (&1)] = poly [])` + (fun th -> REWRITE_TAC[th] THEN ARITH_TAC) THEN + REWRITE_TAC[POLY_NORMALIZE_EQ] THEN + REWRITE_TAC[normalize; CX_INJ; REAL_OF_NUM_EQ; ARITH; NOT_CONS_NIL]);; + +(* ------------------------------------------------------------------------- *) +(* Show that the order of a root (or nonroot!) is bounded by degree. *) +(* ------------------------------------------------------------------------- *) + +let ORDER_DEGREE = prove + (`!a p. ~(poly p = poly []) ==> order a p <= degree p`, + REPEAT STRIP_TAC THEN + FIRST_ASSUM(MP_TAC o SPEC `a:complex` o MATCH_MP ORDER_THM) THEN + DISCH_THEN(MP_TAC o REWRITE_RULE[divides] o CONJUNCT1) THEN + DISCH_THEN(X_CHOOSE_THEN `q:complex list` ASSUME_TAC) THEN + FIRST_ASSUM(SUBST1_TAC o MATCH_MP DEGREE_WELLDEF) THEN + ASM_REWRITE_TAC[LINEAR_POW_MUL_DEGREE] THEN + FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [FUN_EQ_THM]) THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[POLY_MUL] THENL + [UNDISCH_TAC `~(poly p = poly [])` THEN + SIMP_TAC[FUN_EQ_THM; POLY_MUL; poly; COMPLEX_MUL_RZERO]; + DISCH_TAC THEN ARITH_TAC]);; + +(* ------------------------------------------------------------------------- *) +(* Tidier versions of finiteness of roots. *) +(* ------------------------------------------------------------------------- *) + +let POLY_ROOTS_FINITE_SET = prove + (`!p. ~(poly p = poly []) ==> FINITE { x | poly p x = Cx(&0)}`, + GEN_TAC THEN REWRITE_TAC[POLY_ROOTS_FINITE] THEN + DISCH_THEN(X_CHOOSE_THEN `N:num` MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `i:num->complex` ASSUME_TAC) THEN + MATCH_MP_TAC FINITE_SUBSET THEN + EXISTS_TAC `{x:complex | ?n:num. n < N /\ (x = i n)}` THEN + CONJ_TAC THENL + [SPEC_TAC(`N:num`,`N:num`) THEN POP_ASSUM_LIST(K ALL_TAC) THEN + INDUCT_TAC THENL + [SUBGOAL_THEN `{x:complex | ?n. n < 0 /\ (x = i n)} = {}` + (fun th -> REWRITE_TAC[th; FINITE_RULES]) THEN + REWRITE_TAC[EXTENSION; IN_ELIM_THM; NOT_IN_EMPTY; LT]; + SUBGOAL_THEN `{x:complex | ?n. n < SUC N /\ (x = i n)} = + (i N) INSERT {x:complex | ?n. n < N /\ (x = i n)}` + SUBST1_TAC THENL + [REWRITE_TAC[EXTENSION; IN_ELIM_THM; IN_INSERT; LT] THEN MESON_TAC[]; + MATCH_MP_TAC(CONJUNCT2 FINITE_RULES) THEN ASM_REWRITE_TAC[]]]; + ASM_REWRITE_TAC[SUBSET; IN_ELIM_THM]]);; + +(* ------------------------------------------------------------------------- *) +(* Conversions to perform operations if coefficients are rational constants. *) +(* ------------------------------------------------------------------------- *) + +let COMPLEX_RAT_MUL_CONV = + GEN_REWRITE_CONV I [GSYM CX_MUL] THENC RAND_CONV REAL_RAT_MUL_CONV;; + +let COMPLEX_RAT_ADD_CONV = + GEN_REWRITE_CONV I [GSYM CX_ADD] THENC RAND_CONV REAL_RAT_ADD_CONV;; + +let COMPLEX_RAT_EQ_CONV = + GEN_REWRITE_CONV I [CX_INJ] THENC REAL_RAT_EQ_CONV;; + +let POLY_CMUL_CONV = + let cmul_conv0 = GEN_REWRITE_CONV I [CONJUNCT1 poly_cmul] + and cmul_conv1 = GEN_REWRITE_CONV I [CONJUNCT2 poly_cmul] in + let rec POLY_CMUL_CONV tm = + (cmul_conv0 ORELSEC + (cmul_conv1 THENC + LAND_CONV COMPLEX_RAT_MUL_CONV THENC + RAND_CONV POLY_CMUL_CONV)) tm in + POLY_CMUL_CONV;; + +let POLY_ADD_CONV = + let add_conv0 = GEN_REWRITE_CONV I (butlast (CONJUNCTS POLY_ADD_CLAUSES)) + and add_conv1 = GEN_REWRITE_CONV I [last (CONJUNCTS POLY_ADD_CLAUSES)] in + let rec POLY_ADD_CONV tm = + (add_conv0 ORELSEC + (add_conv1 THENC + LAND_CONV COMPLEX_RAT_ADD_CONV THENC + RAND_CONV POLY_ADD_CONV)) tm in + POLY_ADD_CONV;; + +let POLY_MUL_CONV = + let mul_conv0 = GEN_REWRITE_CONV I [CONJUNCT1 POLY_MUL_CLAUSES] + and mul_conv1 = GEN_REWRITE_CONV I [CONJUNCT1(CONJUNCT2 POLY_MUL_CLAUSES)] + and mul_conv2 = GEN_REWRITE_CONV I [CONJUNCT2(CONJUNCT2 POLY_MUL_CLAUSES)] in + let rec POLY_MUL_CONV tm = + (mul_conv0 ORELSEC + (mul_conv1 THENC POLY_CMUL_CONV) ORELSEC + (mul_conv2 THENC + LAND_CONV POLY_CMUL_CONV THENC + RAND_CONV(RAND_CONV POLY_MUL_CONV) THENC + POLY_ADD_CONV)) tm in + POLY_MUL_CONV;; + +let POLY_NORMALIZE_CONV = + let pth = prove + (`normalize (CONS h t) = + (\n. if n = [] then if h = Cx(&0) then [] else [h] else CONS h n) + (normalize t)`, + REWRITE_TAC[normalize]) in + let norm_conv0 = GEN_REWRITE_CONV I [CONJUNCT1 normalize] + and norm_conv1 = GEN_REWRITE_CONV I [pth] + and norm_conv2 = GEN_REWRITE_CONV DEPTH_CONV + [COND_CLAUSES; NOT_CONS_NIL; EQT_INTRO(SPEC_ALL EQ_REFL)] in + let rec POLY_NORMALIZE_CONV tm = + (norm_conv0 ORELSEC + (norm_conv1 THENC + RAND_CONV POLY_NORMALIZE_CONV THENC + BETA_CONV THENC + RATOR_CONV(RAND_CONV(RATOR_CONV(LAND_CONV COMPLEX_RAT_EQ_CONV))) THENC + norm_conv2)) tm in + POLY_NORMALIZE_CONV;; + +(* ------------------------------------------------------------------------- *) +(* Now we're finished with polynomials... *) +(* ------------------------------------------------------------------------- *) + +(************** keep this for now + +do_list reduce_interface + ["divides",`poly_divides:complex list->complex list->bool`; + "exp",`poly_exp:complex list -> num -> complex list`; + "diff",`poly_diff:complex list->complex list`];; + +unparse_as_infix "exp";; + + ****************) diff --git a/Complex/fundamental.ml b/Complex/fundamental.ml new file mode 100644 index 0000000..bc72df1 --- /dev/null +++ b/Complex/fundamental.ml @@ -0,0 +1,683 @@ +(* ========================================================================= *) +(* Fundamental theorem of algebra. *) +(* ========================================================================= *) + +needs "Complex/complex_transc.ml";; +needs "Complex/cpoly.ml";; + +prioritize_complex();; + +(* ------------------------------------------------------------------------- *) +(* A cute trick to reduce magnitude of unimodular number. *) +(* ------------------------------------------------------------------------- *) + +let SQRT_SOS_LT_1 = prove + (`!x y. sqrt(x pow 2 + y pow 2) < &1 <=> x pow 2 + y pow 2 < &1`, + REPEAT GEN_TAC THEN + GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM SQRT_1] THEN + REWRITE_TAC[REAL_POW_2] THEN + SIMP_TAC[SQRT_MONO_LT_EQ; REAL_POS; REAL_LE_ADD; REAL_LE_SQUARE]);; + +let SQRT_SOS_EQ_1 = prove + (`!x y. (sqrt(x pow 2 + y pow 2) = &1) <=> (x pow 2 + y pow 2 = &1)`, + REPEAT GEN_TAC THEN + GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM SQRT_1] THEN + REWRITE_TAC[REAL_POW_2] THEN + SIMP_TAC[SQRT_INJ; REAL_POS; REAL_LE_ADD; REAL_LE_SQUARE]);; + +let UNIMODULAR_REDUCE_NORM = prove + (`!z. (norm(z) = &1) + ==> norm(z + Cx(&1)) < &1 \/ + norm(z - Cx(&1)) < &1 \/ + norm(z + ii) < &1 \/ + norm(z - ii) < &1`, + GEN_TAC THEN + REWRITE_TAC[ii; CX_DEF; complex_add; complex_sub; complex_neg; complex_norm; + RE; IM; REAL_ADD_RID; REAL_NEG_0; SQRT_SOS_LT_1; SQRT_SOS_EQ_1] THEN + SIMP_TAC[REAL_POW_2; + REAL_ARITH `a * a + (b + c) * (b + c) = + (a * a + b * b) + (&2 * b * c + c * c)`; + REAL_ARITH `(b + c) * (b + c) + a * a = + (b * b + a * a) + (&2 * b * c + c * c)`] THEN + DISCH_TAC THEN REWRITE_TAC[REAL_ARITH `&1 + x < &1 <=> &0 < --x`] THEN + REWRITE_TAC[REAL_NEG_ADD; REAL_MUL_LNEG; REAL_MUL_RNEG; REAL_NEG_NEG] THEN + REWRITE_TAC[REAL_MUL_RID] THEN + MATCH_MP_TAC(REAL_ARITH + `~(abs(a) <= &1 /\ abs(b) <= &1) + ==> &0 < --a + --(&1) \/ &0 < a + --(&1) \/ + &0 < --b + --(&1) \/ &0 < b + --(&1)`) THEN + STRIP_TAC THEN UNDISCH_TAC `Re z * Re z + Im z * Im z = &1` THEN + REWRITE_TAC[] THEN + MATCH_MP_TAC(REAL_ARITH + `(&2 * r) * (&2 * r) <= &1 /\ (&2 * i) * (&2 * i) <= &1 + ==> ~(r * r + i * i = &1)`) THEN + REWRITE_TAC[GSYM REAL_POW_2] THEN ONCE_REWRITE_TAC[GSYM REAL_POW2_ABS] THEN + ASM_SIMP_TAC[REAL_POW_1_LE; REAL_ABS_POS]);; + +(* ------------------------------------------------------------------------- *) +(* Hence we can always reduce modulus of 1 + b z^n if nonzero *) +(* ------------------------------------------------------------------------- *) + +let REDUCE_POLY_SIMPLE = prove + (`!b n. ~(b = Cx(&0)) /\ ~(n = 0) + ==> ?z. norm(Cx(&1) + b * z pow n) < &1`, + GEN_TAC THEN MATCH_MP_TAC num_WF THEN REPEAT STRIP_TAC THEN + ASM_CASES_TAC `EVEN(n)` THENL + [FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [EVEN_EXISTS]) THEN + DISCH_THEN(X_CHOOSE_THEN `m:num` SUBST_ALL_TAC) THEN + FIRST_X_ASSUM(MP_TAC o SPEC `m:num`) THEN + ASM_SIMP_TAC[ARITH_RULE `~(2 * m = 0) ==> m < 2 * m /\ ~(m = 0)`] THEN + DISCH_THEN(X_CHOOSE_TAC `w:complex`) THEN EXISTS_TAC `csqrt w` THEN + ASM_REWRITE_TAC[GSYM COMPLEX_POW_POW; CSQRT]; ALL_TAC] THEN + MP_TAC(SPEC `Cx(norm b) / b` UNIMODULAR_REDUCE_NORM) THEN ANTS_TAC THENL + [REWRITE_TAC[COMPLEX_NORM_DIV; COMPLEX_NORM_CX] THEN + ASM_SIMP_TAC[COMPLEX_ABS_NORM; REAL_DIV_REFL; COMPLEX_NORM_ZERO]; + ALL_TAC] THEN DISCH_TAC THEN + SUBGOAL_THEN `?v. norm(Cx(norm b) / b + v pow n) < &1` MP_TAC THENL + [SUBGOAL_THEN `(Cx(&1) pow n = Cx(&1)) /\ + (--Cx(&1) pow n = --Cx(&1)) /\ + (((ii pow n = ii) /\ (--ii pow n = --ii)) \/ + ((ii pow n = --ii) /\ (--ii pow n = ii)))` + MP_TAC THENL + [ALL_TAC; + RULE_ASSUM_TAC(REWRITE_RULE[complex_sub]) THEN ASM_MESON_TAC[]] THEN + FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [NOT_EVEN]) THEN + SIMP_TAC[ODD_EXISTS; LEFT_IMP_EXISTS_THM] THEN + X_GEN_TAC `m:num` THEN DISCH_THEN(K ALL_TAC) THEN + REWRITE_TAC[complex_pow; COMPLEX_POW_NEG; EVEN; EVEN_MULT; ARITH_EVEN] THEN + REWRITE_TAC[GSYM COMPLEX_POW_POW] THEN + REWRITE_TAC[COMPLEX_POW_ONE; COMPLEX_POW_II_2; COMPLEX_MUL_LID; + COMPLEX_POW_NEG] THEN + COND_CASES_TAC THEN + REWRITE_TAC[COMPLEX_MUL_RID; COMPLEX_MUL_RNEG; COMPLEX_NEG_NEG]; + ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `v:complex` ASSUME_TAC) THEN + EXISTS_TAC `v / Cx(root(n) (norm b))` THEN + REWRITE_TAC[COMPLEX_POW_DIV; GSYM CX_POW] THEN + SUBGOAL_THEN `root n (norm b) pow n = norm b` SUBST1_TAC THENL + [UNDISCH_TAC `~(EVEN n)` THEN SPEC_TAC(`n:num`,`n:num`) THEN + INDUCT_TAC THEN SIMP_TAC[EVEN; ROOT_POW_POS; COMPLEX_NORM_POS]; + ALL_TAC] THEN + MATCH_MP_TAC REAL_LT_LCANCEL_IMP THEN EXISTS_TAC `norm(Cx(norm b) / b)` THEN + REWRITE_TAC[GSYM COMPLEX_NORM_MUL; COMPLEX_ADD_LDISTRIB] THEN + REWRITE_TAC[COMPLEX_MUL_RID; REAL_MUL_RID] THEN + SUBGOAL_THEN `norm(Cx(norm b) / b) = &1` SUBST1_TAC THENL + [REWRITE_TAC[COMPLEX_NORM_DIV; COMPLEX_NORM_CX; COMPLEX_ABS_NORM] THEN + ASM_SIMP_TAC[REAL_DIV_REFL; COMPLEX_NORM_ZERO]; ALL_TAC] THEN + REWRITE_TAC[REAL_LT_01; complex_div] THEN + ONCE_REWRITE_TAC[AC COMPLEX_MUL_AC + `(m * b') * b * p * m' = (m * m') * (b * b') * p`] THEN + ASM_SIMP_TAC[COMPLEX_MUL_RINV; COMPLEX_MUL_LID; + CX_INJ; COMPLEX_NORM_ZERO] THEN + ASM_REWRITE_TAC[GSYM complex_div]);; + +(* ------------------------------------------------------------------------- *) +(* Basic lemmas about polynomials. *) +(* ------------------------------------------------------------------------- *) + +let POLY_CMUL_MAP = prove + (`!p c x. poly (MAP (( * ) c) p) x = c * poly p x`, + LIST_INDUCT_TAC THEN REWRITE_TAC[MAP; poly; COMPLEX_MUL_RZERO] THEN + ASM_REWRITE_TAC[COMPLEX_ADD_LDISTRIB] THEN REWRITE_TAC[COMPLEX_MUL_AC]);; + +let POLY_0 = prove + (`!p x. ALL (\b. b = Cx(&0)) p ==> (poly p x = Cx(&0))`, + LIST_INDUCT_TAC THEN + ASM_SIMP_TAC[ALL; poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID]);; + +let POLY_BOUND_EXISTS = prove + (`!p r. ?m. &0 < m /\ !z. norm(z) <= r ==> norm(poly p z) <= m`, + ONCE_REWRITE_TAC[SWAP_FORALL_THM] THEN GEN_TAC THEN + LIST_INDUCT_TAC THENL + [EXISTS_TAC `&1` THEN REWRITE_TAC[poly; COMPLEX_NORM_CX] THEN + REWRITE_TAC[REAL_ABS_NUM; REAL_LT_01; REAL_POS]; ALL_TAC] THEN + POP_ASSUM(X_CHOOSE_THEN `m:real` STRIP_ASSUME_TAC) THEN + EXISTS_TAC `&1 + norm(h) + abs(r * m)` THEN + ASM_SIMP_TAC[REAL_ARITH `&0 <= x /\ &0 <= y ==> &0 < &1 + x + y`; + REAL_ABS_POS; COMPLEX_NORM_POS] THEN + X_GEN_TAC `z:complex` THEN DISCH_TAC THEN + REWRITE_TAC[poly] THEN MATCH_MP_TAC REAL_LE_TRANS THEN + EXISTS_TAC `norm(h) + norm(z * poly t z)` THEN + REWRITE_TAC[COMPLEX_NORM_TRIANGLE] THEN + MATCH_MP_TAC(REAL_ARITH `y <= z ==> x + y <= &1 + x + abs(z)`) THEN + REWRITE_TAC[COMPLEX_NORM_MUL] THEN MATCH_MP_TAC REAL_LE_MUL2 THEN + ASM_SIMP_TAC[COMPLEX_NORM_POS]);; + +(* ------------------------------------------------------------------------- *) +(* Offsetting the variable in a polynomial gives another of same degree. *) +(* ------------------------------------------------------------------------- *) + +let POLY_OFFSET_LEMMA = prove + (`!a p. ?b q. (LENGTH q = LENGTH p) /\ + !x. poly (CONS b q) x = (a + x) * poly p x`, + GEN_TAC THEN LIST_INDUCT_TAC THENL + [EXISTS_TAC `Cx(&0)` THEN EXISTS_TAC `[]:complex list` THEN + REWRITE_TAC[poly; LENGTH; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID]; + ALL_TAC] THEN + POP_ASSUM STRIP_ASSUME_TAC THEN + EXISTS_TAC `a * h` THEN EXISTS_TAC `CONS (b + h) q` THEN + ASM_REWRITE_TAC[LENGTH; poly] THEN X_GEN_TAC `x:complex ` THEN + FIRST_ASSUM(MP_TAC o SPEC `x:complex`) THEN + REWRITE_TAC[poly] THEN DISCH_THEN(MP_TAC o AP_TERM `( * ) x`) THEN + SIMPLE_COMPLEX_ARITH_TAC);; + +let POLY_OFFSET = prove + (`!a p. ?q. (LENGTH q = LENGTH p) /\ !x. poly q x = poly p (a + x)`, + GEN_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[LENGTH; poly] THENL + [EXISTS_TAC `[]:complex list` THEN REWRITE_TAC[poly; LENGTH]; ALL_TAC] THEN + POP_ASSUM(X_CHOOSE_THEN `p:complex list` (STRIP_ASSUME_TAC o GSYM)) THEN + ASM_REWRITE_TAC[] THEN + MP_TAC(SPECL [`a:complex`; `p:complex list`] POLY_OFFSET_LEMMA) THEN + DISCH_THEN(X_CHOOSE_THEN `b:complex` (X_CHOOSE_THEN `r: complex list` + (STRIP_ASSUME_TAC o GSYM))) THEN + EXISTS_TAC `CONS (h + b) r` THEN ASM_REWRITE_TAC[LENGTH] THEN + REWRITE_TAC[poly] THEN SIMPLE_COMPLEX_ARITH_TAC);; + +(* ------------------------------------------------------------------------- *) +(* Bolzano-Weierstrass type property for closed disc in complex plane. *) +(* ------------------------------------------------------------------------- *) + +let METRIC_BOUND_LEMMA = prove + (`!x y. norm(x - y) <= abs(Re(x) - Re(y)) + abs(Im(x) - Im(y))`, + REPEAT GEN_TAC THEN REWRITE_TAC[complex_norm] THEN + MATCH_MP_TAC(REAL_ARITH + `a <= abs(abs x + abs y) ==> a <= abs x + abs y`) THEN + GEN_REWRITE_TAC RAND_CONV [GSYM POW_2_SQRT_ABS] THEN + MATCH_MP_TAC SQRT_MONO_LE THEN + SIMP_TAC[REAL_POW_2; REAL_LE_ADD; REAL_LE_SQUARE] THEN + REWRITE_TAC[complex_add; complex_sub; complex_neg; RE; IM] THEN + REWRITE_TAC[GSYM real_sub] THEN + REWRITE_TAC[REAL_ARITH `(a + b) * (a + b) = a * a + b * b + &2 * a * b`] THEN + REWRITE_TAC[GSYM REAL_ABS_MUL] THEN + REWRITE_TAC[REAL_ARITH `a + b <= abs a + abs b + &2 * abs c`]);; + +let BOLZANO_WEIERSTRASS_COMPLEX_DISC = prove + (`!s r. (!n. norm(s n) <= r) + ==> ?f z. subseq f /\ + !e. &0 < e ==> ?N. !n. n >= N ==> norm(s(f n) - z) < e`, + REPEAT STRIP_TAC THEN + MP_TAC(SPEC `Re o (s:num->complex)` SEQ_MONOSUB) THEN + DISCH_THEN(X_CHOOSE_THEN `f:num->num` MP_TAC) THEN + REWRITE_TAC[o_THM] THEN STRIP_TAC THEN + MP_TAC(SPEC `Im o (s:num->complex) o (f:num->num)` SEQ_MONOSUB) THEN + DISCH_THEN(X_CHOOSE_THEN `g:num->num` MP_TAC) THEN + REWRITE_TAC[o_THM] THEN STRIP_TAC THEN + EXISTS_TAC `(f:num->num) o (g:num->num)` THEN + SUBGOAL_THEN `convergent (\n. Re(s(f n :num))) /\ + convergent (\n. Im(s((f:num->num)(g n))))` + MP_TAC THENL + [CONJ_TAC THEN MATCH_MP_TAC SEQ_BCONV THEN ASM_REWRITE_TAC[bounded] THEN + MAP_EVERY EXISTS_TAC [`r + &1`; `&0`; `0`] THEN + REWRITE_TAC[GE; LE_0; MR1_DEF; REAL_SUB_LZERO; REAL_ABS_NEG] THEN + X_GEN_TAC `n:num` THEN + W(fun (_,w) -> FIRST_ASSUM(MP_TAC o SPEC (funpow 3 rand (lhand w)))) THEN + REWRITE_TAC[complex_norm] THEN + MATCH_MP_TAC(REAL_ARITH `a <= b ==> b <= r ==> a < r + &1`) THEN + GEN_REWRITE_TAC LAND_CONV [GSYM POW_2_SQRT_ABS] THEN + MATCH_MP_TAC SQRT_MONO_LE THEN + REWRITE_TAC[REAL_POW_2; REAL_LE_SQUARE; REAL_LE_ADDR; REAL_LE_ADDL]; + ALL_TAC] THEN + REWRITE_TAC[convergent; SEQ; GE] THEN + DISCH_THEN(CONJUNCTS_THEN2 + (X_CHOOSE_TAC `x:real`) (X_CHOOSE_TAC `y:real`)) THEN + EXISTS_TAC `complex(x,y)` THEN CONJ_TAC THENL + [MAP_EVERY UNDISCH_TAC [`subseq f`; `subseq g`] THEN + REWRITE_TAC[subseq; o_THM] THEN MESON_TAC[]; ALL_TAC] THEN + X_GEN_TAC `e:real` THEN DISCH_TAC THEN + UNDISCH_TAC + `!e. &0 < e + ==> (?N. !n. N <= n ==> abs(Re(s ((f:num->num) n)) - x) < e)` THEN + DISCH_THEN(MP_TAC o SPEC `e / &2`) THEN + FIRST_X_ASSUM(MP_TAC o SPEC `e / &2`) THEN + ASM_SIMP_TAC[REAL_LT_DIV; REAL_OF_NUM_LT; ARITH] THEN + DISCH_THEN(X_CHOOSE_TAC `N2:num`) THEN DISCH_THEN(X_CHOOSE_TAC `N1:num`) THEN + EXISTS_TAC `N1 + N2:num` THEN X_GEN_TAC `n:num` THEN DISCH_TAC THEN + MATCH_MP_TAC REAL_LTE_TRANS THEN EXISTS_TAC `&2 * e / &2` THEN + CONJ_TAC THENL + [ALL_TAC; + SIMP_TAC[REAL_DIV_LMUL; REAL_OF_NUM_EQ; ARITH; REAL_LE_REFL]] THEN + W(MP_TAC o PART_MATCH lhand METRIC_BOUND_LEMMA o lhand o snd) THEN + MATCH_MP_TAC(REAL_ARITH + `a < c /\ b < c ==> x <= a + b ==> x < &2 * c`) THEN + REWRITE_TAC[o_THM; RE; IM] THEN CONJ_TAC THEN FIRST_ASSUM MATCH_MP_TAC THEN + ASM_MESON_TAC[LE_ADD; SEQ_SUBLE; LE_TRANS; ADD_SYM]);; + +(* ------------------------------------------------------------------------- *) +(* Polynomial is continuous. *) +(* ------------------------------------------------------------------------- *) + +let POLY_CONT = prove + (`!p z e. &0 < e + ==> ?d. &0 < d /\ !w. &0 < norm(w - z) /\ norm(w - z) < d + ==> norm(poly p w - poly p z) < e`, + REPEAT STRIP_TAC THEN + MP_TAC(SPECL [`z:complex`; `p:complex list`] POLY_OFFSET) THEN + DISCH_THEN(X_CHOOSE_THEN `q:complex list` (MP_TAC o CONJUNCT2)) THEN + DISCH_THEN(MP_TAC o GEN `w:complex` o SYM o SPEC `w - z`) THEN + REWRITE_TAC[COMPLEX_SUB_ADD2] THEN + DISCH_THEN(fun th -> ONCE_REWRITE_TAC[th]) THEN + REWRITE_TAC[COMPLEX_SUB_REFL] THEN + SPEC_TAC(`q:complex list`,`p:complex list`) THEN + LIST_INDUCT_TAC THEN ASM_REWRITE_TAC[poly] THENL + [EXISTS_TAC `e:real` THEN + ASM_REWRITE_TAC[COMPLEX_SUB_REFL; COMPLEX_NORM_CX; REAL_ABS_NUM]; + ALL_TAC] THEN + REWRITE_TAC[COMPLEX_MUL_LZERO; COMPLEX_ADD_RID; COMPLEX_ADD_SUB] THEN + MP_TAC(SPECL [`t:complex list`; `&1`] POLY_BOUND_EXISTS) THEN + DISCH_THEN(X_CHOOSE_THEN `m:real` STRIP_ASSUME_TAC) THEN + MP_TAC(SPECL [`&1`; `e / m:real`] REAL_DOWN2) THEN + ASM_SIMP_TAC[REAL_LT_DIV; REAL_LT_01] THEN + MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `d:real` THEN + STRIP_TAC THEN ASM_REWRITE_TAC[] THEN X_GEN_TAC `w:complex` THEN + STRIP_TAC THEN REWRITE_TAC[COMPLEX_NORM_MUL] THEN + MATCH_MP_TAC REAL_LET_TRANS THEN EXISTS_TAC `d * m:real` THEN + ASM_SIMP_TAC[GSYM REAL_LT_RDIV_EQ] THEN MATCH_MP_TAC REAL_LE_MUL2 THEN + ASM_MESON_TAC[REAL_LT_TRANS; REAL_LT_IMP_LE; COMPLEX_NORM_POS]);; + +(* ------------------------------------------------------------------------- *) +(* Hence a polynomial attains minimum on a closed disc in the complex plane. *) +(* ------------------------------------------------------------------------- *) + +let POLY_MINIMUM_MODULUS_DISC = prove + (`!p r. ?z. !w. norm(w) <= r ==> norm(poly p z) <= norm(poly p w)`, + let lemma = prove + (`P /\ (m = --x) /\ --x < y <=> (--x = m) /\ P /\ m < y`, + MESON_TAC[]) in + REPEAT GEN_TAC THEN + ASM_CASES_TAC `&0 <= r` THENL + [ALL_TAC; ASM_MESON_TAC[COMPLEX_NORM_POS; REAL_LE_TRANS]] THEN + MP_TAC(SPEC `\x. ?z. norm(z) <= r /\ (norm(poly p z) = --x)` + REAL_SUP_EXISTS) THEN + REWRITE_TAC[] THEN ANTS_TAC THENL + [CONJ_TAC THENL + [MAP_EVERY EXISTS_TAC [`--norm(poly p (Cx(&0)))`; `Cx(&0)`] THEN + ASM_REWRITE_TAC[REAL_NEG_NEG; COMPLEX_NORM_CX; REAL_ABS_NUM]; + EXISTS_TAC `&1` THEN + REWRITE_TAC[REAL_ARITH `(a = --b) <=> (--b = a:real)`] THEN + REWRITE_TAC[REAL_ARITH `x < &1 <=> --(&1) < --x`] THEN + SIMP_TAC[LEFT_IMP_EXISTS_THM] THEN + SIMP_TAC[REAL_ARITH `&0 <= x ==> --(&1) < x`; COMPLEX_NORM_POS]]; + ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `s:real` MP_TAC) THEN + ONCE_REWRITE_TAC[REAL_ARITH `a < b <=> --b < --a:real`] THEN + ABBREV_TAC `m = --s:real` THEN + DISCH_THEN(MP_TAC o GEN `y:real` o SPEC `--y:real`) THEN + REWRITE_TAC[REAL_NEG_NEG] THEN + REWRITE_TAC[LEFT_AND_EXISTS_THM; GSYM CONJ_ASSOC; lemma] THEN + ONCE_REWRITE_TAC[SWAP_EXISTS_THM] THEN + ONCE_REWRITE_TAC[REAL_ARITH `(--a = b) <=> (a = --b:real)`] THEN + REWRITE_TAC[LEFT_EXISTS_AND_THM; EXISTS_REFL] THEN + DISCH_THEN(fun th -> MP_TAC th THEN MP_TAC(SPEC `m:real` th)) THEN + REWRITE_TAC[REAL_LT_REFL; NOT_EXISTS_THM] THEN + REWRITE_TAC[TAUT `~(a /\ b) <=> a ==> ~b`] THEN + REWRITE_TAC[REAL_NOT_LT] THEN DISCH_TAC THEN + DISCH_THEN(MP_TAC o GEN `n:num` o SPEC `m + inv(&(SUC n))`) THEN + REWRITE_TAC[REAL_LT_ADDR; REAL_LT_INV_EQ; REAL_OF_NUM_LT; LT_0] THEN + REWRITE_TAC[SKOLEM_THM; FORALL_AND_THM] THEN + DISCH_THEN(X_CHOOSE_THEN `s:num->complex` STRIP_ASSUME_TAC) THEN + MP_TAC(SPECL [`s:num->complex`; `r:real`] + BOLZANO_WEIERSTRASS_COMPLEX_DISC) THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_THEN `f:num->num` (X_CHOOSE_THEN `z:complex` + STRIP_ASSUME_TAC)) THEN + EXISTS_TAC `z:complex` THEN X_GEN_TAC `w:complex` THEN DISCH_TAC THEN + SUBGOAL_THEN `norm(poly p z) = m` (fun th -> ASM_SIMP_TAC[th]) THEN + MATCH_MP_TAC(REAL_ARITH `~(&0 < abs(a - b)) ==> (a = b)`) THEN DISCH_TAC THEN + ABBREV_TAC `e = abs(norm(poly p z) - m)` THEN + MP_TAC(SPECL [`p:complex list`; `z:complex`; `e / &2`] POLY_CONT) THEN + ASM_SIMP_TAC[REAL_LT_DIV; REAL_OF_NUM_LT; ARITH] THEN + DISCH_THEN(X_CHOOSE_THEN `d:real` STRIP_ASSUME_TAC) THEN + SUBGOAL_THEN `!w. norm(w - z) < d ==> norm(poly p w - poly p z) < e / &2` + MP_TAC THENL + [X_GEN_TAC `u:complex` THEN + ASM_CASES_TAC `u = z:complex` THEN + ASM_SIMP_TAC[COMPLEX_SUB_REFL; REAL_LT_DIV; REAL_OF_NUM_LT; ARITH; + COMPLEX_NORM_CX; REAL_ABS_NUM] THEN + DISCH_TAC THEN FIRST_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[] THEN + ASM_REWRITE_TAC[COMPLEX_NORM_NZ; COMPLEX_SUB_0]; ALL_TAC] THEN + FIRST_ASSUM(K ALL_TAC o check (is_conj o lhand o + snd o dest_forall o concl)) THEN + DISCH_TAC THEN + FIRST_X_ASSUM(MP_TAC o SPEC `d:real`) THEN ASM_REWRITE_TAC[GE] THEN + DISCH_THEN(X_CHOOSE_THEN `N1:num` ASSUME_TAC) THEN + MP_TAC(SPEC `&2 / e` REAL_ARCH_SIMPLE) THEN + DISCH_THEN(X_CHOOSE_THEN `N2:num` ASSUME_TAC) THEN + SUBGOAL_THEN `norm(poly p (s((f:num->num) (N1 + N2))) - poly p z) < e / &2` + MP_TAC THENL + [FIRST_ASSUM MATCH_MP_TAC THEN ASM_SIMP_TAC[LE_ADD]; ALL_TAC] THEN + MATCH_MP_TAC(REAL_ARITH + `!m. abs(norm(psfn) - m) < e2 /\ + &2 * e2 <= abs(norm(psfn) - m) + norm(psfn - pz) + ==> norm(psfn - pz) < e2 ==> F`) THEN + EXISTS_TAC `m:real` THEN CONJ_TAC THENL + [MATCH_MP_TAC REAL_LTE_TRANS THEN EXISTS_TAC `inv(&(SUC(N1 + N2)))` THEN + CONJ_TAC THENL + [MATCH_MP_TAC(REAL_ARITH + `m <= x /\ x < m + e ==> abs(x - m:real) < e`) THEN + ASM_SIMP_TAC[] THEN + MATCH_MP_TAC REAL_LTE_TRANS THEN + EXISTS_TAC `m + inv(&(SUC(f(N1 + N2:num))))` THEN + ASM_REWRITE_TAC[REAL_LE_LADD] THEN MATCH_MP_TAC REAL_LE_INV2 THEN + ASM_SIMP_TAC[REAL_OF_NUM_LT; REAL_OF_NUM_LE; LT_0; LE_SUC; SEQ_SUBLE]; + GEN_REWRITE_TAC RAND_CONV [GSYM REAL_INV_DIV] THEN + MATCH_MP_TAC REAL_LE_INV2 THEN + ASM_SIMP_TAC[REAL_LT_DIV; REAL_OF_NUM_LT; ARITH] THEN + MATCH_MP_TAC REAL_LE_TRANS THEN EXISTS_TAC `&N2` THEN + ASM_REWRITE_TAC[REAL_OF_NUM_LE] THEN ARITH_TAC]; ALL_TAC] THEN + SIMP_TAC[REAL_DIV_LMUL; REAL_OF_NUM_EQ; ARITH] THEN + EXPAND_TAC "e" THEN + MATCH_MP_TAC(REAL_ARITH + `abs(norm(psfn) - norm(pz)) <= norm(psfn - pz) + ==> abs(norm(pz) - m) <= abs(norm(psfn) - m) + norm(psfn - pz)`) THEN + REWRITE_TAC[COMPLEX_NORM_ABS_NORM]);; + +(* ------------------------------------------------------------------------- *) +(* Nonzero polynomial in z goes to infinity as z does. *) +(* ------------------------------------------------------------------------- *) + +let POLY_INFINITY = prove + (`!p a. EX (\b. ~(b = Cx(&0))) p + ==> !d. ?r. !z. r <= norm(z) ==> d <= norm(poly (CONS a p) z)`, + LIST_INDUCT_TAC THEN REWRITE_TAC[EX] THEN X_GEN_TAC `a:complex` THEN + ASM_CASES_TAC `EX (\b. ~(b = Cx(&0))) t` THEN ASM_REWRITE_TAC[] THENL + [X_GEN_TAC `d:real` THEN + FIRST_X_ASSUM(MP_TAC o SPEC `h:complex`) THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_TAC `r:real` o SPEC `d + norm(a)`) THEN + EXISTS_TAC `&1 + abs(r)` THEN X_GEN_TAC `z:complex` THEN DISCH_TAC THEN + ONCE_REWRITE_TAC[poly] THEN + MATCH_MP_TAC REAL_LE_TRANS THEN + EXISTS_TAC `norm(z * poly (CONS h t) z) - norm(a)` THEN CONJ_TAC THENL + [ALL_TAC; + ONCE_REWRITE_TAC[COMPLEX_ADD_SYM] THEN + REWRITE_TAC[REAL_LE_SUB_RADD; COMPLEX_NORM_TRIANGLE_SUB]] THEN + REWRITE_TAC[REAL_LE_SUB_LADD] THEN + MATCH_MP_TAC REAL_LE_TRANS THEN + EXISTS_TAC `&1 * norm(poly (CONS h t) z)` THEN CONJ_TAC THENL + [REWRITE_TAC[REAL_MUL_LID] THEN FIRST_ASSUM MATCH_MP_TAC THEN + ASM_SIMP_TAC[REAL_ARITH `&1 + abs(r) <= x ==> r <= x`]; + REWRITE_TAC[COMPLEX_NORM_MUL] THEN MATCH_MP_TAC REAL_LE_RMUL THEN + REWRITE_TAC[COMPLEX_NORM_POS] THEN + ASM_MESON_TAC[REAL_ARITH `&1 + abs(r) <= x ==> &1 <= x`]]; + RULE_ASSUM_TAC(REWRITE_RULE[NOT_EX]) THEN + ASM_SIMP_TAC[poly; POLY_0; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + DISCH_TAC THEN X_GEN_TAC `d:real` THEN + EXISTS_TAC `(abs(d) + norm(a)) / norm(h)` THEN X_GEN_TAC `z:complex` THEN + ASM_SIMP_TAC[REAL_LE_LDIV_EQ; COMPLEX_NORM_NZ; GSYM COMPLEX_NORM_MUL] THEN + MATCH_MP_TAC(REAL_ARITH + `mzh <= mazh + ma ==> abs(d) + ma <= mzh ==> d <= mazh`) THEN + ONCE_REWRITE_TAC[COMPLEX_ADD_SYM] THEN + REWRITE_TAC[COMPLEX_NORM_TRIANGLE_SUB]]);; + +(* ------------------------------------------------------------------------- *) +(* Hence polynomial's modulus attains its minimum somewhere. *) +(* ------------------------------------------------------------------------- *) + +let POLY_MINIMUM_MODULUS = prove + (`!p. ?z. !w. norm(poly p z) <= norm(poly p w)`, + LIST_INDUCT_TAC THEN REWRITE_TAC[poly; REAL_LE_REFL] THEN + ASM_CASES_TAC `EX (\b. ~(b = Cx(&0))) t` THENL + [FIRST_ASSUM(MP_TAC o SPEC `h:complex` o MATCH_MP POLY_INFINITY) THEN + DISCH_THEN(MP_TAC o SPEC `norm(poly (CONS h t) (Cx(&0)))`) THEN + DISCH_THEN(X_CHOOSE_THEN `r:real` ASSUME_TAC) THEN + MP_TAC(SPECL [`CONS (h:complex) t`; `abs(r)`] + POLY_MINIMUM_MODULUS_DISC) THEN + REWRITE_TAC[GSYM(CONJUNCT2 poly)] THEN + ASM_MESON_TAC[REAL_ARITH `r <= z \/ z <= abs(r)`; REAL_LE_TRANS; + COMPLEX_NORM_CX; REAL_ABS_NUM; REAL_ABS_POS]; + FIRST_ASSUM(MP_TAC o GEN_REWRITE_RULE I [NOT_EX]) THEN + REWRITE_TAC[] THEN DISCH_THEN(ASSUME_TAC o MATCH_MP POLY_0) THEN + ASM_REWRITE_TAC[COMPLEX_MUL_RZERO; REAL_LE_REFL]]);; + +(* ------------------------------------------------------------------------- *) +(* Constant function (non-syntactic characterization). *) +(* ------------------------------------------------------------------------- *) + +let constant = new_definition + `constant f = !w z. f(w) = f(z)`;; + +let NONCONSTANT_LENGTH = prove + (`!p. ~constant(poly p) ==> 2 <= LENGTH p`, + REWRITE_TAC[constant] THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[poly] THEN + REWRITE_TAC[LENGTH; ARITH_RULE `2 <= SUC n <=> ~(n = 0)`] THEN + SIMP_TAC[TAUT `~a ==> ~b <=> b ==> a`; LENGTH_EQ_NIL; poly] THEN + REWRITE_TAC[COMPLEX_MUL_RZERO]);; + +(* ------------------------------------------------------------------------- *) +(* Decomposition of polynomial, skipping zero coefficients after the first. *) +(* ------------------------------------------------------------------------- *) + +let POLY_DECOMPOSE_LEMMA = prove + (`!p. ~(!z. ~(z = Cx(&0)) ==> (poly p z = Cx(&0))) + ==> ?k a q. ~(a = Cx(&0)) /\ + (SUC(LENGTH q + k) = LENGTH p) /\ + !z. poly p z = z pow k * poly (CONS a q) z`, + LIST_INDUCT_TAC THENL [REWRITE_TAC[poly]; ALL_TAC] THEN + ASM_CASES_TAC `h = Cx(&0)` THENL + [GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [poly] THEN + ASM_SIMP_TAC[COMPLEX_ADD_LID; COMPLEX_ENTIRE] THEN + DISCH_THEN(ANTE_RES_THEN MP_TAC) THEN + DISCH_THEN(X_CHOOSE_THEN `k:num` (X_CHOOSE_THEN `a:complex` + (X_CHOOSE_THEN `q:complex list` STRIP_ASSUME_TAC))) THEN + MAP_EVERY EXISTS_TAC [`k + 1`; `a:complex`; `q:complex list`] THEN + ASM_REWRITE_TAC[poly; LENGTH; GSYM ADD1; ADD_CLAUSES] THEN + REWRITE_TAC[COMPLEX_ADD_LID; complex_pow; GSYM COMPLEX_MUL_ASSOC]; + DISCH_THEN(K ALL_TAC) THEN + MAP_EVERY EXISTS_TAC [`0`; `h:complex`; `t:complex list`] THEN + ASM_REWRITE_TAC[complex_pow; COMPLEX_MUL_LID; ADD_CLAUSES; LENGTH]]);; + +let POLY_DECOMPOSE = prove + (`!p. ~constant(poly p) + ==> ?k a q. ~(a = Cx(&0)) /\ ~(k = 0) /\ + (LENGTH q + k + 1 = LENGTH p) /\ + !z. poly p z = poly p (Cx(&0)) + + z pow k * poly (CONS a q) z`, + LIST_INDUCT_TAC THENL [REWRITE_TAC[constant; poly]; ALL_TAC] THEN + POP_ASSUM(K ALL_TAC) THEN DISCH_TAC THEN + MP_TAC(SPEC `t:complex list` POLY_DECOMPOSE_LEMMA) THEN ANTS_TAC THENL + [POP_ASSUM MP_TAC THEN REWRITE_TAC[constant; poly] THEN + REWRITE_TAC[TAUT `~b ==> ~a <=> a ==> b`; COMPLEX_EQ_ADD_LCANCEL] THEN + SIMP_TAC[TAUT `~a ==> b <=> a \/ b`; GSYM COMPLEX_ENTIRE]; ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `k:num` (X_CHOOSE_THEN `a:complex` + (X_CHOOSE_THEN `q:complex list` STRIP_ASSUME_TAC))) THEN + MAP_EVERY EXISTS_TAC [`SUC k`; `a:complex`; `q:complex list`] THEN + ASM_REWRITE_TAC[ADD_CLAUSES; GSYM ADD1; LENGTH; NOT_SUC] THEN + ASM_REWRITE_TAC[poly; COMPLEX_MUL_LZERO; COMPLEX_ADD_RID; complex_pow] THEN + REWRITE_TAC[GSYM COMPLEX_MUL_ASSOC]);; + +let POLY_REPLICATE_APPEND = prove + (`!n p x. poly (APPEND (REPLICATE n (Cx(&0))) p) x = x pow n * poly p x`, + INDUCT_TAC THEN + REWRITE_TAC[REPLICATE; APPEND; complex_pow; COMPLEX_MUL_LID] THEN + ASM_REWRITE_TAC[poly; COMPLEX_ADD_LID] THEN REWRITE_TAC[COMPLEX_MUL_ASSOC]);; + +(* ------------------------------------------------------------------------- *) +(* Fundamental theorem. *) +(* ------------------------------------------------------------------------- *) + +let FUNDAMENTAL_THEOREM_OF_ALGEBRA = prove + (`!p. ~constant(poly p) ==> ?z. poly p z = Cx(&0)`, + SUBGOAL_THEN + `!n p. (LENGTH p = n) /\ ~constant(poly p) ==> ?z. poly p z = Cx(&0)` + (fun th -> MESON_TAC[th]) THEN + MATCH_MP_TAC num_WF THEN + X_GEN_TAC `n:num` THEN STRIP_TAC THEN + X_GEN_TAC `p:complex list` THEN STRIP_TAC THEN + FIRST_ASSUM(MP_TAC o MATCH_MP NONCONSTANT_LENGTH) THEN + ASM_REWRITE_TAC[] THEN DISCH_TAC THEN + X_CHOOSE_TAC `c:complex` (SPEC `p:complex list` POLY_MINIMUM_MODULUS) THEN + ASM_CASES_TAC `poly p c = Cx(&0)` THENL [ASM_MESON_TAC[]; ALL_TAC] THEN + MP_TAC(SPECL [`c:complex`; `p:complex list`] POLY_OFFSET) THEN + DISCH_THEN(X_CHOOSE_THEN `q:complex list` MP_TAC) THEN + DISCH_THEN(CONJUNCTS_THEN2 (SUBST_ALL_TAC o SYM) ASSUME_TAC) THEN + SUBGOAL_THEN `~constant(poly q)` ASSUME_TAC THENL + [UNDISCH_TAC `~(constant(poly p))` THEN + SUBGOAL_THEN `!z. poly q (z - c) = poly p z` + (fun th -> MESON_TAC[th; constant]) THEN + ASM_MESON_TAC[SIMPLE_COMPLEX_ARITH `a + (x - a) = x`]; ALL_TAC] THEN + SUBGOAL_THEN `poly p c = poly q (Cx(&0))` SUBST_ALL_TAC THENL + [ASM_MESON_TAC[COMPLEX_ADD_RID]; ALL_TAC] THEN + SUBGOAL_THEN `!w. norm(poly q (Cx(&0))) <= norm(poly q w)` MP_TAC THENL + [ASM_MESON_TAC[]; ALL_TAC] THEN + POP_ASSUM_LIST(MAP_EVERY (fun th -> + if free_in `p:complex list` (concl th) + then ALL_TAC else ASSUME_TAC th)) THEN + MATCH_MP_TAC(TAUT `~a ==> a ==> b`) THEN + REWRITE_TAC[NOT_FORALL_THM; REAL_NOT_LE] THEN + ABBREV_TAC `a0 = poly q (Cx(&0))` THEN + SUBGOAL_THEN + `!z. poly q z = poly (MAP (( * ) (inv(a0))) q) z * a0` + ASSUME_TAC THENL + [REWRITE_TAC[POLY_CMUL_MAP] THEN + ONCE_REWRITE_TAC[AC COMPLEX_MUL_AC `(a * b) * c = b * c * a`] THEN + ASM_SIMP_TAC[COMPLEX_MUL_RINV; COMPLEX_MUL_RID]; + ALL_TAC] THEN + ABBREV_TAC `r = MAP (( * ) (inv(a0))) q` THEN + SUBGOAL_THEN `LENGTH(q:complex list) = LENGTH(r:complex list)` + SUBST_ALL_TAC THENL + [EXPAND_TAC "r" THEN REWRITE_TAC[LENGTH_MAP]; ALL_TAC] THEN + SUBGOAL_THEN `~(constant(poly r))` ASSUME_TAC THENL + [UNDISCH_TAC `~constant(poly q)` THEN + ASM_REWRITE_TAC[constant; COMPLEX_EQ_MUL_RCANCEL] THEN MESON_TAC[]; + ALL_TAC] THEN + SUBGOAL_THEN `poly r (Cx(&0)) = Cx(&1)` ASSUME_TAC THENL + [FIRST_X_ASSUM(MP_TAC o SPEC `Cx(&0)`) THEN + ASM_REWRITE_TAC[] THEN + GEN_REWRITE_TAC (LAND_CONV o LAND_CONV) [GSYM COMPLEX_MUL_LID] THEN + ASM_SIMP_TAC[COMPLEX_EQ_MUL_RCANCEL]; ALL_TAC] THEN + ASM_REWRITE_TAC[] THEN + POP_ASSUM_LIST(MAP_EVERY (fun th -> + if free_in `q:complex list` (concl th) + then ALL_TAC else ASSUME_TAC th)) THEN + ASM_SIMP_TAC[GSYM REAL_LT_RDIV_EQ; COMPLEX_NORM_NZ; COMPLEX_NORM_MUL; + REAL_DIV_REFL; COMPLEX_NORM_ZERO] THEN + FIRST_ASSUM(MP_TAC o MATCH_MP POLY_DECOMPOSE) THEN + DISCH_THEN(X_CHOOSE_THEN `k:num` (X_CHOOSE_THEN `a:complex` + (X_CHOOSE_THEN `s:complex list` MP_TAC))) THEN + DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN + DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN + DISCH_THEN(CONJUNCTS_THEN2 (SUBST_ALL_TAC o SYM) MP_TAC) THEN + DISCH_THEN(fun th -> ONCE_REWRITE_TAC[th]) THEN ASM_REWRITE_TAC[] THEN + ASM_CASES_TAC `k + 1 = n` THENL + [UNDISCH_TAC `LENGTH(s:complex list) + k + 1 = n` THEN + ASM_REWRITE_TAC[ARITH_RULE `(s + m = m) <=> (s = 0)`; LENGTH_EQ_NIL] THEN + REWRITE_TAC[LENGTH_EQ_NIL] THEN DISCH_THEN SUBST1_TAC THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + ONCE_REWRITE_TAC[COMPLEX_MUL_SYM] THEN + MATCH_MP_TAC REDUCE_POLY_SIMPLE THEN ASM_REWRITE_TAC[] THEN + MAP_EVERY UNDISCH_TAC [`k + 1 = n`; `2 <= n`] THEN ARITH_TAC; ALL_TAC] THEN + FIRST_X_ASSUM(MP_TAC o SPEC `k + 1`) THEN ANTS_TAC THENL + [UNDISCH_TAC `~(k + 1 = n)` THEN + UNDISCH_TAC `LENGTH(s:complex list) + k + 1 = n` THEN ARITH_TAC; + ALL_TAC] THEN + DISCH_THEN(MP_TAC o SPEC + `CONS (Cx(&1)) (APPEND (REPLICATE (k - 1) (Cx(&0))) [a])`) THEN + ANTS_TAC THENL + [CONJ_TAC THENL + [REWRITE_TAC[LENGTH; LENGTH_APPEND; LENGTH_REPLICATE] THEN + UNDISCH_TAC `~(k = 0)` THEN ARITH_TAC; ALL_TAC] THEN + REWRITE_TAC[constant; POLY_REPLICATE_APPEND; poly] THEN + REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + DISCH_THEN(MP_TAC o SPECL [`Cx(&0)`; `Cx(&1)`]) THEN + REWRITE_TAC[COMPLEX_MUL_LID; COMPLEX_MUL_LZERO; COMPLEX_ADD_RID] THEN + ASM_REWRITE_TAC[COMPLEX_ENTIRE; COMPLEX_POW_ONE; SIMPLE_COMPLEX_ARITH + `(a = a + b) <=> (b = Cx(&0))`] THEN + REWRITE_TAC[CX_INJ; REAL_OF_NUM_EQ; ARITH_EQ]; ALL_TAC] THEN + REWRITE_TAC[constant; POLY_REPLICATE_APPEND; poly] THEN + REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + ONCE_REWRITE_TAC[AC COMPLEX_MUL_AC `a * b * c = (a * b) * c`] THEN + REWRITE_TAC[GSYM(CONJUNCT2 complex_pow)] THEN + ASM_SIMP_TAC[ARITH_RULE `~(k = 0) ==> (SUC(k - 1) = k)`] THEN + DISCH_THEN(X_CHOOSE_TAC `w:complex`) THEN + MP_TAC(SPECL [`s:complex list`; `norm(w)`] POLY_BOUND_EXISTS) THEN + DISCH_THEN(X_CHOOSE_THEN `m:real` STRIP_ASSUME_TAC) THEN + SUBGOAL_THEN `~(w = Cx(&0))` ASSUME_TAC THENL + [UNDISCH_TAC `Cx(&1) + w pow k * a = Cx(&0)` THEN + ONCE_REWRITE_TAC[TAUT `a ==> ~b <=> b ==> ~a`] THEN + DISCH_THEN SUBST1_TAC THEN + UNDISCH_TAC `~(k = 0)` THEN SPEC_TAC(`k:num`,`k:num`) THEN + INDUCT_TAC THEN REWRITE_TAC[complex_pow; COMPLEX_MUL_LZERO] THEN + REWRITE_TAC[COMPLEX_ADD_RID; CX_INJ; REAL_OF_NUM_EQ; ARITH_EQ]; + ALL_TAC] THEN + MP_TAC(SPECL [`&1`; `inv(norm(w) pow (k + 1) * m)`] REAL_DOWN2) THEN + ASM_SIMP_TAC[REAL_LT_01; REAL_LT_INV_EQ; REAL_LT_MUL; REAL_POW_LT; + COMPLEX_NORM_NZ] THEN + DISCH_THEN(X_CHOOSE_THEN `t:real` STRIP_ASSUME_TAC) THEN + EXISTS_TAC `Cx(t) * w` THEN REWRITE_TAC[COMPLEX_POW_MUL] THEN + REWRITE_TAC[COMPLEX_ADD_LDISTRIB; GSYM COMPLEX_MUL_ASSOC] THEN + FIRST_ASSUM(SUBST1_TAC o MATCH_MP (SIMPLE_COMPLEX_ARITH + `(a + w = Cx(&0)) ==> (w = --a)`)) THEN + REWRITE_TAC[GSYM CX_NEG; GSYM CX_POW; GSYM CX_MUL] THEN + REWRITE_TAC[COMPLEX_ADD_ASSOC; GSYM CX_ADD] THEN + MATCH_MP_TAC REAL_LET_TRANS THEN + EXISTS_TAC `norm(Cx(&1 + t pow k * -- &1)) + + norm(Cx(t pow k) * w pow k * Cx t * w * poly s (Cx t * w))` THEN + REWRITE_TAC[COMPLEX_NORM_TRIANGLE] THEN REWRITE_TAC[COMPLEX_NORM_CX] THEN + MATCH_MP_TAC(REAL_ARITH + `&0 <= x /\ x < t /\ t <= &1 ==> abs(&1 + t * --(&1)) + x < &1`) THEN + REWRITE_TAC[COMPLEX_NORM_POS] THEN + ASM_SIMP_TAC[REAL_POW_1_LE; REAL_LT_IMP_LE] THEN + ONCE_REWRITE_TAC[COMPLEX_NORM_MUL] THEN REWRITE_TAC[COMPLEX_NORM_CX] THEN + ASM_SIMP_TAC[real_abs; REAL_LT_IMP_LE; REAL_POW_LE] THEN + GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_RID] THEN + MATCH_MP_TAC REAL_LT_LMUL THEN ASM_SIMP_TAC[REAL_POW_LT] THEN + ONCE_REWRITE_TAC[AC COMPLEX_MUL_AC `a * b * c * d = b * (c * a) * d`] THEN + REWRITE_TAC[GSYM(CONJUNCT2 complex_pow)] THEN + REWRITE_TAC[COMPLEX_NORM_MUL; ADD1; COMPLEX_NORM_CX] THEN + MATCH_MP_TAC REAL_LET_TRANS THEN + EXISTS_TAC `abs t * norm(w pow (k + 1)) * m` THEN CONJ_TAC THENL + [MATCH_MP_TAC REAL_LE_LMUL THEN REWRITE_TAC[REAL_ABS_POS] THEN + MATCH_MP_TAC REAL_LE_LMUL THEN REWRITE_TAC[COMPLEX_NORM_POS] THEN + FIRST_ASSUM MATCH_MP_TAC THEN REWRITE_TAC[COMPLEX_NORM_MUL] THEN + GEN_REWRITE_TAC RAND_CONV [GSYM REAL_MUL_LID] THEN + MATCH_MP_TAC REAL_LE_RMUL THEN REWRITE_TAC[COMPLEX_NORM_POS] THEN + ASM_SIMP_TAC[COMPLEX_NORM_CX; REAL_ARITH + `&0 < t /\ t < &1 ==> abs(t) <= &1`]; ALL_TAC] THEN + ASM_SIMP_TAC[GSYM REAL_LT_RDIV_EQ; REAL_LT_MUL; COMPLEX_NORM_POW; + REAL_POW_LT; COMPLEX_NORM_NZ] THEN + ASM_SIMP_TAC[real_div; REAL_MUL_LID; + REAL_ARITH `&0 < t /\ t < x ==> abs(t) < x`]);; + +(* ------------------------------------------------------------------------- *) +(* Alternative version with a syntactic notion of constant polynomial. *) +(* ------------------------------------------------------------------------- *) + +let FUNDAMENTAL_THEOREM_OF_ALGEBRA_ALT = prove + (`!p. ~(?a l. ~(a = Cx(&0)) /\ ALL (\b. b = Cx(&0)) l /\ (p = CONS a l)) + ==> ?z. poly p z = Cx(&0)`, + LIST_INDUCT_TAC THEN REWRITE_TAC[poly; CONS_11] THEN + POP_ASSUM_LIST(K ALL_TAC) THEN + ONCE_REWRITE_TAC[AC CONJ_ACI `a /\ b /\ c /\ d <=> c /\ d /\ a /\ b`] THEN + REWRITE_TAC[RIGHT_EXISTS_AND_THM; UNWIND_THM1] THEN + ASM_CASES_TAC `h = Cx(&0)` THEN ASM_REWRITE_TAC[COMPLEX_ADD_LID] THENL + [EXISTS_TAC `Cx(&0)` THEN ASM_REWRITE_TAC[COMPLEX_MUL_LZERO]; ALL_TAC] THEN + DISCH_TAC THEN REWRITE_TAC[GSYM(CONJUNCT2 poly)] THEN + MATCH_MP_TAC FUNDAMENTAL_THEOREM_OF_ALGEBRA THEN + UNDISCH_TAC `~ALL (\b. b = Cx (&0)) t` THEN + REWRITE_TAC[TAUT `~b ==> ~a <=> a ==> b`] THEN POP_ASSUM(K ALL_TAC) THEN + REWRITE_TAC[constant; poly; REAL_EQ_LADD] THEN + DISCH_THEN(MP_TAC o SPEC `Cx(&0)` o ONCE_REWRITE_RULE[SWAP_FORALL_THM]) THEN + REWRITE_TAC[COMPLEX_MUL_LZERO; COMPLEX_EQ_ADD_LCANCEL] THEN + REWRITE_TAC[COMPLEX_ENTIRE; TAUT `a \/ b <=> ~a ==> b`] THEN + SPEC_TAC(`t:complex list`,`p:complex list`) THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[ALL] THEN + ASM_CASES_TAC `h = Cx(&0)` THEN + ASM_SIMP_TAC[poly; COMPLEX_ADD_LID; COMPLEX_ENTIRE] THEN + MP_TAC(SPECL [`t:complex list`; `&1`] POLY_BOUND_EXISTS) THEN + DISCH_THEN(X_CHOOSE_THEN `m:real` STRIP_ASSUME_TAC) THEN + MP_TAC(SPECL [`norm(h) / m`; `&1`] REAL_DOWN2) THEN + ASM_SIMP_TAC[REAL_LT_01; REAL_LT_DIV; COMPLEX_NORM_NZ] THEN + DISCH_THEN(X_CHOOSE_THEN `x:real` STRIP_ASSUME_TAC) THEN + DISCH_THEN(MP_TAC o SPEC `Cx(x)`) THEN + ASM_SIMP_TAC[CX_INJ; REAL_LT_IMP_NZ] THEN + REWRITE_TAC[SIMPLE_COMPLEX_ARITH `(x + y = Cx(&0)) <=> (y = --x)`] THEN + DISCH_THEN(MP_TAC o AP_TERM `norm`) THEN REWRITE_TAC[COMPLEX_NORM_NEG] THEN + MATCH_MP_TAC(REAL_ARITH `abs(a) < abs(b) ==> ~(a = b)`) THEN + REWRITE_TAC[real_abs; COMPLEX_NORM_POS] THEN + REWRITE_TAC[COMPLEX_NORM_MUL; COMPLEX_NORM_CX] THEN + MATCH_MP_TAC REAL_LTE_TRANS THEN EXISTS_TAC `norm(h) / m * m` THEN + CONJ_TAC THENL + [ALL_TAC; ASM_SIMP_TAC[REAL_LE_REFL; REAL_DIV_RMUL; REAL_LT_IMP_NZ]] THEN + MATCH_MP_TAC REAL_LET_TRANS THEN EXISTS_TAC `abs(x) * m` THEN + ASM_SIMP_TAC[REAL_LT_RMUL; REAL_ARITH `&0 < x /\ x < a ==> abs(x) < a`] THEN + ASM_MESON_TAC[REAL_LE_LMUL; REAL_ABS_POS; COMPLEX_NORM_CX; + REAL_ARITH `&0 < x /\ x < &1 ==> abs(x) <= &1`]);; diff --git a/Complex/grobner_examples.ml b/Complex/grobner_examples.ml new file mode 100644 index 0000000..468e64d --- /dev/null +++ b/Complex/grobner_examples.ml @@ -0,0 +1,589 @@ +(* ========================================================================= *) +(* Examples of using the Grobner basis procedure. *) +(* ========================================================================= *) + +time COMPLEX_ARITH + `!a b c. + (a * x pow 2 + b * x + c = Cx(&0)) /\ + (a * y pow 2 + b * y + c = Cx(&0)) /\ + ~(x = y) + ==> (a * (x + y) + b = Cx(&0))`;; + +time COMPLEX_ARITH + `!a b c. + (a * x pow 2 + b * x + c = Cx(&0)) /\ + (Cx(&2) * a * y pow 2 + Cx(&2) * b * y + Cx(&2) * c = Cx(&0)) /\ + ~(x = y) + ==> (a * (x + y) + b = Cx(&0))`;; + +(* ------------------------------------------------------------------------- *) +(* Another example. *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `~((y_1 = Cx(&2) * y_3) /\ + (y_2 = Cx(&2) * y_4) /\ + (y_1 * y_3 = y_2 * y_4) /\ + ((y_1 pow 2 - y_2 pow 2) * z = Cx(&1)))`;; + +time COMPLEX_ARITH + `!y_1 y_2 y_3 y_4. + (y_1 = Cx(&2) * y_3) /\ + (y_2 = Cx(&2) * y_4) /\ + (y_1 * y_3 = y_2 * y_4) + ==> (y_1 pow 2 = y_2 pow 2)`;; + +(* ------------------------------------------------------------------------- *) +(* Angle at centre vs. angle at circumference. *) +(* Formulation from "Real quantifier elimination in practice" paper. *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `~((c pow 2 = a pow 2 + b pow 2) /\ + (c pow 2 = x0 pow 2 + (y0 - b) pow 2) /\ + (y0 * t1 = a + x0) /\ + (y0 * t2 = a - x0) /\ + ((Cx(&1) - t1 * t2) * t = t1 + t2) /\ + (u * (b * t - a) = Cx(&1)) /\ + (v1 * a + v2 * x0 + v3 * y0 = Cx(&1)))`;; + +time COMPLEX_ARITH + `(c pow 2 = a pow 2 + b pow 2) /\ + (c pow 2 = x0 pow 2 + (y0 - b) pow 2) /\ + (y0 * t1 = a + x0) /\ + (y0 * t2 = a - x0) /\ + ((Cx(&1) - t1 * t2) * t = t1 + t2) /\ + (~(a = Cx(&0)) \/ ~(x0 = Cx(&0)) \/ ~(y0 = Cx(&0))) + ==> (b * t = a)`;; + +time COMPLEX_ARITH + `(c pow 2 = a pow 2 + b pow 2) /\ + (c pow 2 = x0 pow 2 + (y0 - b) pow 2) /\ + (y0 * t1 = a + x0) /\ + (y0 * t2 = a - x0) /\ + ((Cx(&1) - t1 * t2) * t = t1 + t2) /\ + (~(a = Cx(&0)) /\ ~(x0 = Cx(&0)) /\ ~(y0 = Cx(&0))) + ==> (b * t = a)`;; + +(* ------------------------------------------------------------------------- *) +(* Another example (note we rule out points 1, 2 or 3 coinciding). *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `((x1 - x0) pow 2 + (y1 - y0) pow 2 = + (x2 - x0) pow 2 + (y2 - y0) pow 2) /\ + ((x2 - x0) pow 2 + (y2 - y0) pow 2 = + (x3 - x0) pow 2 + (y3 - y0) pow 2) /\ + ((x1 - x0') pow 2 + (y1 - y0') pow 2 = + (x2 - x0') pow 2 + (y2 - y0') pow 2) /\ + ((x2 - x0') pow 2 + (y2 - y0') pow 2 = + (x3 - x0') pow 2 + (y3 - y0') pow 2) /\ + (a12 * (x1 - x2) + b12 * (y1 - y2) = Cx(&1)) /\ + (a13 * (x1 - x3) + b13 * (y1 - y3) = Cx(&1)) /\ + (a23 * (x2 - x3) + b23 * (y2 - y3) = Cx(&1)) /\ + ~((x1 - x0) pow 2 + (y1 - y0) pow 2 = Cx(&0)) + ==> (x0' = x0) /\ (y0' = y0)`;; + +time COMPLEX_ARITH + `~(((x1 - x0) pow 2 + (y1 - y0) pow 2 = + (x2 - x0) pow 2 + (y2 - y0) pow 2) /\ + ((x2 - x0) pow 2 + (y2 - y0) pow 2 = + (x3 - x0) pow 2 + (y3 - y0) pow 2) /\ + ((x1 - x0') pow 2 + (y1 - y0') pow 2 = + (x2 - x0') pow 2 + (y2 - y0') pow 2) /\ + ((x2 - x0') pow 2 + (y2 - y0') pow 2 = + (x3 - x0') pow 2 + (y3 - y0') pow 2) /\ + (a12 * (x1 - x2) + b12 * (y1 - y2) = Cx(&1)) /\ + (a13 * (x1 - x3) + b13 * (y1 - y3) = Cx(&1)) /\ + (a23 * (x2 - x3) + b23 * (y2 - y3) = Cx(&1)) /\ + (z * (x0' - x0) = Cx(&1)) /\ + (z' * (y0' - y0) = Cx(&1)) /\ + (z'' * ((x1 - x0) pow 2 + (y1 - y0) pow 2) = Cx(&1)) /\ + (z''' * ((x1 - x09) pow 2 + (y1 - y09) pow 2) = Cx(&1)))`;; + +(* ------------------------------------------------------------------------- *) +(* These are pure algebraic simplification. *) +(* ------------------------------------------------------------------------- *) + +let LAGRANGE_4 = time COMPLEX_ARITH + `(((x1 pow 2) + (x2 pow 2) + (x3 pow 2) + (x4 pow 2)) * + ((y1 pow 2) + (y2 pow 2) + (y3 pow 2) + (y4 pow 2))) = + ((((((x1*y1) - (x2*y2)) - (x3*y3)) - (x4*y4)) pow 2) + + (((((x1*y2) + (x2*y1)) + (x3*y4)) - (x4*y3)) pow 2) + + (((((x1*y3) - (x2*y4)) + (x3*y1)) + (x4*y2)) pow 2) + + (((((x1*y4) + (x2*y3)) - (x3*y2)) + (x4*y1)) pow 2))`;; + +let LAGRANGE_8 = time COMPLEX_ARITH + `((p1 pow 2 + q1 pow 2 + r1 pow 2 + s1 pow 2 + t1 pow 2 + u1 pow 2 + v1 pow 2 + w1 pow 2) * + (p2 pow 2 + q2 pow 2 + r2 pow 2 + s2 pow 2 + t2 pow 2 + u2 pow 2 + v2 pow 2 + w2 pow 2)) = + ((p1 * p2 - q1 * q2 - r1 * r2 - s1 * s2 - t1 * t2 - u1 * u2 - v1 * v2 - w1* w2) pow 2 + + (p1 * q2 + q1 * p2 + r1 * s2 - s1 * r2 + t1 * u2 - u1 * t2 - v1 * w2 + w1* v2) pow 2 + + (p1 * r2 - q1 * s2 + r1 * p2 + s1 * q2 + t1 * v2 + u1 * w2 - v1 * t2 - w1* u2) pow 2 + + (p1 * s2 + q1 * r2 - r1 * q2 + s1 * p2 + t1 * w2 - u1 * v2 + v1 * u2 - w1* t2) pow 2 + + (p1 * t2 - q1 * u2 - r1 * v2 - s1 * w2 + t1 * p2 + u1 * q2 + v1 * r2 + w1* s2) pow 2 + + (p1 * u2 + q1 * t2 - r1 * w2 + s1 * v2 - t1 * q2 + u1 * p2 - v1 * s2 + w1* r2) pow 2 + + (p1 * v2 + q1 * w2 + r1 * t2 - s1 * u2 - t1 * r2 + u1 * s2 + v1 * p2 - w1* q2) pow 2 + + (p1 * w2 - q1 * v2 + r1 * u2 + s1 * t2 - t1 * s2 - u1 * r2 + v1 * q2 + w1* p2) pow 2)`;; + +let LIOUVILLE = time COMPLEX_ARITH + `((x1 pow 2 + x2 pow 2 + x3 pow 2 + x4 pow 2) pow 2) = + (Cx(&1 / &6) * ((x1 + x2) pow 4 + (x1 + x3) pow 4 + (x1 + x4) pow 4 + + (x2 + x3) pow 4 + (x2 + x4) pow 4 + (x3 + x4) pow 4) + + Cx(&1 / &6) * ((x1 - x2) pow 4 + (x1 - x3) pow 4 + (x1 - x4) pow 4 + + (x2 - x3) pow 4 + (x2 - x4) pow 4 + (x3 - x4) pow 4))`;; + +let FLECK = time COMPLEX_ARITH + `((x1 pow 2 + x2 pow 2 + x3 pow 2 + x4 pow 2) pow 3) = + (Cx(&1 / &60) * ((x1 + x2 + x3) pow 6 + (x1 + x2 - x3) pow 6 + + (x1 - x2 + x3) pow 6 + (x1 - x2 - x3) pow 6 + + (x1 + x2 + x4) pow 6 + (x1 + x2 - x4) pow 6 + + (x1 - x2 + x4) pow 6 + (x1 - x2 - x4) pow 6 + + (x1 + x3 + x4) pow 6 + (x1 + x3 - x4) pow 6 + + (x1 - x3 + x4) pow 6 + (x1 - x3 - x4) pow 6 + + (x2 + x3 + x4) pow 6 + (x2 + x3 - x4) pow 6 + + (x2 - x3 + x4) pow 6 + (x2 - x3 - x4) pow 6) + + Cx(&1 / &30) * ((x1 + x2) pow 6 + (x1 - x2) pow 6 + + (x1 + x3) pow 6 + (x1 - x3) pow 6 + + (x1 + x4) pow 6 + (x1 - x4) pow 6 + + (x2 + x3) pow 6 + (x2 - x3) pow 6 + + (x2 + x4) pow 6 + (x2 - x4) pow 6 + + (x3 + x4) pow 6 + (x3 - x4) pow 6) + + Cx(&3 / &5) * (x1 pow 6 + x2 pow 6 + x3 pow 6 + x4 pow 6))`;; + +let HURWITZ = time COMPLEX_ARITH + `!x1 x2 x3 x4. + (x1 pow 2 + x2 pow 2 + x3 pow 2 + x4 pow 2) pow 4 = + Cx(&1 / &840) * ((x1 + x2 + x3 + x4) pow 8 + + (x1 + x2 + x3 - x4) pow 8 + + (x1 + x2 - x3 + x4) pow 8 + + (x1 + x2 - x3 - x4) pow 8 + + (x1 - x2 + x3 + x4) pow 8 + + (x1 - x2 + x3 - x4) pow 8 + + (x1 - x2 - x3 + x4) pow 8 + + (x1 - x2 - x3 - x4) pow 8) + + Cx(&1 / &5040) * ((Cx(&2) * x1 + x2 + x3) pow 8 + + (Cx(&2) * x1 + x2 - x3) pow 8 + + (Cx(&2) * x1 - x2 + x3) pow 8 + + (Cx(&2) * x1 - x2 - x3) pow 8 + + (Cx(&2) * x1 + x2 + x4) pow 8 + + (Cx(&2) * x1 + x2 - x4) pow 8 + + (Cx(&2) * x1 - x2 + x4) pow 8 + + (Cx(&2) * x1 - x2 - x4) pow 8 + + (Cx(&2) * x1 + x3 + x4) pow 8 + + (Cx(&2) * x1 + x3 - x4) pow 8 + + (Cx(&2) * x1 - x3 + x4) pow 8 + + (Cx(&2) * x1 - x3 - x4) pow 8 + + (Cx(&2) * x2 + x3 + x4) pow 8 + + (Cx(&2) * x2 + x3 - x4) pow 8 + + (Cx(&2) * x2 - x3 + x4) pow 8 + + (Cx(&2) * x2 - x3 - x4) pow 8 + + (x1 + Cx(&2) * x2 + x3) pow 8 + + (x1 + Cx(&2) * x2 - x3) pow 8 + + (x1 - Cx(&2) * x2 + x3) pow 8 + + (x1 - Cx(&2) * x2 - x3) pow 8 + + (x1 + Cx(&2) * x2 + x4) pow 8 + + (x1 + Cx(&2) * x2 - x4) pow 8 + + (x1 - Cx(&2) * x2 + x4) pow 8 + + (x1 - Cx(&2) * x2 - x4) pow 8 + + (x1 + Cx(&2) * x3 + x4) pow 8 + + (x1 + Cx(&2) * x3 - x4) pow 8 + + (x1 - Cx(&2) * x3 + x4) pow 8 + + (x1 - Cx(&2) * x3 - x4) pow 8 + + (x2 + Cx(&2) * x3 + x4) pow 8 + + (x2 + Cx(&2) * x3 - x4) pow 8 + + (x2 - Cx(&2) * x3 + x4) pow 8 + + (x2 - Cx(&2) * x3 - x4) pow 8 + + (x1 + x2 + Cx(&2) * x3) pow 8 + + (x1 + x2 - Cx(&2) * x3) pow 8 + + (x1 - x2 + Cx(&2) * x3) pow 8 + + (x1 - x2 - Cx(&2) * x3) pow 8 + + (x1 + x2 + Cx(&2) * x4) pow 8 + + (x1 + x2 - Cx(&2) * x4) pow 8 + + (x1 - x2 + Cx(&2) * x4) pow 8 + + (x1 - x2 - Cx(&2) * x4) pow 8 + + (x1 + x3 + Cx(&2) * x4) pow 8 + + (x1 + x3 - Cx(&2) * x4) pow 8 + + (x1 - x3 + Cx(&2) * x4) pow 8 + + (x1 - x3 - Cx(&2) * x4) pow 8 + + (x2 + x3 + Cx(&2) * x4) pow 8 + + (x2 + x3 - Cx(&2) * x4) pow 8 + + (x2 - x3 + Cx(&2) * x4) pow 8 + + (x2 - x3 - Cx(&2) * x4) pow 8) + + Cx(&1 / &84) * ((x1 + x2) pow 8 + (x1 - x2) pow 8 + + (x1 + x3) pow 8 + (x1 - x3) pow 8 + + (x1 + x4) pow 8 + (x1 - x4) pow 8 + + (x2 + x3) pow 8 + (x2 - x3) pow 8 + + (x2 + x4) pow 8 + (x2 - x4) pow 8 + + (x3 + x4) pow 8 + (x3 - x4) pow 8) + + Cx(&1 / &840) * ((Cx(&2) * x1) pow 8 + (Cx(&2) * x2) pow 8 + + (Cx(&2) * x3) pow 8 + (Cx(&2) * x4) pow 8)`;; + +let SCHUR = time COMPLEX_ARITH + `Cx(&22680) * (x1 pow 2 + x2 pow 2 + x3 pow 2 + x4 pow 2) pow 5 = + Cx(&9) * ((Cx(&2) * x1) pow 10 + + (Cx(&2) * x2) pow 10 + + (Cx(&2) * x3) pow 10 + + (Cx(&2) * x4) pow 10) + + Cx(&180) * ((x1 + x2) pow 10 + (x1 - x2) pow 10 + + (x1 + x3) pow 10 + (x1 - x3) pow 10 + + (x1 + x4) pow 10 + (x1 - x4) pow 10 + + (x2 + x3) pow 10 + (x2 - x3) pow 10 + + (x2 + x4) pow 10 + (x2 - x4) pow 10 + + (x3 + x4) pow 10 + (x3 - x4) pow 10) + + ((Cx(&2) * x1 + x2 + x3) pow 10 + + (Cx(&2) * x1 + x2 - x3) pow 10 + + (Cx(&2) * x1 - x2 + x3) pow 10 + + (Cx(&2) * x1 - x2 - x3) pow 10 + + (Cx(&2) * x1 + x2 + x4) pow 10 + + (Cx(&2) * x1 + x2 - x4) pow 10 + + (Cx(&2) * x1 - x2 + x4) pow 10 + + (Cx(&2) * x1 - x2 - x4) pow 10 + + (Cx(&2) * x1 + x3 + x4) pow 10 + + (Cx(&2) * x1 + x3 - x4) pow 10 + + (Cx(&2) * x1 - x3 + x4) pow 10 + + (Cx(&2) * x1 - x3 - x4) pow 10 + + (Cx(&2) * x2 + x3 + x4) pow 10 + + (Cx(&2) * x2 + x3 - x4) pow 10 + + (Cx(&2) * x2 - x3 + x4) pow 10 + + (Cx(&2) * x2 - x3 - x4) pow 10 + + (x1 + Cx(&2) * x2 + x3) pow 10 + + (x1 + Cx(&2) * x2 - x3) pow 10 + + (x1 - Cx(&2) * x2 + x3) pow 10 + + (x1 - Cx(&2) * x2 - x3) pow 10 + + (x1 + Cx(&2) * x2 + x4) pow 10 + + (x1 + Cx(&2) * x2 - x4) pow 10 + + (x1 - Cx(&2) * x2 + x4) pow 10 + + (x1 - Cx(&2) * x2 - x4) pow 10 + + (x1 + Cx(&2) * x3 + x4) pow 10 + + (x1 + Cx(&2) * x3 - x4) pow 10 + + (x1 - Cx(&2) * x3 + x4) pow 10 + + (x1 - Cx(&2) * x3 - x4) pow 10 + + (x2 + Cx(&2) * x3 + x4) pow 10 + + (x2 + Cx(&2) * x3 - x4) pow 10 + + (x2 - Cx(&2) * x3 + x4) pow 10 + + (x2 - Cx(&2) * x3 - x4) pow 10 + + (x1 + x2 + Cx(&2) * x3) pow 10 + + (x1 + x2 - Cx(&2) * x3) pow 10 + + (x1 - x2 + Cx(&2) * x3) pow 10 + + (x1 - x2 - Cx(&2) * x3) pow 10 + + (x1 + x2 + Cx(&2) * x4) pow 10 + + (x1 + x2 - Cx(&2) * x4) pow 10 + + (x1 - x2 + Cx(&2) * x4) pow 10 + + (x1 - x2 - Cx(&2) * x4) pow 10 + + (x1 + x3 + Cx(&2) * x4) pow 10 + + (x1 + x3 - Cx(&2) * x4) pow 10 + + (x1 - x3 + Cx(&2) * x4) pow 10 + + (x1 - x3 - Cx(&2) * x4) pow 10 + + (x2 + x3 + Cx(&2) * x4) pow 10 + + (x2 + x3 - Cx(&2) * x4) pow 10 + + (x2 - x3 + Cx(&2) * x4) pow 10 + + (x2 - x3 - Cx(&2) * x4) pow 10) + + Cx(&9) * ((x1 + x2 + x3 + x4) pow 10 + + (x1 + x2 + x3 - x4) pow 10 + + (x1 + x2 - x3 + x4) pow 10 + + (x1 + x2 - x3 - x4) pow 10 + + (x1 - x2 + x3 + x4) pow 10 + + (x1 - x2 + x3 - x4) pow 10 + + (x1 - x2 - x3 + x4) pow 10 + + (x1 - x2 - x3 - x4) pow 10)`;; + +(* ------------------------------------------------------------------------- *) +(* Intersection of diagonals of a parallelogram is their midpoint. *) +(* Kapur "...Dixon resultants, Groebner Bases, and Characteristic Sets", 3.1 *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `(x1 = u3) /\ + (x1 * (u2 - u1) = x2 * u3) /\ + (x4 * (x2 - u1) = x1 * (x3 - u1)) /\ + (x3 * u3 = x4 * u2) /\ + ~(u1 = Cx(&0)) /\ + ~(u3 = Cx(&0)) + ==> (x3 pow 2 + x4 pow 2 = (u2 - x3) pow 2 + (u3 - x4) pow 2)`;; + +(* ------------------------------------------------------------------------- *) +(* Chou's formulation of same property. *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `(u1 * x1 - u1 * u3 = Cx(&0)) /\ + (u3 * x2 - (u2 - u1) * x1 = Cx(&0)) /\ + (x1 * x4 - (x2 - u1) * x3 - u1 * x1 = Cx(&0)) /\ + (u3 * x4 - u2 * x3 = Cx(&0)) /\ + ~(u1 = Cx(&0)) /\ + ~(u3 = Cx(&0)) + ==> (Cx(&2) * u2 * x4 + Cx(&2) * u3 * x3 - u3 pow 2 - u2 pow 2 = Cx(&0))`;; + +(* ------------------------------------------------------------------------- *) +(* Perpendicular lines property; from Kapur's earlier paper. *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `(y1 * y3 + x1 * x3 = Cx(&0)) /\ + (y3 * (y2 - y3) + (x2 - x3) * x3 = Cx(&0)) /\ + ~(x3 = Cx(&0)) /\ + ~(y3 = Cx(&0)) + ==> (y1 * (x2 - x3) = x1 * (y2 - y3))`;; + +(* ------------------------------------------------------------------------- *) +(* Simson's theorem (Chou, p7). *) +(* ------------------------------------------------------------------------- *) + +time COMPLEX_ARITH + `(Cx(&2) * u2 * x2 + Cx(&2) * u3 * x1 - u3 pow 2 - u2 pow 2 = Cx(&0)) /\ + (Cx(&2) * u1 * x2 - u1 pow 2 = Cx(&0)) /\ + (--(x3 pow 2) + Cx(&2) * x2 * x3 + Cx(&2) * u4 * x1 - u4 pow 2 = Cx(&0)) /\ + (u3 * x5 + (--u2 + u1) * x4 - u1 * u3 = Cx(&0)) /\ + ((u2 - u1) * x5 + u3 * x4 + (--u2 + u1) * x3 - u3 * u4 = Cx(&0)) /\ + (u3 * x7 - u2 * x6 = Cx(&0)) /\ + (u2 * x7 + u3 * x6 - u2 * x3 - u3 * u4 = Cx(&0)) /\ + ~(Cx(&4) * u1 * u3 = Cx(&0)) /\ + ~(Cx(&2) * u1 = Cx(&0)) /\ + ~(--(u3 pow 2) - u2 pow 2 + Cx(&2) * u1 * u2 - u1 pow 2 = Cx(&0)) /\ + ~(u3 = Cx(&0)) /\ + ~(--(u3 pow 2) - u2 pow 2 = Cx(&0)) /\ + ~(u2 = Cx(&0)) + ==> (x4 * x7 + (--x5 + x3) * x6 - x3 * x4 = Cx(&0))`;; + +(* ------------------------------------------------------------------------- *) +(* Determinants from Coq convex hull paper (some require reals or order). *) +(* ------------------------------------------------------------------------- *) + +let det3 = new_definition + `det3(a11,a12,a13, + a21,a22,a23, + a31,a32,a33) = + a11 * (a22 * a33 - a32 * a23) - + a12 * (a21 * a33 - a31 * a23) + + a13 * (a21 * a32 - a31 * a22)`;; + +let DET_TRANSPOSE = prove + (`det3(a1,b1,c1,a2,b2,c2,a3,b3,c3) = + det3(a1,a2,a3,b1,b2,b3,c1,c2,c3)`, + REWRITE_TAC[det3] THEN CONV_TAC(time COMPLEX_ARITH));; + +let sdet3 = new_definition + `sdet3(p,q,r) = det3(FST p,SND p,Cx(&1), + FST q,SND q,Cx(&1), + FST r,SND r,Cx(&1))`;; + +let SDET3_PERMUTE_1 = prove + (`sdet3(p,q,r) = sdet3(q,r,p)`, + REWRITE_TAC[sdet3; det3] THEN CONV_TAC(time COMPLEX_ARITH));; + +let SDET3_PERMUTE_2 = prove + (`sdet3(p,q,r) = --(sdet3(p,r,q))`, + REWRITE_TAC[sdet3; det3] THEN CONV_TAC(time COMPLEX_ARITH));; + +let SDET_SUM = prove + (`sdet3(p,q,r) - sdet3(t,q,r) - sdet3(p,t,r) - sdet3(p,q,t) = Cx(&0)`, + REWRITE_TAC[sdet3; det3] THEN CONV_TAC(time COMPLEX_ARITH));; + +let SDET_CRAMER = prove + (`sdet3(s,t,q) * sdet3(t,p,r) = sdet3(t,q,r) * sdet3(s,t,p) + + sdet3(t,p,q) * sdet3(s,t,r)`, + REWRITE_TAC[sdet3; det3] THEN CONV_TAC(time COMPLEX_ARITH));; + +let SDET_NZ = prove + (`!p q r. ~(sdet3(p,q,r) = Cx(&0)) + ==> ~(p = q) /\ ~(q = r) /\ ~(r = p)`, + REWRITE_TAC[FORALL_PAIR_THM; PAIR_EQ; sdet3; det3] THEN + CONV_TAC(time COMPLEX_ARITH));; + +let SDET_LINCOMB = prove + (`(FST p * sdet3(i,j,k) = + FST i * sdet3(j,k,p) + FST j * sdet3(k,i,p) + FST k * sdet3(i,j,p)) /\ + (SND p * sdet3(i,j,k) = + SND i * sdet3(j,k,p) + SND j * sdet3(k,i,p) + SND k * sdet3(i,j,p))`, + REWRITE_TAC[sdet3; det3] THEN CONV_TAC(time COMPLEX_ARITH));; + +(***** I'm not sure if this is true; there must be some + sufficient degenerate conditions.... + +let th = prove + (`~(~(xp pow 2 + yp pow 2 = Cx(&0)) /\ + ~(xq pow 2 + yq pow 2 = Cx(&0)) /\ + ~(xr pow 2 + yr pow 2 = Cx(&0)) /\ + (det3(xp,yp,Cx(&1), + xq,yq,Cx(&1), + xr,yr,Cx(&1)) = Cx(&0)) /\ + (det3(yp,xp pow 2 + yp pow 2,Cx(&1), + yq,xq pow 2 + yq pow 2,Cx(&1), + yr,xr pow 2 + yr pow 2,Cx(&1)) = Cx(&0)) /\ + (det3(xp,xp pow 2 + yp pow 2,Cx(&1), + xq,xq pow 2 + yq pow 2,Cx(&1), + xr,xr pow 2 + yr pow 2,Cx(&1)) = Cx(&0)))`, + REWRITE_TAC[det3] THEN + CONV_TAC(time COMPLEX_ARITH));; + +***************) + +(* ------------------------------------------------------------------------- *) +(* Some geometry concepts (just "axiomatic" in this file). *) +(* ------------------------------------------------------------------------- *) + +prioritize_real();; + +let collinear = new_definition + `collinear (a:real#real) b c <=> + ((FST a - FST b) * (SND b - SND c) = + (SND a - SND b) * (FST b - FST c))`;; + +let parallel = new_definition + `parallel (a,b) (c,d) <=> + ((FST a - FST b) * (SND c - SND d) = (SND a - SND b) * (FST c - FST d))`;; + +let perpendicular = new_definition + `perpendicular (a,b) (c,d) <=> + ((FST a - FST b) * (FST c - FST d) + (SND a - SND b) * (SND c - SND d) = + &0)`;; + +let oncircle_with_diagonal = new_definition + `oncircle_with_diagonal a (b,c) = perpendicular (b,a) (c,a)`;; + +let length = new_definition + `length (a,b) = sqrt((FST a - FST b) pow 2 + (SND a - SND b) pow 2)`;; + +let lengths_eq = new_definition + `lengths_eq (a,b) (c,d) <=> + ((FST a - FST b) pow 2 + (SND a - SND b) pow 2 = + (FST c - FST d) pow 2 + (SND c - SND d) pow 2)`;; + +let is_midpoint = new_definition + `is_midpoint b (a,c) <=> + (&2 * FST b = FST a + FST c) /\ + (&2 * SND b = SND a + SND c)`;; + +(* ------------------------------------------------------------------------- *) +(* Chou isn't explicit about this. *) +(* ------------------------------------------------------------------------- *) + +let is_intersection = new_definition + `is_intersection p (a,b) (c,d) <=> + collinear a p b /\ collinear c p d`;; + +(* ------------------------------------------------------------------------- *) +(* This is used in some degenerate conditions. See Chou, p18. *) +(* ------------------------------------------------------------------------- *) + +let isotropic = new_definition + `isotropic (a,b) = perpendicular (a,b) (a,b)`;; + +(* ------------------------------------------------------------------------- *) +(* This increases degree, but sometimes makes complex assertion useful. *) +(* ------------------------------------------------------------------------- *) + +let distinctpairs = new_definition + `distinctpairs pprs <=> + ~(ITLIST (\(a,b) pr. ((FST a - FST b) pow 2 + (SND a - SND b) pow 2) * pr) + pprs (&1) = &0)`;; + +(* ------------------------------------------------------------------------- *) +(* Simple tactic to remove defined concepts and expand coordinates. *) +(* ------------------------------------------------------------------------- *) + +let (EXPAND_COORDS_TAC:tactic) = + let complex2_ty = `:real#real` in + fun (asl,w) -> + (let fvs = filter (fun v -> type_of v = complex2_ty) (frees w) in + MAP_EVERY (fun v -> SPEC_TAC(v,v)) fvs THEN + GEN_REWRITE_TAC DEPTH_CONV [FORALL_PAIR_THM; EXISTS_PAIR_THM] THEN + REPEAT GEN_TAC) (asl,w);; + +let PAIR_BETA_THM = prove + (`(\(x,y). P x y) (a,b) = P a b`, + CONV_TAC(LAND_CONV GEN_BETA_CONV) THEN REFL_TAC);; + +let GEOM_TAC = + EXPAND_COORDS_TAC THEN + GEN_REWRITE_TAC TOP_DEPTH_CONV + [collinear; parallel; perpendicular; oncircle_with_diagonal; + length; lengths_eq; is_midpoint; is_intersection; distinctpairs; + isotropic; ITLIST; PAIR_BETA_THM; BETA_THM; PAIR_EQ; FST; SND];; + +(* ------------------------------------------------------------------------- *) +(* Centroid (Chou, example 142). *) +(* ------------------------------------------------------------------------- *) + +let CENTROID = time prove + (`is_midpoint d (b,c) /\ + is_midpoint e (a,c) /\ + is_midpoint f (a,b) /\ + is_intersection m (b,e) (a,d) + ==> collinear c f m`, + GEOM_TAC THEN CONV_TAC GROBNER_REAL_ARITH);; + +(* ------------------------------------------------------------------------- *) +(* Gauss's theorem (Chou, example 15). *) +(* ------------------------------------------------------------------------- *) + +let GAUSS = time prove + (`collinear x a0 a3 /\ + collinear x a1 a2 /\ + collinear y a2 a3 /\ + collinear y a1 a0 /\ + is_midpoint m1 (a1,a3) /\ + is_midpoint m2 (a0,a2) /\ + is_midpoint m3 (x,y) + ==> collinear m1 m2 m3`, + GEOM_TAC THEN CONV_TAC GROBNER_REAL_ARITH);; + +(* ------------------------------------------------------------------------- *) +(* Simson's theorem (Chou, example 288). *) +(* ------------------------------------------------------------------------- *) + +(**** These are all hideously slow. At least the first one works. + I haven't had the patience to try the rest. + +let SIMSON = time prove + (`lengths_eq (O,a) (O,b) /\ + lengths_eq (O,a) (O,c) /\ + lengths_eq (d,O) (O,a) /\ + perpendicular (e,d) (b,c) /\ + collinear e b c /\ + perpendicular (f,d) (a,c) /\ + collinear f a c /\ + perpendicular (g,d) (a,b) /\ + collinear g a b /\ + ~(collinear a c b) /\ + ~(lengths_eq (a,b) (a,a)) /\ + ~(lengths_eq (a,c) (a,a)) /\ + ~(lengths_eq (b,c) (a,a)) + ==> collinear e f g`, + GEOM_TAC THEN CONV_TAC GROBNER_REAL_ARITH);; + +let SIMSON = time prove + (`lengths_eq (O,a) (O,b) /\ + lengths_eq (O,a) (O,c) /\ + lengths_eq (d,O) (O,a) /\ + perpendicular (e,d) (b,c) /\ + collinear e b c /\ + perpendicular (f,d) (a,c) /\ + collinear f a c /\ + perpendicular (g,d) (a,b) /\ + collinear g a b /\ + ~(a = b) /\ ~(a = c) /\ ~(a = d) /\ ~(b = c) /\ ~(b = d) /\ ~(c = d) + ==> collinear e f g`, + GEOM_TAC THEN CONV_TAC GROBNER_REAL_ARITH);; + +let SIMSON = time prove + (`lengths_eq (O,a) (O,b) /\ + lengths_eq (O,a) (O,c) /\ + lengths_eq (d,O) (O,a) /\ + perpendicular (e,d) (b,c) /\ + collinear e b c /\ + perpendicular (f,d) (a,c) /\ + collinear f a c /\ + perpendicular (g,d) (a,b) /\ + collinear g a b /\ + ~(collinear a c b) /\ + ~(isotropic (a,b)) /\ + ~(isotropic (a,c)) /\ + ~(isotropic (b,c)) /\ + ~(isotropic (a,d)) /\ + ~(isotropic (b,d)) /\ + ~(isotropic (c,d)) + ==> collinear e f g`, + GEOM_TAC THEN CONV_TAC GROBNER_REAL_ARITH);; + +****************) diff --git a/Complex/quelim.ml b/Complex/quelim.ml new file mode 100644 index 0000000..a628f29 --- /dev/null +++ b/Complex/quelim.ml @@ -0,0 +1,923 @@ +(* ========================================================================= *) +(* Naive quantifier elimination for complex numbers. *) +(* ========================================================================= *) + +needs "Complex/fundamental.ml";; + +let NULLSTELLENSATZ_LEMMA = prove + (`!n p q. (!x. (poly p x = Cx(&0)) ==> (poly q x = Cx(&0))) /\ + (degree p = n) /\ ~(n = 0) + ==> p divides (q exp n)`, + MATCH_MP_TAC num_WF THEN X_GEN_TAC `n:num` THEN DISCH_TAC THEN + MAP_EVERY X_GEN_TAC [`p:complex list`; `q:complex list`] THEN + ASM_CASES_TAC `?a. poly p a = Cx(&0)` THENL + [ALL_TAC; + DISCH_THEN(K ALL_TAC) THEN + FIRST_ASSUM(MP_TAC o MATCH_MP + (ONCE_REWRITE_RULE[TAUT `a ==> b <=> ~b ==> ~a`] + FUNDAMENTAL_THEOREM_OF_ALGEBRA_ALT)) THEN + REWRITE_TAC[LEFT_IMP_EXISTS_THM] THEN + MAP_EVERY X_GEN_TAC [`k:complex`; `zeros:complex list`] THEN + STRIP_TAC THEN REWRITE_TAC[divides] THEN + EXISTS_TAC `[inv(k)] ** q exp n` THEN + ASM_REWRITE_TAC[FUN_EQ_THM; POLY_MUL] THEN X_GEN_TAC `z:complex` THEN + ASM_SIMP_TAC[COMPLEX_MUL_ASSOC; COMPLEX_MUL_RINV; COMPLEX_MUL_LID; + poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; POLY_0]] THEN + FIRST_X_ASSUM(X_CHOOSE_THEN `a:complex` MP_TAC) THEN + DISCH_THEN(fun th -> ASSUME_TAC th THEN MP_TAC th) THEN + GEN_REWRITE_TAC LAND_CONV [ORDER_ROOT] THEN + ASM_CASES_TAC `poly p = poly []` THEN ASM_REWRITE_TAC[] THENL + [ASM_SIMP_TAC[DEGREE_ZERO] THEN MESON_TAC[]; ALL_TAC] THEN + STRIP_TAC THEN STRIP_TAC THEN + MP_TAC(SPECL [`p:complex list`; `a:complex`; `order a p`] ORDER) THEN + ASM_REWRITE_TAC[] THEN STRIP_TAC THEN + FIRST_ASSUM(MP_TAC o SPEC `a:complex` o MATCH_MP ORDER_DEGREE) THEN + ASM_REWRITE_TAC[] THEN DISCH_TAC THEN + FIRST_ASSUM(MP_TAC o SPEC `a:complex`) THEN + REWRITE_TAC[ASSUME `poly p a = Cx(&0)`] THEN + REWRITE_TAC[POLY_LINEAR_DIVIDES] THEN + ASM_CASES_TAC `q:complex list = []` THENL + [DISCH_TAC THEN MATCH_MP_TAC POLY_DIVIDES_ZERO THEN + UNDISCH_TAC `~(n = 0)` THEN SPEC_TAC(`n:num`,`n:num`) THEN + INDUCT_TAC THEN ASM_REWRITE_TAC[poly_exp] THEN DISCH_TAC THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL; COMPLEX_MUL_LZERO; poly]; + ALL_TAC] THEN + ASM_REWRITE_TAC[] THEN + DISCH_THEN(X_CHOOSE_THEN `r:complex list` SUBST_ALL_TAC) THEN + UNDISCH_TAC `[--a; Cx (&1)] exp (order a p) divides p` THEN + GEN_REWRITE_TAC LAND_CONV [divides] THEN + DISCH_THEN(X_CHOOSE_THEN `s:complex list` ASSUME_TAC) THEN + SUBGOAL_THEN `~(poly s = poly [])` ASSUME_TAC THENL + [DISCH_TAC THEN UNDISCH_TAC `~(poly p = poly [])` THEN + ASM_REWRITE_TAC[POLY_ENTIRE]; ALL_TAC] THEN + ASM_CASES_TAC `degree s = 0` THENL + [SUBGOAL_THEN `?k. ~(k = Cx(&0)) /\ (poly s = poly [k])` MP_TAC THENL + [EXISTS_TAC `LAST(normalize s)` THEN + ASM_SIMP_TAC[NORMAL_NORMALIZE; GSYM POLY_NORMALIZE_ZERO] THEN + GEN_REWRITE_TAC LAND_CONV [GSYM POLY_NORMALIZE] THEN + UNDISCH_TAC `degree s = 0` THEN + FIRST_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV + [POLY_NORMALIZE_ZERO]) THEN + REWRITE_TAC[degree] THEN + SPEC_TAC(`normalize s`,`s:complex list`) THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[NOT_CONS_NIL] THEN + REWRITE_TAC[LENGTH; PRE; poly; LAST] THEN + COND_CASES_TAC THEN ASM_REWRITE_TAC[] THEN + ASM_REWRITE_TAC[LENGTH_EQ_NIL]; ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `k:complex` STRIP_ASSUME_TAC) THEN + REWRITE_TAC[divides] THEN + EXISTS_TAC `[inv(k)] ** [--a; Cx (&1)] exp (n - order a p) ** r exp n` THEN + ASM_REWRITE_TAC[FUN_EQ_THM; POLY_MUL; POLY_EXP; COMPLEX_POW_MUL] THEN + X_GEN_TAC `z:complex` THEN + ONCE_REWRITE_TAC[AC COMPLEX_MUL_AC + `(a * b) * c * d * e = ((d * a) * (c * b)) * e`] THEN + AP_THM_TAC THEN AP_TERM_TAC THEN + REWRITE_TAC[GSYM COMPLEX_POW_ADD] THEN ASM_SIMP_TAC[SUB_ADD] THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID; COMPLEX_MUL_RID] THEN + ASM_SIMP_TAC[COMPLEX_MUL_LINV; COMPLEX_MUL_RID]; ALL_TAC] THEN + SUBGOAL_THEN `degree s < n` ASSUME_TAC THENL + [EXPAND_TAC "n" THEN + FIRST_ASSUM(SUBST1_TAC o MATCH_MP DEGREE_WELLDEF) THEN + REWRITE_TAC[LINEAR_POW_MUL_DEGREE] THEN + ASM_REWRITE_TAC[] THEN UNDISCH_TAC `~(order a p = 0)` THEN ARITH_TAC; + ALL_TAC] THEN + FIRST_X_ASSUM(MP_TAC o SPEC `degree s`) THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN(MP_TAC o SPECL [`s:complex list`; `r:complex list`]) THEN + ASM_REWRITE_TAC[] THEN ANTS_TAC THENL + [X_GEN_TAC `z:complex` THEN DISCH_TAC THEN + UNDISCH_TAC + `!x. (poly p x = Cx(&0)) ==> (poly([--a; Cx (&1)] ** r) x = Cx(&0))` THEN + DISCH_THEN(MP_TAC o SPEC `z:complex`) THEN ASM_REWRITE_TAC[] THEN + ASM_REWRITE_TAC[POLY_MUL; COMPLEX_MUL_RID] THEN + REWRITE_TAC[COMPLEX_ENTIRE] THEN + MATCH_MP_TAC(TAUT `~a ==> (a \/ b ==> b)`) THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + REWRITE_TAC[SIMPLE_COMPLEX_ARITH + `(--a + z * Cx(&1) = Cx(&0)) <=> (z = a)`] THEN + DISCH_THEN SUBST_ALL_TAC THEN + UNDISCH_TAC `poly s a = Cx (&0)` THEN + ASM_REWRITE_TAC[POLY_LINEAR_DIVIDES; DE_MORGAN_THM] THEN + CONJ_TAC THENL [ASM_MESON_TAC[]; ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `u:complex list` SUBST_ALL_TAC) THEN + UNDISCH_TAC `~([--a; Cx (&1)] exp SUC (order a p) divides p)` THEN + REWRITE_TAC[divides] THEN + EXISTS_TAC `u:complex list` THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[POLY_MUL; poly_exp; COMPLEX_MUL_AC; FUN_EQ_THM]; + ALL_TAC] THEN + REWRITE_TAC[divides] THEN + DISCH_THEN(X_CHOOSE_THEN `u:complex list` ASSUME_TAC) THEN + EXISTS_TAC + `u ** [--a; Cx(&1)] exp (n - order a p) ** r exp (n - degree s)` THEN + ASM_REWRITE_TAC[FUN_EQ_THM; POLY_MUL; POLY_EXP; COMPLEX_POW_MUL] THEN + X_GEN_TAC `z:complex` THEN + ONCE_REWRITE_TAC[AC COMPLEX_MUL_AC + `(ap * s) * u * anp * rns = (anp * ap) * rns * s * u`] THEN + REWRITE_TAC[GSYM COMPLEX_POW_ADD] THEN + ASM_SIMP_TAC[SUB_ADD] THEN AP_TERM_TAC THEN + GEN_REWRITE_TAC (RAND_CONV o RAND_CONV) [GSYM POLY_MUL] THEN + SUBST1_TAC(SYM(ASSUME `poly (r exp degree s) = poly (s ** u)`)) THEN + REWRITE_TAC[POLY_EXP; GSYM COMPLEX_POW_ADD] THEN + ASM_SIMP_TAC[SUB_ADD; LT_IMP_LE]);; + +let NULLSTELLENSATZ_UNIVARIATE = prove + (`!p q. (!x. (poly p x = Cx(&0)) ==> (poly q x = Cx(&0))) <=> + p divides (q exp (degree p)) \/ + ((poly p = poly []) /\ (poly q = poly []))`, + REPEAT GEN_TAC THEN ASM_CASES_TAC `poly p = poly []` THENL + [ASM_REWRITE_TAC[poly] THEN + FIRST_ASSUM(SUBST1_TAC o MATCH_MP DEGREE_WELLDEF) THEN + REWRITE_TAC[degree; normalize; LENGTH; ARITH; poly_exp] THEN + ASM_REWRITE_TAC[divides; FUN_EQ_THM; POLY_MUL; poly; + COMPLEX_MUL_LZERO; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + REWRITE_TAC[CX_INJ; REAL_OF_NUM_EQ; ARITH]; ALL_TAC] THEN + ASM_CASES_TAC `degree p = 0` THENL + [ALL_TAC; + MP_TAC(SPECL [`degree p`; `p:complex list`; `q:complex list`] + NULLSTELLENSATZ_LEMMA) THEN + ASM_REWRITE_TAC[] THEN DISCH_TAC THEN EQ_TAC THEN ASM_REWRITE_TAC[] THEN + DISCH_THEN(fun th -> + X_GEN_TAC `z:complex` THEN DISCH_TAC THEN MP_TAC th) THEN + ASM_REWRITE_TAC[divides; FUN_EQ_THM; POLY_MUL] THEN + DISCH_THEN(CHOOSE_THEN (MP_TAC o SPEC `z:complex`)) THEN + ASM_REWRITE_TAC[POLY_EXP; COMPLEX_MUL_LZERO; COMPLEX_POW_EQ_0]] THEN + ASM_REWRITE_TAC[poly_exp] THEN + SUBGOAL_THEN `?k. ~(k = Cx(&0)) /\ (poly p = poly [k])` MP_TAC THENL + [SUBST1_TAC(SYM(SPEC `p:complex list` POLY_NORMALIZE)) THEN + EXISTS_TAC `LAST(normalize p)` THEN + ASM_SIMP_TAC[NORMAL_NORMALIZE; GSYM POLY_NORMALIZE_ZERO] THEN + GEN_REWRITE_TAC LAND_CONV [GSYM POLY_NORMALIZE] THEN + UNDISCH_TAC `degree p = 0` THEN + FIRST_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV + [POLY_NORMALIZE_ZERO]) THEN + REWRITE_TAC[degree] THEN + SPEC_TAC(`normalize p`,`p:complex list`) THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[NOT_CONS_NIL] THEN + REWRITE_TAC[LENGTH; PRE; poly; LAST] THEN + SIMP_TAC[LENGTH_EQ_NIL; POLY_NORMALIZE]; ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `k:complex` STRIP_ASSUME_TAC) THEN + ASM_REWRITE_TAC[divides; poly; FUN_EQ_THM; POLY_MUL] THEN + ASM_REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + EXISTS_TAC `[inv(k)]` THEN + ASM_REWRITE_TAC[divides; poly; FUN_EQ_THM; POLY_MUL] THEN + ASM_REWRITE_TAC[COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + ASM_SIMP_TAC[COMPLEX_MUL_RINV]);; + +(* ------------------------------------------------------------------------- *) +(* Useful lemma I should have proved ages ago. *) +(* ------------------------------------------------------------------------- *) + +let CONSTANT_DEGREE = prove + (`!p. constant(poly p) <=> (degree p = 0)`, + GEN_TAC THEN REWRITE_TAC[constant] THEN EQ_TAC THENL + [DISCH_THEN(ASSUME_TAC o GSYM o SPEC `Cx(&0)`) THEN + SUBGOAL_THEN `degree [poly p (Cx(&0))] = 0` MP_TAC THENL + [REWRITE_TAC[degree; normalize] THEN + COND_CASES_TAC THEN REWRITE_TAC[LENGTH] THEN CONV_TAC NUM_REDUCE_CONV; + ALL_TAC] THEN + MATCH_MP_TAC(ARITH_RULE `(x = y) ==> (x = 0) ==> (y = 0)`) THEN + MATCH_MP_TAC DEGREE_WELLDEF THEN + REWRITE_TAC[FUN_EQ_THM; poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + FIRST_ASSUM(ACCEPT_TAC o GSYM); + ONCE_REWRITE_TAC[GSYM POLY_NORMALIZE] THEN REWRITE_TAC[degree] THEN + SPEC_TAC(`normalize p`,`l:complex list`) THEN + MATCH_MP_TAC list_INDUCT THEN REWRITE_TAC[poly] THEN + SIMP_TAC[LENGTH; PRE; LENGTH_EQ_NIL; poly; COMPLEX_MUL_RZERO]]);; + +(* ------------------------------------------------------------------------- *) +(* It would be nicer to prove this without using algebraic closure... *) +(* ------------------------------------------------------------------------- *) + +let DIVIDES_DEGREE_LEMMA = prove + (`!n p q. (degree(p) = n) + ==> n <= degree(p ** q) \/ (poly(p ** q) = poly [])`, + INDUCT_TAC THEN REWRITE_TAC[LE_0] THEN REPEAT STRIP_TAC THEN + MP_TAC(SPEC `p:complex list` FUNDAMENTAL_THEOREM_OF_ALGEBRA) THEN + ASM_REWRITE_TAC[CONSTANT_DEGREE; NOT_SUC] THEN + DISCH_THEN(X_CHOOSE_THEN `a:complex` MP_TAC) THEN + GEN_REWRITE_TAC LAND_CONV [POLY_LINEAR_DIVIDES] THEN + DISCH_THEN(DISJ_CASES_THEN2 SUBST_ALL_TAC MP_TAC) THENL + [REWRITE_TAC[POLY_MUL; poly; COMPLEX_MUL_LZERO; FUN_EQ_THM]; + ALL_TAC] THEN + DISCH_THEN(X_CHOOSE_THEN `r:complex list` SUBST_ALL_TAC) THEN + SUBGOAL_THEN `poly (([--a; Cx (&1)] ** r) ** q) = + poly ([--a; Cx (&1)] ** (r ** q))` + ASSUME_TAC THENL + [REWRITE_TAC[FUN_EQ_THM; POLY_MUL; COMPLEX_MUL_ASSOC]; ALL_TAC] THEN + FIRST_ASSUM(SUBST1_TAC o MATCH_MP DEGREE_WELLDEF) THEN + ASM_REWRITE_TAC[] THEN + MP_TAC(SPECL [`r ** q`; `--a`] LINEAR_MUL_DEGREE) THEN + ASM_CASES_TAC `poly (r ** q) = poly []` THENL + [REWRITE_TAC[FUN_EQ_THM] THEN + ONCE_REWRITE_TAC[POLY_MUL] THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO]; ALL_TAC] THEN + ASM_REWRITE_TAC[] THEN DISCH_TAC THEN ASM_REWRITE_TAC[] THEN + SUBGOAL_THEN `n <= degree(r ** q) \/ (poly(r ** q) = poly [])` MP_TAC THENL + [ALL_TAC; + REWRITE_TAC[ARITH_RULE `SUC n <= m + 1 <=> n <= m`] THEN + STRIP_TAC THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[FUN_EQ_THM] THEN + ONCE_REWRITE_TAC[POLY_MUL] THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO]] THEN + MP_TAC(SPECL [`r:complex list`; `--a`] LINEAR_MUL_DEGREE) THEN ANTS_TAC THENL + [UNDISCH_TAC `~(poly (r ** q) = poly [])` THEN + REWRITE_TAC[TAUT `~b ==> ~a <=> a ==> b`] THEN + SIMP_TAC[poly; FUN_EQ_THM; POLY_MUL; COMPLEX_ENTIRE]; ALL_TAC] THEN + DISCH_THEN SUBST_ALL_TAC THEN FIRST_ASSUM MATCH_MP_TAC THEN + UNDISCH_TAC `degree r + 1 = SUC n` THEN ARITH_TAC);; + +let DIVIDES_DEGREE = prove + (`!p q. p divides q ==> degree(p) <= degree(q) \/ (poly q = poly [])`, + REPEAT GEN_TAC THEN REWRITE_TAC[divides; LEFT_IMP_EXISTS_THM] THEN + X_GEN_TAC `r:complex list` THEN DISCH_TAC THEN + FIRST_ASSUM(SUBST1_TAC o MATCH_MP DEGREE_WELLDEF) THEN ASM_REWRITE_TAC[] THEN + ASM_MESON_TAC[DIVIDES_DEGREE_LEMMA]);; + +(* ------------------------------------------------------------------------- *) +(* Arithmetic operations on multivariate polynomials. *) +(* ------------------------------------------------------------------------- *) + +let MPOLY_BASE_CONV = + let pth_0 = prove + (`Cx(&0) = poly [] x`, + REWRITE_TAC[poly]) + and pth_1 = prove + (`c = poly [c] x`, + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID]) + and pth_var = prove + (`x = poly [Cx(&0); Cx(&1)] x`, + REWRITE_TAC[poly; COMPLEX_ADD_LID; COMPLEX_MUL_RZERO] THEN + REWRITE_TAC[COMPLEX_ADD_RID; COMPLEX_MUL_RID]) + and zero_tm = `Cx(&0)` + and c_tm = `c:complex` + and x_tm = `x:complex` in + let rec MPOLY_BASE_CONV avs tm = + if avs = [] then REFL tm + else if tm = zero_tm then INST [hd avs,x_tm] pth_0 + else if tm = hd avs then + let th1 = INST [tm,x_tm] pth_var in + let th2 = + (LAND_CONV + (COMB2_CONV (RAND_CONV (MPOLY_BASE_CONV (tl avs))) + (LAND_CONV (MPOLY_BASE_CONV (tl avs))))) + (rand(concl th1)) in + TRANS th1 th2 + else + let th1 = MPOLY_BASE_CONV (tl avs) tm in + let th2 = INST [hd avs,x_tm; rand(concl th1),c_tm] pth_1 in + TRANS th1 th2 in + MPOLY_BASE_CONV;; + +let MPOLY_NORM_CONV = + let pth_0 = prove + (`poly [Cx(&0)] x = poly [] x`, + REWRITE_TAC[poly; COMPLEX_ADD_RID; COMPLEX_MUL_RZERO]) + and pth_1 = prove + (`poly [poly [] y] x = poly [] x`, + REWRITE_TAC[poly; COMPLEX_ADD_RID; COMPLEX_MUL_RZERO]) in + let conv_fwd = REWR_CONV(CONJUNCT2 poly) + and conv_bck = REWR_CONV(GSYM(CONJUNCT2 poly)) + and conv_0 = GEN_REWRITE_CONV I [pth_0] + and conv_1 = GEN_REWRITE_CONV I [pth_1] in + let rec NORM0_CONV tm = + (conv_0 ORELSEC + (conv_fwd THENC RAND_CONV(RAND_CONV NORM0_CONV) THENC conv_bck THENC + TRY_CONV NORM0_CONV)) tm + and NORM1_CONV tm = + (conv_1 ORELSEC + (conv_fwd THENC RAND_CONV(RAND_CONV NORM1_CONV) THENC conv_bck THENC + TRY_CONV NORM1_CONV)) tm in + fun avs -> TRY_CONV(if avs = [] then NORM0_CONV else NORM1_CONV);; + +let MPOLY_ADD_CONV,MPOLY_TADD_CONV = + let add_conv0 = GEN_REWRITE_CONV I (butlast (CONJUNCTS POLY_ADD_CLAUSES)) + and add_conv1 = GEN_REWRITE_CONV I [last (CONJUNCTS POLY_ADD_CLAUSES)] + and add_conv = REWR_CONV(GSYM POLY_ADD) in + let rec MPOLY_ADD_CONV avs tm = + if avs = [] then COMPLEX_RAT_ADD_CONV tm else + (add_conv THENC LAND_CONV(MPOLY_TADD_CONV avs) THENC + MPOLY_NORM_CONV (tl avs)) tm + and MPOLY_TADD_CONV avs tm = + (add_conv0 ORELSEC + (add_conv1 THENC + LAND_CONV (MPOLY_ADD_CONV (tl avs)) THENC + RAND_CONV (MPOLY_TADD_CONV avs))) tm in + MPOLY_ADD_CONV,MPOLY_TADD_CONV;; + +let MPOLY_CMUL_CONV,MPOLY_TCMUL_CONV,MPOLY_MUL_CONV,MPOLY_TMUL_CONV = + let cmul_conv0 = GEN_REWRITE_CONV I [CONJUNCT1 poly_cmul] + and cmul_conv1 = GEN_REWRITE_CONV I [CONJUNCT2 poly_cmul] + and cmul_conv = REWR_CONV(GSYM POLY_CMUL) + and mul_conv0 = GEN_REWRITE_CONV I [CONJUNCT1 POLY_MUL_CLAUSES] + and mul_conv1 = GEN_REWRITE_CONV I [CONJUNCT1(CONJUNCT2 POLY_MUL_CLAUSES)] + and mul_conv2 = GEN_REWRITE_CONV I [CONJUNCT2(CONJUNCT2 POLY_MUL_CLAUSES)] + and mul_conv = REWR_CONV(GSYM POLY_MUL) in + let rec MPOLY_CMUL_CONV avs tm = + (cmul_conv THENC LAND_CONV(MPOLY_TCMUL_CONV avs)) tm + and MPOLY_TCMUL_CONV avs tm = + (cmul_conv0 ORELSEC + (cmul_conv1 THENC + LAND_CONV (MPOLY_MUL_CONV (tl avs)) THENC + RAND_CONV (MPOLY_TCMUL_CONV avs))) tm + and MPOLY_MUL_CONV avs tm = + if avs = [] then COMPLEX_RAT_MUL_CONV tm else + (mul_conv THENC LAND_CONV(MPOLY_TMUL_CONV avs)) tm + and MPOLY_TMUL_CONV avs tm = + (mul_conv0 ORELSEC + (mul_conv1 THENC MPOLY_TCMUL_CONV avs) ORELSEC + (mul_conv2 THENC + COMB2_CONV (RAND_CONV(MPOLY_TCMUL_CONV avs)) + (COMB2_CONV (RAND_CONV(MPOLY_BASE_CONV (tl avs))) + (MPOLY_TMUL_CONV avs)) THENC + MPOLY_TADD_CONV avs)) tm in + MPOLY_CMUL_CONV,MPOLY_TCMUL_CONV,MPOLY_MUL_CONV,MPOLY_TMUL_CONV;; + +let MPOLY_SUB_CONV = + let pth = prove + (`(poly p x - poly q x) = (poly p x + Cx(--(&1)) * poly q x)`, + SIMPLE_COMPLEX_ARITH_TAC) in + let APPLY_PTH_CONV = REWR_CONV pth in + fun avs -> + APPLY_PTH_CONV THENC + RAND_CONV(LAND_CONV (MPOLY_BASE_CONV (tl avs)) THENC + MPOLY_CMUL_CONV avs) THENC + MPOLY_ADD_CONV avs;; + +let MPOLY_POW_CONV = + let cnv_0 = GEN_REWRITE_CONV I [CONJUNCT1 complex_pow] + and cnv_1 = GEN_REWRITE_CONV I [CONJUNCT2 complex_pow] in + let rec MPOLY_POW_CONV avs tm = + try (cnv_0 THENC MPOLY_BASE_CONV avs) tm with Failure _ -> + (RAND_CONV num_CONV THENC + cnv_1 THENC (RAND_CONV (MPOLY_POW_CONV avs)) THENC + MPOLY_MUL_CONV avs) tm in + MPOLY_POW_CONV;; + +(* ------------------------------------------------------------------------- *) +(* Recursive conversion to polynomial form. *) +(* ------------------------------------------------------------------------- *) + +let POLYNATE_CONV = + let ELIM_SUB_CONV = REWR_CONV + (SIMPLE_COMPLEX_ARITH `x - y = x + Cx(--(&1)) * y`) + and ELIM_NEG_CONV = REWR_CONV + (SIMPLE_COMPLEX_ARITH `--x = Cx(--(&1)) * x`) + and ELIM_POW_0_CONV = GEN_REWRITE_CONV I [CONJUNCT1 complex_pow] + and ELIM_POW_1_CONV = + RAND_CONV num_CONV THENC GEN_REWRITE_CONV I [CONJUNCT2 complex_pow] in + let rec ELIM_POW_CONV tm = + (ELIM_POW_0_CONV ORELSEC (ELIM_POW_1_CONV THENC RAND_CONV ELIM_POW_CONV)) + tm in + let polynet = itlist (uncurry net_of_conv) + [`x pow n`,(fun cnv avs -> LAND_CONV (cnv avs) THENC MPOLY_POW_CONV avs); + `x * y`,(fun cnv avs -> BINOP_CONV (cnv avs) THENC MPOLY_MUL_CONV avs); + `x + y`,(fun cnv avs -> BINOP_CONV (cnv avs) THENC MPOLY_ADD_CONV avs); + `x - y`,(fun cnv avs -> BINOP_CONV (cnv avs) THENC MPOLY_SUB_CONV avs); + `--x`,(fun cnv avs -> ELIM_NEG_CONV THENC (cnv avs))] + empty_net in + let rec POLYNATE_CONV avs tm = + try snd(hd(lookup tm polynet)) POLYNATE_CONV avs tm + with Failure _ -> MPOLY_BASE_CONV avs tm in + POLYNATE_CONV;; + +(* ------------------------------------------------------------------------- *) +(* Cancellation conversion. *) +(* ------------------------------------------------------------------------- *) + +let POLY_PAD_RULE = + let pth = prove + (`(poly p x = Cx(&0)) ==> (poly (CONS (Cx(&0)) p) x = Cx(&0))`, + SIMP_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_LID]) in + let MATCH_pth = MATCH_MP pth in + fun avs th -> + let th1 = MATCH_pth th in + CONV_RULE(funpow 3 LAND_CONV (MPOLY_BASE_CONV (tl avs))) th1;; + +let POLY_CANCEL_EQ_CONV = + let pth_1 = prove + (`(p = Cx(&0)) /\ ~(a = Cx(&0)) + ==> !q b. (q = Cx(&0)) <=> (a * q - b * p = Cx(&0))`, + SIMP_TAC[COMPLEX_MUL_RZERO; COMPLEX_SUB_RZERO; COMPLEX_ENTIRE]) in + let MATCH_CANCEL_THM = MATCH_MP pth_1 in + let rec POLY_CANCEL_EQ_CONV avs n ath eth tm = + let m = length(dest_list(lhand(lhand tm))) in + if m < n then REFL tm else + let th1 = funpow (m - n) (POLY_PAD_RULE avs) eth in + let th2 = MATCH_CANCEL_THM (CONJ th1 ath) in + let th3 = SPECL [lhs tm; last(dest_list(lhand(lhs tm)))] th2 in + let th4 = CONV_RULE(RAND_CONV(LAND_CONV + (BINOP_CONV(MPOLY_CMUL_CONV avs)))) th3 in + let th5 = CONV_RULE(RAND_CONV(LAND_CONV (MPOLY_SUB_CONV avs))) th4 in + TRANS th5 (POLY_CANCEL_EQ_CONV avs n ath eth (rand(concl th5))) in + POLY_CANCEL_EQ_CONV;; + +let RESOLVE_EQ_RAW = + let pth = prove + (`(poly [] x = Cx(&0)) /\ + (poly [c] x = c)`, + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID]) in + let REWRITE_pth = GEN_REWRITE_CONV LAND_CONV [pth] in + let rec RESOLVE_EQ asm tm = + try EQT_INTRO(find (fun th -> concl th = tm) asm) with Failure _ -> + let tm' = mk_neg tm in + try EQF_INTRO(find (fun th -> concl th = tm') asm) with Failure _ -> + try let th1 = REWRITE_pth tm in + TRANS th1 (RESOLVE_EQ asm (rand(concl th1))) + with Failure _ -> COMPLEX_RAT_EQ_CONV tm in + RESOLVE_EQ;; + +let RESOLVE_EQ asm tm = + let th = RESOLVE_EQ_RAW asm tm in + try EQF_ELIM th with Failure _ -> EQT_ELIM th;; + +let RESOLVE_EQ_THEN = + let MATCH_pth = MATCH_MP + (TAUT `(p ==> (q <=> q1)) /\ (~p ==> (q <=> q2)) + ==> (q <=> (p /\ q1 \/ ~p /\ q2))`) in + fun asm tm yfn nfn -> + try let th = RESOLVE_EQ asm tm in + if is_neg(concl th) then nfn (th::asm) th else yfn (th::asm) th + with Failure _ -> + let tm' = mk_neg tm in + let yth = DISCH tm (yfn (ASSUME tm :: asm) (ASSUME tm)) + and nth = DISCH tm' (nfn (ASSUME tm' :: asm) (ASSUME tm')) in + MATCH_pth (CONJ yth nth);; + +let POLY_CANCEL_ENE_CONV avs n ath eth tm = + if is_neg tm then RAND_CONV(POLY_CANCEL_EQ_CONV avs n ath eth) tm + else POLY_CANCEL_EQ_CONV avs n ath eth tm;; + +let RESOLVE_NE = + let NEGATE_NEGATE_RULE = GEN_REWRITE_RULE I [TAUT `p <=> (~p <=> F)`] in + fun asm tm -> + try let th = RESOLVE_EQ asm (rand tm) in + if is_neg(concl th) then EQT_INTRO th + else NEGATE_NEGATE_RULE th + with Failure _ -> REFL tm;; + +(* ------------------------------------------------------------------------- *) +(* Conversion for division of polynomials. *) +(* ------------------------------------------------------------------------- *) + +let LAST_CONV = GEN_REWRITE_CONV REPEATC [LAST_CLAUSES];; + +let LENGTH_CONV = + let cnv_0 = GEN_REWRITE_CONV I [CONJUNCT1 LENGTH] + and cnv_1 = GEN_REWRITE_CONV I [CONJUNCT2 LENGTH] in + let rec LENGTH_CONV tm = + try cnv_0 tm with Failure _ -> + (cnv_1 THENC RAND_CONV LENGTH_CONV THENC NUM_SUC_CONV) tm in + LENGTH_CONV;; + +let EXPAND_EX_BETA_CONV = + let pth = prove(`EX P [c] = P c`,REWRITE_TAC[EX]) in + let cnv_0 = GEN_REWRITE_CONV I [CONJUNCT1 EX] + and cnv_1 = GEN_REWRITE_CONV I [pth] + and cnv_2 = GEN_REWRITE_CONV I [CONJUNCT2 EX] in + let rec EXPAND_EX_BETA_CONV tm = + try (cnv_1 THENC BETA_CONV) tm with Failure _ -> try + (cnv_2 THENC COMB2_CONV (RAND_CONV BETA_CONV) + EXPAND_EX_BETA_CONV) tm + with Failure _ -> cnv_0 tm in + EXPAND_EX_BETA_CONV;; + +let POLY_DIVIDES_PAD_RULE = + let pth = prove + (`p divides q ==> p divides (CONS (Cx(&0)) q)`, + REWRITE_TAC[divides; FUN_EQ_THM; POLY_MUL; poly; COMPLEX_ADD_LID] THEN + DISCH_THEN(X_CHOOSE_THEN `r:complex list` ASSUME_TAC) THEN + EXISTS_TAC `[Cx(&0); Cx(&1)] ** r` THEN + ASM_REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_LID; COMPLEX_ADD_RID; + COMPLEX_MUL_RID; POLY_MUL] THEN + REWRITE_TAC[COMPLEX_MUL_AC]) in + let APPLY_pth = MATCH_MP pth in + fun avs n tm -> + funpow n + (CONV_RULE(RAND_CONV(LAND_CONV(MPOLY_BASE_CONV (tl avs)))) o APPLY_pth) + (SPEC tm POLY_DIVIDES_REFL);; + +let POLY_DIVIDES_PAD_CONST_RULE = + let pth = prove + (`p divides q ==> !a. p divides (a ## q)`, + REWRITE_TAC[FUN_EQ_THM; divides; POLY_CMUL; POLY_MUL] THEN + DISCH_THEN(X_CHOOSE_THEN `r:complex list` ASSUME_TAC) THEN + X_GEN_TAC `a:complex` THEN EXISTS_TAC `[a] ** r` THEN + ASM_REWRITE_TAC[POLY_MUL; poly] THEN SIMPLE_COMPLEX_ARITH_TAC) in + let APPLY_pth = MATCH_MP pth in + fun avs n a tm -> + let th1 = POLY_DIVIDES_PAD_RULE avs n tm in + let th2 = SPEC a (APPLY_pth th1) in + CONV_RULE(RAND_CONV(MPOLY_TCMUL_CONV avs)) th2;; + +let EXPAND_EX_BETA_RESOLVE_CONV asm tm = + let th1 = EXPAND_EX_BETA_CONV tm in + let djs = disjuncts(rand(concl th1)) in + let th2 = end_itlist MK_DISJ (map (RESOLVE_NE asm) djs) in + TRANS th1 th2;; + +let POLY_DIVIDES_CONV = + let pth_0 = prove + (`LENGTH q < LENGTH p + ==> ~(LAST p = Cx(&0)) + ==> (p divides q <=> ~(EX (\c. ~(c = Cx(&0))) q))`, + REPEAT STRIP_TAC THEN REWRITE_TAC[NOT_EX; GSYM POLY_ZERO] THEN EQ_TAC THENL + [ALL_TAC; + SIMP_TAC[divides; POLY_MUL; FUN_EQ_THM] THEN + DISCH_TAC THEN EXISTS_TAC `[]:complex list` THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO]] THEN + DISCH_THEN(MP_TAC o MATCH_MP DIVIDES_DEGREE) THEN + MATCH_MP_TAC(TAUT `(~b ==> ~a) ==> (a \/ b ==> b)`) THEN + DISCH_TAC THEN REWRITE_TAC[NOT_LE] THEN ASM_SIMP_TAC[NORMAL_DEGREE] THEN + REWRITE_TAC[degree] THEN + FIRST_ASSUM(MATCH_MP_TAC o MATCH_MP (ARITH_RULE + `lq < lp ==> ~(lq = 0) /\ dq <= lq - 1 ==> dq < lp - 1`)) THEN + CONJ_TAC THENL [ASM_MESON_TAC[LENGTH_EQ_NIL]; ALL_TAC] THEN + MATCH_MP_TAC(ARITH_RULE `m <= n ==> PRE m <= n - 1`) THEN + REWRITE_TAC[LENGTH_NORMALIZE_LE]) in + let APPLY_pth0 = PART_MATCH (lhand o rand o rand) pth_0 in + let pth_1 = prove + (`~(a = Cx(&0)) + ==> p divides p' + ==> (!x. a * poly q x - poly p' x = poly r x) + ==> (p divides q <=> p divides r)`, + DISCH_TAC THEN REWRITE_TAC[divides; LEFT_IMP_EXISTS_THM] THEN + X_GEN_TAC `t:complex list` THEN DISCH_THEN SUBST1_TAC THEN + REWRITE_TAC[FUN_EQ_THM; POLY_MUL] THEN + DISCH_THEN(fun th -> REWRITE_TAC[GSYM th]) THEN EQ_TAC THEN + DISCH_THEN(X_CHOOSE_THEN `s:complex list` MP_TAC) THENL + [DISCH_THEN(fun th -> REWRITE_TAC[th]) THEN + EXISTS_TAC `a ## s ++ --(Cx(&1)) ## t` THEN + REWRITE_TAC[POLY_MUL; POLY_ADD; POLY_CMUL] THEN + REWRITE_TAC[poly] THEN SIMPLE_COMPLEX_ARITH_TAC; + REWRITE_TAC[POLY_MUL] THEN DISCH_TAC THEN + EXISTS_TAC `[inv(a)] ** (t ++ s)` THEN + X_GEN_TAC `z:complex` THEN + ONCE_REWRITE_TAC[COMPLEX_MUL_SYM] THEN + REWRITE_TAC[POLY_MUL; POLY_ADD; GSYM COMPLEX_MUL_ASSOC] THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID] THEN + SUBGOAL_THEN `a * poly q z = (poly t z + poly s z) * poly p z` + MP_TAC THENL + [FIRST_ASSUM(MP_TAC o SPEC `z:complex`) THEN SIMPLE_COMPLEX_ARITH_TAC; + ALL_TAC] THEN + DISCH_THEN(MP_TAC o AP_TERM `( * ) (inv a)`) THEN + ASM_SIMP_TAC[COMPLEX_MUL_ASSOC; COMPLEX_MUL_LINV; COMPLEX_MUL_LID]]) in + let MATCH_pth1 = MATCH_MP pth_1 in + let rec DIVIDE_STEP_CONV avs sfn n tm = + let m = length(dest_list(rand tm)) in + if m < n then REFL tm else + let th1 = POLY_DIVIDES_PAD_CONST_RULE avs (m - n) + (last(dest_list(rand tm))) (lhand tm) in + let th2 = MATCH_MP (sfn tm) th1 in + let av,bod = dest_forall(lhand(concl th2)) in + let tm1 = vsubst [hd avs,av] (lhand bod) in + let th3 = (LAND_CONV (MPOLY_CMUL_CONV avs) THENC MPOLY_SUB_CONV avs) tm1 in + let th4 = MATCH_MP th2 (GEN (hd avs) th3) in + TRANS th4 (DIVIDE_STEP_CONV avs sfn n (rand(concl th4))) in + let zero_tm = `Cx(&0)` in + fun asm avs tm -> + let ath = RESOLVE_EQ asm (mk_eq(last(dest_list(lhand tm)),zero_tm)) in + let sfn = PART_MATCH (lhand o rand o rand) (MATCH_pth1 ath) + and n = length(dest_list(lhand tm)) in + let th1 = DIVIDE_STEP_CONV avs sfn n tm in + let th2 = APPLY_pth0 (rand(concl th1)) in + let th3 = (BINOP_CONV LENGTH_CONV THENC NUM_LT_CONV) (lhand(concl th2)) in + let th4 = MP th2 (EQT_ELIM th3) in + let th5 = CONV_RULE(LAND_CONV(RAND_CONV(LAND_CONV LAST_CONV))) th4 in + let th6 = TRANS th1 (MP th5 ath) in + CONV_RULE(RAND_CONV(RAND_CONV(EXPAND_EX_BETA_RESOLVE_CONV asm))) th6;; + +(* ------------------------------------------------------------------------- *) +(* Apply basic Nullstellensatz principle. *) +(* ------------------------------------------------------------------------- *) + +let BASIC_QUELIM_CONV = + let pth_1 = prove + (`((?x. (poly p x = Cx(&0)) /\ ~(poly [] x = Cx(&0))) <=> F) /\ + ((?x. ~(poly [] x = Cx(&0))) <=> F) /\ + ((?x. ~(poly [c] x = Cx(&0))) <=> ~(c = Cx(&0))) /\ + ((?x. (poly [] x = Cx(&0))) <=> T) /\ + ((?x. (poly [c] x = Cx(&0))) <=> (c = Cx(&0)))`, + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID]) in + let APPLY_pth1 = GEN_REWRITE_CONV I [pth_1] in + let pth_2 = prove + (`~(LAST(CONS a (CONS b p)) = Cx(&0)) + ==> ((?x. poly (CONS a (CONS b p)) x = Cx(&0)) <=> T)`, + REPEAT STRIP_TAC THEN + MP_TAC(SPEC `CONS (a:complex) (CONS b p)` + FUNDAMENTAL_THEOREM_OF_ALGEBRA_ALT) THEN + REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN + REWRITE_TAC[NOT_EXISTS_THM; CONS_11] THEN + REPEAT STRIP_TAC THEN + SUBGOAL_THEN `~(ALL (\c. c = Cx(&0)) (CONS b p))` + (fun th -> MP_TAC th THEN ASM_REWRITE_TAC[]) THEN + UNDISCH_TAC `~(LAST (CONS a (CONS b p)) = Cx (&0))` THEN + ONCE_REWRITE_TAC[LAST] THEN REWRITE_TAC[NOT_CONS_NIL] THEN + REWRITE_TAC[TAUT `~a ==> ~b <=> b ==> a`] THEN + SPEC_TAC(`p:complex list`,`p:complex list`) THEN + LIST_INDUCT_TAC THEN ONCE_REWRITE_TAC[LAST] THEN + REWRITE_TAC[ALL; NOT_CONS_NIL] THEN + STRIP_TAC THEN FIRST_ASSUM(UNDISCH_TAC o check is_imp o concl) THEN + REWRITE_TAC[LAST] THEN COND_CASES_TAC THEN ASM_REWRITE_TAC[ALL]) in + let APPLY_pth2 = PART_MATCH (lhand o rand) pth_2 in + let pth_2b = prove + (`(?x. ~(poly p x = Cx(&0))) <=> EX (\c. ~(c = Cx(&0))) p`, + REWRITE_TAC[GSYM NOT_FORALL_THM] THEN + ONCE_REWRITE_TAC[TAUT `(~a <=> b) <=> (a <=> ~b)`] THEN + REWRITE_TAC[NOT_EX; GSYM POLY_ZERO; poly; FUN_EQ_THM]) in + let APPLY_pth2b = GEN_REWRITE_CONV I [pth_2b] in + let pth_3 = prove + (`~(LAST(CONS a p) = Cx(&0)) + ==> ((?x. (poly (CONS a p) x = Cx(&0)) /\ ~(poly q x = Cx(&0))) <=> + ~((CONS a p) divides (q exp (LENGTH p))))`, + REPEAT STRIP_TAC THEN + MP_TAC(SPECL [`CONS (a:complex) p`; `q:complex list`] + NULLSTELLENSATZ_UNIVARIATE) THEN + ASM_SIMP_TAC[degree; NORMALIZE_EQ; LENGTH; PRE] THEN + SUBGOAL_THEN `~(poly (CONS a p) = poly [])` + (fun th -> REWRITE_TAC[th] THEN MESON_TAC[]) THEN + REWRITE_TAC[POLY_ZERO] THEN POP_ASSUM MP_TAC THEN + SPEC_TAC(`p:complex list`,`p:complex list`) THEN + REWRITE_TAC[LAST] THEN + LIST_INDUCT_TAC THEN REWRITE_TAC[LAST; ALL; NOT_CONS_NIL] THEN + POP_ASSUM MP_TAC THEN COND_CASES_TAC THEN ASM_SIMP_TAC[ALL] THEN + CONV_TAC TAUT) in + let APPLY_pth3 = PART_MATCH (lhand o rand) pth_3 in + let POLY_EXP_DIVIDES_CONV = + let pth_4 = prove + (`(!x. (poly (q exp n) x = poly r x)) + ==> (p divides (q exp n) <=> p divides r)`, + SIMP_TAC[divides; POLY_EXP; FUN_EQ_THM]) in + let APPLY_pth4 = MATCH_MP pth_4 + and poly_tm = `poly` + and REWR_POLY_EXP_CONV = REWR_CONV POLY_EXP in + let POLY_EXP_DIVIDES_CONV avs tm = + let tm1 = mk_comb(mk_comb(poly_tm,rand tm),hd avs) in + let th1 = REWR_POLY_EXP_CONV tm1 in + let th2 = TRANS th1 (MPOLY_POW_CONV avs (rand(concl th1))) in + PART_MATCH lhand (APPLY_pth4 (GEN (hd avs) th2)) tm in + POLY_EXP_DIVIDES_CONV in + fun asm avs tm -> + try APPLY_pth1 tm with Failure _ -> + try let th1 = APPLY_pth2 tm in + let th2 = CONV_RULE(LAND_CONV(RAND_CONV(LAND_CONV LAST_CONV))) th1 in + let th3 = try MATCH_MP th2 (RESOLVE_EQ asm (rand(lhand(concl th2)))) + with Failure _ -> failwith "Sanity failure (2a)" in + th3 + with Failure _ -> try + let th1 = APPLY_pth2b tm in + TRANS th1 (EXPAND_EX_BETA_RESOLVE_CONV asm (rand(concl th1))) + with Failure _ -> + let th1 = APPLY_pth3 tm in + let th2 = CONV_RULE(LAND_CONV(RAND_CONV(LAND_CONV LAST_CONV))) th1 in + let th3 = try MATCH_MP th2 (RESOLVE_EQ asm (rand(lhand(concl th2)))) + with Failure _ -> failwith "Sanity failure (2b)" in + let th4 = CONV_RULE (funpow 4 RAND_CONV LENGTH_CONV) th3 in + let th5 = + CONV_RULE(RAND_CONV(RAND_CONV(POLY_EXP_DIVIDES_CONV avs))) th4 in + CONV_RULE(RAND_CONV(RAND_CONV(POLY_DIVIDES_CONV asm avs))) th5;; + +(* ------------------------------------------------------------------------- *) +(* Put into canonical form by multiplying inequalities. *) +(* ------------------------------------------------------------------------- *) + +let POLY_NE_MULT_CONV = + let pth = prove + (`~(poly p x = Cx(&0)) /\ ~(poly q x = Cx(&0)) <=> + ~(poly p x * poly q x = Cx(&0))`, + REWRITE_TAC[COMPLEX_ENTIRE; DE_MORGAN_THM]) in + let APPLY_pth = REWR_CONV pth in + let rec POLY_NE_MULT_CONV avs tm = + if not(is_conj tm) then REFL tm else + let l,r = dest_conj tm in + let th1 = MK_COMB(AP_TERM (rator(rator tm)) (POLY_NE_MULT_CONV avs l), + POLY_NE_MULT_CONV avs r) in + let th2 = TRANS th1 (APPLY_pth (rand(concl th1))) in + CONV_RULE(RAND_CONV(RAND_CONV(LAND_CONV(MPOLY_MUL_CONV avs)))) th2 in + POLY_NE_MULT_CONV;; + +let CORE_QUELIM_CONV = + let CONJ_AC_RULE = AC CONJ_ACI in + let CORE_QUELIM_CONV asm avs tm = + let ev,bod = dest_exists tm in + let cjs = conjuncts bod in + let eqs,neqs = partition is_eq cjs in + if eqs = [] then + let th1 = MK_EXISTS ev (POLY_NE_MULT_CONV avs bod) in + TRANS th1 (BASIC_QUELIM_CONV asm avs (rand(concl th1))) + else if length eqs > 1 then failwith "CORE_QUELIM_CONV: Sanity failure" + else if neqs = [] then BASIC_QUELIM_CONV asm avs tm else + let tm1 = mk_conj(hd eqs,list_mk_conj neqs) in + let th1 = CONJ_AC_RULE(mk_eq(bod,tm1)) in + let th2 = CONV_RULE(funpow 2 RAND_CONV(POLY_NE_MULT_CONV avs)) th1 in + let th3 = MK_EXISTS ev th2 in + TRANS th3 (BASIC_QUELIM_CONV asm avs (rand(concl th3))) in + CORE_QUELIM_CONV;; + +(* ------------------------------------------------------------------------- *) +(* Main elimination coversion (for a single quantifier). *) +(* ------------------------------------------------------------------------- *) + +let RESOLVE_EQ_NE = + let DNE_RULE = GEN_REWRITE_RULE I + [TAUT `((p <=> T) <=> (~p <=> F)) /\ ((p <=> F) <=> (~p <=> T))`] in + fun asm tm -> + if is_neg tm then DNE_RULE(RESOLVE_EQ_RAW asm (rand tm)) + else RESOLVE_EQ_RAW asm tm;; + +let COMPLEX_QUELIM_CONV = + let pth_0 = prove + (`((poly [] x = Cx(&0)) <=> T) /\ + ((poly [] x = Cx(&0)) /\ p <=> p)`, + REWRITE_TAC[poly]) + and pth_1 = prove + (`(~(poly [] x = Cx(&0)) <=> F) /\ + (~(poly [] x = Cx(&0)) /\ p <=> F)`, + REWRITE_TAC[poly]) + and pth_2 = prove + (`(p ==> (q <=> r)) ==> (p /\ q <=> p /\ r)`, + CONV_TAC TAUT) + and zero_tm = `Cx(&0)` + and true_tm = `T` in + let ELIM_ZERO_RULE = GEN_REWRITE_RULE RAND_CONV [pth_0] + and ELIM_NONZERO_RULE = GEN_REWRITE_RULE RAND_CONV [pth_1] + and INCORP_ASSUM_THM = MATCH_MP pth_2 + and CONJ_AC_RULE = AC CONJ_ACI in + let POLY_CONST_CONV = + let pth = prove + (`((poly [c] x = y) <=> (c = y)) /\ + (~(poly [c] x = y) <=> ~(c = y))`, + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_RID]) in + TRY_CONV(GEN_REWRITE_CONV I [pth]) in + let EXISTS_TRIV_CONV = REWR_CONV EXISTS_SIMP + and EXISTS_PUSH_CONV = REWR_CONV RIGHT_EXISTS_AND_THM + and AND_SIMP_CONV = GEN_REWRITE_CONV DEPTH_CONV + [TAUT `(p /\ F <=> F) /\ (p /\ T <=> p) /\ + (F /\ p <=> F) /\ (T /\ p <=> p)`] + and RESOLVE_OR_CONST_CONV asm tm = + try RESOLVE_EQ_NE asm tm with Failure _ -> POLY_CONST_CONV tm + and false_tm = `F` in + let rec COMPLEX_QUELIM_CONV asm avs tm = + let ev,bod = dest_exists tm in + let cjs = conjuncts bod in + let cjs_set = setify cjs in + if length cjs_set < length cjs then + let th1 = CONJ_AC_RULE(mk_eq(bod,list_mk_conj cjs_set)) in + let th2 = MK_EXISTS ev th1 in + TRANS th2 (COMPLEX_QUELIM_CONV asm avs (rand(concl th2))) + else + let eqs,neqs = partition is_eq cjs in + let lens = map (length o dest_list o lhand o lhs) eqs + and nens = map (length o dest_list o lhand o lhs o rand) neqs in + try let zeq = el (index 0 lens) eqs in + if cjs = [zeq] then BASIC_QUELIM_CONV asm avs tm else + let cjs' = zeq::(subtract cjs [zeq]) in + let th1 = ELIM_ZERO_RULE(CONJ_AC_RULE(mk_eq(bod,list_mk_conj cjs'))) in + let th2 = MK_EXISTS ev th1 in + TRANS th2 (COMPLEX_QUELIM_CONV asm avs (rand(concl th2))) + with Failure _ -> try + let zne = el (index 0 nens) neqs in + if cjs = [zne] then BASIC_QUELIM_CONV asm avs tm else + let cjs' = zne::(subtract cjs [zne]) in + let th1 = ELIM_NONZERO_RULE + (CONJ_AC_RULE(mk_eq(bod,list_mk_conj cjs'))) in + CONV_RULE (RAND_CONV EXISTS_TRIV_CONV) (MK_EXISTS ev th1) + with Failure _ -> try + let ones = map snd (filter (fun (n,_) -> n = 1) + (zip lens eqs @ zip nens neqs)) in + if ones = [] then failwith "" else + let cjs' = subtract cjs ones in + if cjs' = [] then + let th1 = MK_EXISTS ev (SUBS_CONV(map POLY_CONST_CONV cjs) bod) in + TRANS th1 (EXISTS_TRIV_CONV (rand(concl th1))) + else + let tha = SUBS_CONV (map (RESOLVE_OR_CONST_CONV asm) ones) + (list_mk_conj ones) in + let thb = CONV_RULE (RAND_CONV AND_SIMP_CONV) tha in + if rand(concl thb) = false_tm then + let thc = MK_CONJ thb (REFL(list_mk_conj cjs')) in + let thd = CONV_RULE(RAND_CONV AND_SIMP_CONV) thc in + let the = CONJ_AC_RULE(mk_eq(bod,lhand(concl thd))) in + let thf = MK_EXISTS ev (TRANS the thd) in + CONV_RULE(RAND_CONV EXISTS_TRIV_CONV) thf + else + let thc = MK_CONJ thb (REFL(list_mk_conj cjs')) in + let thd = CONJ_AC_RULE(mk_eq(bod,lhand(concl thc))) in + let the = MK_EXISTS ev (TRANS thd thc) in + let th4 = TRANS the(EXISTS_PUSH_CONV(rand(concl the))) in + let tm4 = rand(concl th4) in + let th5 = COMPLEX_QUELIM_CONV asm avs (rand tm4) in + TRANS th4 (AP_TERM (rator tm4) th5) + with Failure _ -> + if eqs = [] or + (length eqs = 1 & + (let ceq = mk_eq(last(dest_list(lhand(lhs(hd eqs)))),zero_tm) in + try concl(RESOLVE_EQ asm ceq) = mk_neg ceq with Failure _ -> false) & + (let h = hd lens in forall (fun n -> n < h) nens)) + then + CORE_QUELIM_CONV asm avs tm + else + let n = end_itlist min lens in + let eq = el (index n lens) eqs in + let pol = lhand(lhand eq) in + let atm = last(dest_list pol) in + let zeq = mk_eq(atm,zero_tm) in + RESOLVE_EQ_THEN asm zeq + (fun asm' yth -> + let th0 = TRANS yth (MPOLY_BASE_CONV (tl avs) zero_tm) in + let th1 = + GEN_REWRITE_CONV + (LAND_CONV o LAND_CONV o funpow (n - 1) RAND_CONV o LAND_CONV) + [th0] eq in + let th2 = LAND_CONV(MPOLY_NORM_CONV avs) (rand(concl th1)) in + let th3 = MK_EXISTS ev (SUBS_CONV[TRANS th1 th2] bod) in + TRANS th3 (COMPLEX_QUELIM_CONV asm' avs (rand(concl th3)))) + (fun asm' nth -> + let oth = subtract cjs [eq] in + if oth = [] then COMPLEX_QUELIM_CONV asm' avs tm else + let eth = ASSUME eq in + let ths = map (POLY_CANCEL_ENE_CONV avs n nth eth) oth in + let th1 = DISCH eq (end_itlist MK_CONJ ths) in + let th2 = INCORP_ASSUM_THM th1 in + let th3 = TRANS (CONJ_AC_RULE(mk_eq(bod,lhand(concl th2)))) th2 in + let th4 = MK_EXISTS ev th3 in + TRANS th4 (COMPLEX_QUELIM_CONV asm' avs (rand(concl th4)))) in + fun asm avs -> time(COMPLEX_QUELIM_CONV asm avs);; + +(* ------------------------------------------------------------------------- *) +(* NNF conversion doing "conditionals" ~(p /\ q \/ ~p /\ r) intelligently. *) +(* ------------------------------------------------------------------------- *) + +let NNF_COND_CONV = + let NOT_EXISTS_UNIQUE_THM = prove + (`~(?!x. P x) <=> (!x. ~P x) \/ ?x x'. P x /\ P x' /\ ~(x = x')`, + REWRITE_TAC[EXISTS_UNIQUE_THM; DE_MORGAN_THM; NOT_EXISTS_THM] THEN + REWRITE_TAC[NOT_FORALL_THM; NOT_IMP; CONJ_ASSOC]) in + let tauts = + [TAUT `~(~p) <=> p`; + TAUT `~(p /\ q) <=> ~p \/ ~q`; + TAUT `~(p \/ q) <=> ~p /\ ~q`; + TAUT `~(p ==> q) <=> p /\ ~q`; + TAUT `p ==> q <=> ~p \/ q`; + NOT_FORALL_THM; + NOT_EXISTS_THM; + EXISTS_UNIQUE_THM; + NOT_EXISTS_UNIQUE_THM; + TAUT `~(p <=> q) <=> (p /\ ~q) \/ (~p /\ q)`; + TAUT `(p <=> q) <=> (p /\ q) \/ (~p /\ ~q)`; + TAUT `~(p /\ q \/ ~p /\ r) <=> p /\ ~q \/ ~p /\ ~r`] in + GEN_REWRITE_CONV TOP_SWEEP_CONV tauts;; + +(* ------------------------------------------------------------------------- *) +(* Overall procedure for multiple quantifiers in any first order formula. *) +(* ------------------------------------------------------------------------- *) + +let FULL_COMPLEX_QUELIM_CONV = + let ELIM_FORALL_CONV = + let pth = prove(`(!x. P x) <=> ~(?x. ~(P x))`,MESON_TAC[]) in + REWR_CONV pth in + let ELIM_EQ_CONV = + let pth = SIMPLE_COMPLEX_ARITH `(x = y) <=> (x - y = Cx(&0))` + and zero_tm = `Cx(&0)` in + let REWR_pth = REWR_CONV pth in + fun avs tm -> + if rand tm = zero_tm then LAND_CONV(POLYNATE_CONV avs) tm + else (REWR_pth THENC LAND_CONV(POLYNATE_CONV avs)) tm in + let SIMP_DNF_CONV = + GEN_REWRITE_CONV TOP_DEPTH_CONV (basic_rewrites()) THENC + NNF_COND_CONV THENC DNF_CONV in + let DISTRIB_EXISTS_CONV = GEN_REWRITE_CONV I [EXISTS_OR_THM] in + let TRIV_EXISTS_CONV = GEN_REWRITE_CONV I [EXISTS_SIMP] in + let complex_ty = `:complex` in + let FINAL_SIMP_CONV = + GEN_REWRITE_CONV DEPTH_CONV [CX_INJ] THENC + REAL_RAT_REDUCE_CONV THENC + GEN_REWRITE_CONV TOP_DEPTH_CONV (basic_rewrites()) in + let rec FULL_COMPLEX_QUELIM_CONV avs tm = + if is_forall tm then + let th1 = ELIM_FORALL_CONV tm in + let th2 = FULL_COMPLEX_QUELIM_CONV avs (rand(concl th1)) in + TRANS th1 th2 + else if is_neg tm then + AP_TERM (rator tm) (FULL_COMPLEX_QUELIM_CONV avs (rand tm)) + else if is_conj tm or is_disj tm or is_imp tm or is_iff tm then + let lop,r = dest_comb tm in + let op,l = dest_comb lop in + let thl = FULL_COMPLEX_QUELIM_CONV avs l + and thr = FULL_COMPLEX_QUELIM_CONV avs r in + MK_COMB(AP_TERM(rator(rator tm)) thl,thr) + else if is_exists tm then + let ev,bod = dest_exists tm in + let th0 = FULL_COMPLEX_QUELIM_CONV (ev::avs) bod in + let th1 = MK_EXISTS ev (CONV_RULE(RAND_CONV SIMP_DNF_CONV) th0) in + TRANS th1 (DISTRIB_AND_COMPLEX_QUELIM_CONV (ev::avs) (rand(concl th1))) + else if is_eq tm then + ELIM_EQ_CONV avs tm + else failwith "unexpected type of formula" + and DISTRIB_AND_COMPLEX_QUELIM_CONV avs tm = + try TRIV_EXISTS_CONV tm + with Failure _ -> try + (DISTRIB_EXISTS_CONV THENC + BINOP_CONV (DISTRIB_AND_COMPLEX_QUELIM_CONV avs)) tm + with Failure _ -> COMPLEX_QUELIM_CONV [] avs tm in + fun tm -> + let avs = filter (fun t -> type_of t = complex_ty) (frees tm) in + (FULL_COMPLEX_QUELIM_CONV avs THENC FINAL_SIMP_CONV) tm;; diff --git a/Complex/quelim_examples.ml b/Complex/quelim_examples.ml new file mode 100644 index 0000000..8e5d021 --- /dev/null +++ b/Complex/quelim_examples.ml @@ -0,0 +1,185 @@ +(* ========================================================================= *) +(* Some examples of full complex quantifier elimination. *) +(* ========================================================================= *) + +let th = time prove + (`!x y. (x pow 2 = Cx(&2)) /\ (y pow 2 = Cx(&3)) + ==> ((x * y) pow 2 = Cx(&6))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`!x a. (a pow 2 = Cx(&2)) /\ (x pow 2 + a * x + Cx(&1) = Cx(&0)) + ==> (x pow 4 + Cx(&1) = Cx(&0))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`!a x. (a pow 2 = Cx(&2)) /\ (x pow 2 + a * x + Cx(&1) = Cx(&0)) + ==> (x pow 4 + Cx(&1) = Cx(&0))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`~(?a x y. (a pow 2 = Cx(&2)) /\ + (x pow 2 + a * x + Cx(&1) = Cx(&0)) /\ + (y * (x pow 4 + Cx(&1)) + Cx(&1) = Cx(&0)))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`!x. ?y. x pow 2 = y pow 3`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`!x y z a b. (a + b) * (x - y + z) - (a - b) * (x + y + z) = + Cx(&2) * (b * x + b * z - a * y)`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`!a b. ~(a = b) ==> ?x y. (y * x pow 2 = a) /\ (y * x pow 2 + x = b)`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`!a b c x y. (a * x pow 2 + b * x + c = Cx(&0)) /\ + (a * y pow 2 + b * y + c = Cx(&0)) /\ + ~(x = y) + ==> (a * x * y = c) /\ (a * (x + y) + b = Cx(&0))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +let th = time prove + (`~(!a b c x y. (a * x pow 2 + b * x + c = Cx(&0)) /\ + (a * y pow 2 + b * y + c = Cx(&0)) + ==> (a * x * y = c) /\ (a * (x + y) + b = Cx(&0)))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +(** geometric example from ``Algorithms for Computer Algebra'': + right triangle where perp. bisector of hypotenuse passes through the + right angle is isoseles. + **) + +let th = time prove + (`!y_1 y_2 y_3 y_4. + (y_1 = Cx(&2) * y_3) /\ + (y_2 = Cx(&2) * y_4) /\ + (y_1 * y_3 = y_2 * y_4) + ==> (y_1 pow 2 = y_2 pow 2)`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +(** geometric example: gradient condition for two lines to be non-parallel. + **) + +let th = time prove + (`!a1 b1 c1 a2 b2 c2. + ~(a1 * b2 = a2 * b1) + ==> ?x y. (a1 * x + b1 * y = c1) /\ (a2 * x + b2 * y = c2)`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +(*********** Apparently takes too long + +let th = time prove + (`!a b c x y. (a * x pow 2 + b * x + c = Cx(&0)) /\ + (a * y pow 2 + b * y + c = Cx(&0)) /\ + (!z. (a * z pow 2 + b * z + c = Cx(&0)) + ==> (z = x) \/ (z = y)) + ==> (a * x * y = c) /\ (a * (x + y) + b = Cx(&0))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + +*************) + +(* ------------------------------------------------------------------------- *) +(* Any three points determine a circle. Not true in complex number version! *) +(* ------------------------------------------------------------------------- *) + +(******** And it takes a lot of memory! + +let th = time prove + (`~(!x1 y1 x2 y2 x3 y3. + ?x0 y0. ((x1 - x0) pow 2 + (y1 - y0) pow 2 = + (x2 - x0) pow 2 + (y2 - y0) pow 2) /\ + ((x2 - x0) pow 2 + (y2 - y0) pow 2 = + (x3 - x0) pow 2 + (y3 - y0) pow 2))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + + **************) + +(* ------------------------------------------------------------------------- *) +(* To show we don't need to consider only closed formulas. *) +(* Can eliminate some, then deal with the rest manually and painfully. *) +(* ------------------------------------------------------------------------- *) + +let th = time prove + (`(?x y. (a * x pow 2 + b * x + c = Cx(&0)) /\ + (a * y pow 2 + b * y + c = Cx(&0)) /\ + ~(x = y)) <=> + (a = Cx(&0)) /\ (b = Cx(&0)) /\ (c = Cx(&0)) \/ + ~(a = Cx(&0)) /\ ~(b pow 2 = Cx(&4) * a * c)`, + CONV_TAC(LAND_CONV FULL_COMPLEX_QUELIM_CONV) THEN + REWRITE_TAC[poly; COMPLEX_MUL_RZERO; COMPLEX_ADD_LID; COMPLEX_ADD_RID] THEN + REWRITE_TAC[COMPLEX_ENTIRE; CX_INJ; REAL_OF_NUM_EQ; ARITH] THEN + ASM_CASES_TAC `a = Cx(&0)` THEN + ASM_REWRITE_TAC[COMPLEX_MUL_LZERO; COMPLEX_MUL_RZERO] THENL + [ASM_CASES_TAC `b = Cx(&0)` THEN + ASM_REWRITE_TAC[COMPLEX_MUL_LZERO; COMPLEX_MUL_RZERO]; + ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH + `b * b * c * Cx(--(&1)) + a * c * c * Cx(&4) = + c * (Cx(&4) * a * c - b * b)`] THEN + ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH + `b * b * b * Cx(--(&1)) + a * b * c * Cx (&4) = + b * (Cx(&4) * a * c - b * b)`] THEN + ONCE_REWRITE_TAC[SIMPLE_COMPLEX_ARITH + `b * b * Cx (&1) + a * c * Cx(--(&4)) = + Cx(--(&1)) * (Cx(&4) * a * c - b * b)`] THEN + REWRITE_TAC[COMPLEX_ENTIRE; COMPLEX_SUB_0; CX_INJ] THEN + CONV_TAC REAL_RAT_REDUCE_CONV THEN + ASM_CASES_TAC `b = Cx(&0)` THEN ASM_REWRITE_TAC[] THEN + ASM_CASES_TAC `c = Cx(&0)` THEN ASM_REWRITE_TAC[] THEN + REWRITE_TAC[COMPLEX_POW_2; COMPLEX_MUL_RZERO; COMPLEX_MUL_LZERO] THEN + REWRITE_TAC[EQ_SYM_EQ]]);; + +(* ------------------------------------------------------------------------- *) +(* Do the same thing directly. *) +(* ------------------------------------------------------------------------- *) + +(**** This seems barely feasible + +let th = time prove + (`!a b c. + (?x y. (a * x pow 2 + b * x + c = Cx(&0)) /\ + (a * y pow 2 + b * y + c = Cx(&0)) /\ + ~(x = y)) <=> + (a = Cx(&0)) /\ (b = Cx(&0)) /\ (c = Cx(&0)) \/ + ~(a = Cx(&0)) /\ ~(b pow 2 = Cx(&4) * a * c)`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + + ****) + +(* ------------------------------------------------------------------------- *) +(* More ambitious: determine a unique circle. Also not true over complexes. *) +(* (consider the points (k, k i) where i is the imaginary unit...) *) +(* ------------------------------------------------------------------------- *) + +(********** Takes too long, I think, and too much memory too + +let th = prove + (`~(!x1 y1 x2 y2 x3 y3 x0 y0 x0' y0'. + ((x1 - x0) pow 2 + (y1 - y0) pow 2 = + (x2 - x0) pow 2 + (y2 - y0) pow 2) /\ + ((x2 - x0) pow 2 + (y2 - y0) pow 2 = + (x3 - x0) pow 2 + (y3 - y0) pow 2) /\ + ((x1 - x0') pow 2 + (y1 - y0') pow 2 = + (x2 - x0') pow 2 + (y2 - y0') pow 2) /\ + ((x2 - x0') pow 2 + (y2 - y0') pow 2 = + (x3 - x0') pow 2 + (y3 - y0') pow 2) + ==> (x0 = x0') /\ (y0 = y0'))`, + CONV_TAC FULL_COMPLEX_QUELIM_CONV);; + + *************) + +(* ------------------------------------------------------------------------- *) +(* Side of a triangle in terms of its bisectors; Kapur survey 5.1. *) +(* ------------------------------------------------------------------------- *) + +(************* +let th = time FULL_COMPLEX_QUELIM_CONV + `?b c. (p1 = ai pow 2 * (b + c) pow 2 - c * b * (c + b - a) * (c + b + a)) /\ + (p2 = ae pow 2 * (c - b) pow 2 - c * b * (a + b - c) * (a - b + a)) /\ + (p3 = be pow 2 * (c - a) pow 2 - a * c * (a + b - c) * (c + b - a))`;; + + *************) diff --git a/make.ml b/make.ml new file mode 100644 index 0000000..266ecf3 --- /dev/null +++ b/make.ml @@ -0,0 +1,15 @@ +needs "Library/analysis.ml";; (* Basic real analysis *) +needs "Library/transc.ml";; (* Real transcendental functions *) +needs "Library/floor.ml";; (* Floor and frac functions *) + +needs "Complex/complexnumbers.ml";; (* Basic complex number defs *) +needs "Complex/complex_transc.ml";; (* Complex transcendental functions *) + +needs "Complex/cpoly.ml";; (* Complex polynomials *) +needs "Complex/fundamental.ml";; (* Fundamental theorem of algebra *) +needs "Complex/quelim.ml";; (* Quantifier elimination algorithm *) +needs "Complex/complex_grobner.ml";; (* Grobner bases with HOL proofs *) +needs "Complex/complex_real.ml";; (* Special case of reals *) + +needs "Complex/quelim_examples.ml";; (* Examples of using quantifier elim *) +needs "Complex/grobner_examples.ml";; (* Examples of using Grobner bases *)