(* ========================================================================= *)
(* Simple universal variant of Bob Solovay's procedure for vector spaces.    *)
(* ========================================================================= *)

needs "Multivariate/misc.ml";;
needs "Multivariate/vectors.ml";;

(* ------------------------------------------------------------------------- *)
(* Initial simplification so we just use dot products between vectors.       *)
(* ------------------------------------------------------------------------- *)

let VECTOR_SUB_ELIM_THM = 
prove (`(--x = --(&1) % x) /\ (x - y = x + --(&1) % y)`,
VECTOR_ARITH_TAC);;
let NORM_ELIM_THM = 
prove (`!P t. P (norm t) = !x. &0 <= x /\ (x pow 2 = (t:real^N) dot t) ==> P x`,
GEN_TAC THEN REWRITE_TAC[vector_norm] THEN MESON_TAC[DOT_POS_LE; SQRT_POW2; SQRT_UNIQUE; REAL_POW_2; REAL_POW2_ABS; REAL_ABS_POS]);;
let NORM_ELIM_CONV = let dest_norm tm = let nm,v = dest_comb tm in if fst(dest_const nm) <> "vector_norm" then failwith "dest_norm" else v in let is_norm = can dest_norm in fun tm -> let t = find_term (fun t -> is_norm t & free_in t tm) tm in let v = dest_norm t in let w = genvar(type_of t) in let th1 = ISPECL [mk_abs(w,subst[w,t] tm); v] NORM_ELIM_THM in CONV_RULE(COMB2_CONV (RAND_CONV BETA_CONV) (BINDER_CONV(RAND_CONV BETA_CONV))) th1;; let NORM_ELIM_TAC = CONV_TAC NORM_ELIM_CONV THEN GEN_TAC;; let SOLOVAY_TAC = REWRITE_TAC[orthogonal; GSYM DOT_EQ_0] THEN REWRITE_TAC[VECTOR_EQ] THEN REWRITE_TAC[VECTOR_SUB_ELIM_THM] THEN REWRITE_TAC[NORM_EQ; NORM_LE; NORM_LT; real_gt; real_ge] THEN REPEAT NORM_ELIM_TAC THEN REWRITE_TAC[DOT_LADD; DOT_RADD; DOT_LMUL; DOT_RMUL];; (* ------------------------------------------------------------------------- *) (* Iterative Gram-Schmidt type process. *) (* ------------------------------------------------------------------------- *)
let component = new_definition
  `component (b:real^N) x = (b dot x) / (b dot b)`;;
let COMPONENT_ORTHOGONAL = 
prove (`!b:real^N x. orthogonal b (x - (component b x) % b)`,
REPEAT GEN_TAC THEN ASM_CASES_TAC `b = vec 0 :real^N` THENL [ASM_REWRITE_TAC[orthogonal; DOT_LZERO]; ALL_TAC] THEN ASM_SIMP_TAC[orthogonal; component] THEN REWRITE_TAC[DOT_RSUB; DOT_RMUL] THEN ASM_SIMP_TAC[REAL_SUB_REFL; REAL_DIV_RMUL; DOT_EQ_0]);;
let ORTHOGONAL_SUM_LEMMA = 
prove (`!cs vs. ALL (orthogonal x) vs /\ orthogonal x z /\ (LENGTH cs = LENGTH vs) ==> orthogonal x (ITLIST2 (\a v s. a % v + s) cs vs z)`,
LIST_INDUCT_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[NOT_CONS_NIL; NOT_SUC; ITLIST2; LENGTH; ALL] THEN ASM_SIMP_TAC[ORTHOGONAL_CLAUSES; SUC_INJ]);;
let GRAM_SCHMIDT_LEMMA = 
prove (`!w:real^N vs. ?u as. ALL (orthogonal u) vs /\ (LENGTH as = LENGTH vs) /\ (w = ITLIST2 (\a v s. a % v + s) as vs u)`,
ONCE_REWRITE_TAC[SWAP_FORALL_THM] THEN MATCH_MP_TAC list_INDUCT THEN SIMP_TAC[ALL; LENGTH; ITLIST2; LENGTH_EQ_NIL] THEN CONJ_TAC THENL [X_GEN_TAC `w:real^N` THEN EXISTS_TAC `w:real^N` THEN EXISTS_TAC `[]:real list` THEN REWRITE_TAC[ITLIST2]; ALL_TAC] THEN MAP_EVERY X_GEN_TAC [`v:real^N`; `vs:(real^N)list`] THEN REWRITE_TAC[LENGTH_EQ_CONS] THEN DISCH_TAC THEN X_GEN_TAC `w:real^N` THEN FIRST_X_ASSUM(fun th -> MP_TAC(SPEC `w:real^N` th) THEN MP_TAC(SPEC `v:real^N` th)) THEN DISCH_THEN(X_CHOOSE_THEN `z:real^N` (X_CHOOSE_THEN `cs:real list` (STRIP_ASSUME_TAC o GSYM))) THEN DISCH_THEN(X_CHOOSE_THEN `u:real^N` (X_CHOOSE_THEN `as:real list` (STRIP_ASSUME_TAC o GSYM))) THEN MP_TAC(ISPECL [`z:real^N`; `u:real^N`] COMPONENT_ORTHOGONAL) THEN ABBREV_TAC `k = component z (u:real^N)` THEN ABBREV_TAC `x = u - k % z :real^N` THEN DISCH_TAC THEN MAP_EVERY EXISTS_TAC [`x:real^N`; `CONS k (MAP2 (\a c. a - k * c) as cs)`] THEN REWRITE_TAC[CONS_11; RIGHT_EXISTS_AND_THM; GSYM CONJ_ASSOC; UNWIND_THM1] THEN SUBGOAL_THEN `ALL (orthogonal(x:real^N)) vs` ASSUME_TAC THENL [UNDISCH_TAC `ALL (orthogonal(z:real^N)) vs` THEN UNDISCH_TAC `ALL (orthogonal(u:real^N)) vs` THEN REWRITE_TAC[IMP_IMP; AND_ALL] THEN MATCH_MP_TAC MONO_ALL THEN REWRITE_TAC[] THEN EXPAND_TAC "x" THEN SIMP_TAC[ORTHOGONAL_CLAUSES]; ALL_TAC] THEN REPEAT CONJ_TAC THENL [EXPAND_TAC "v" THEN MATCH_MP_TAC ORTHOGONAL_SUM_LEMMA THEN ASM_REWRITE_TAC[] THEN ASM_MESON_TAC[ORTHOGONAL_SYM]; FIRST_ASSUM ACCEPT_TAC; ASM_MESON_TAC[LENGTH_MAP2]; ALL_TAC] THEN REWRITE_TAC[ITLIST2; VECTOR_ARITH `(a = b + c:real^N) = (c = a - b)`] THEN MAP_EVERY EXPAND_TAC ["v";
"w"; "x"] THEN UNDISCH_TAC `LENGTH(vs:(real^N)list) = LENGTH(cs:real list)` THEN UNDISCH_TAC `LENGTH(vs:(real^N)list) = LENGTH(as:real list)` THEN REWRITE_TAC[IMP_CONJ] THEN MAP_EVERY (fun v -> SPEC_TAC(v,v)) [`vs:(real^N)list`; `cs:real list`; `as:real list`] THEN POP_ASSUM_LIST(K ALL_TAC) THEN LIST_INDUCT_TAC THEN LIST_INDUCT_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[NOT_CONS_NIL; NOT_SUC; ITLIST2; LENGTH; ALL; SUC_INJ; MAP2] THEN ASM_SIMP_TAC[] THEN REPEAT DISCH_TAC THEN VECTOR_ARITH_TAC);; (* ------------------------------------------------------------------------- *) (* Hence this is a simple equality. *) (* ------------------------------------------------------------------------- *)
let SOLOVAY_LEMMA = 
prove (`!P vs. (!w:real^N. P w vs) = (!as u. ALL (orthogonal u) vs /\ (LENGTH as = LENGTH vs) ==> P (ITLIST2 (\a v s. a % v + s) as vs u) vs)`,
REPEAT GEN_TAC THEN EQ_TAC THENL [MESON_TAC[]; ALL_TAC] THEN DISCH_TAC THEN X_GEN_TAC `w:real^N` THEN MP_TAC(ISPECL [`w:real^N`; `vs:(real^N)list`] GRAM_SCHMIDT_LEMMA) THEN ASM_MESON_TAC[]);;
(* ------------------------------------------------------------------------- *) (* Set up the specific instances to get rid of list stuff. *) (* ------------------------------------------------------------------------- *)
let FORALL_LENGTH_CLAUSES = 
prove (`((!l. (LENGTH l = 0) ==> P l) = P []) /\ ((!l. (LENGTH l = SUC n) ==> P l) = (!h t. (LENGTH t = n) ==> P (CONS h t)))`,
let ORTHOGONAL_SIMP_CLAUSES = 
prove (`orthogonal u x ==> (u dot x = &0) /\ (x dot u = &0) /\ (u dot (a % x) = &0) /\ ((a % x) dot u = &0) /\ (u dot (a % x + y) = u dot y) /\ ((a % x + y) dot u = y dot u) /\ (u dot (y + a % x) = u dot y) /\ ((y + a % x) dot u = y dot u)`,
SIMP_TAC[orthogonal; DOT_LADD; DOT_RADD; DOT_LMUL; DOT_RMUL] THEN GEN_REWRITE_TAC (LAND_CONV o LAND_CONV) [DOT_SYM] THEN SIMP_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO; REAL_ADD_LID; REAL_ADD_RID]);;
(* ------------------------------------------------------------------------- *) (* A nicer proforma version. *) (* ------------------------------------------------------------------------- *)
let ITLIST2_0_LEMMA = 
prove (`!u as vs. ITLIST2 (\a v s. a % v + s) as vs u = ITLIST2 (\a v s. a % v + s) as vs (vec 0) + u`,
GEN_TAC THEN LIST_INDUCT_TAC THEN REWRITE_TAC[ITLIST2_DEF; VECTOR_ADD_LID] THEN ASM_REWRITE_TAC[VECTOR_ADD_ASSOC]);;
let SOLOVAY_PROFORMA_EQ = 
prove (`(!w:real^N. P (MAP ((dot) w) (CONS w vs)) vs) = (!u. ALL (orthogonal u) vs ==> !as. (LENGTH as = LENGTH vs) ==> P (CONS ((ITLIST2 (\a v s. a % v + s) as vs (vec 0)) dot (ITLIST2 (\a v s. a % v + s) as vs (vec 0)) + u dot u) (MAP ((dot) (ITLIST2 (\a v s. a % v + s) as vs (vec 0))) vs)) vs)`,
MP_TAC(ISPEC `\w:real^N vs. P (MAP ((dot) w) (CONS w vs)) vs :bool` SOLOVAY_LEMMA) THEN REWRITE_TAC[] THEN DISCH_THEN(fun th -> GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [th]) THEN GEN_REWRITE_TAC LAND_CONV [SWAP_FORALL_THM] THEN AP_TERM_TAC THEN GEN_REWRITE_TAC I [FUN_EQ_THM] THEN X_GEN_TAC `u:real^N` THEN REWRITE_TAC[] THEN REWRITE_TAC[RIGHT_IMP_FORALL_THM] THEN AP_TERM_TAC THEN GEN_REWRITE_TAC I [FUN_EQ_THM] THEN X_GEN_TAC `as:(real)list` THEN REWRITE_TAC[IMP_IMP] THEN MATCH_MP_TAC(TAUT `(a ==> (b = c)) ==> (a ==> b <=> a ==> c)`) THEN STRIP_TAC THEN REWRITE_TAC[MAP] THEN BINOP_TAC THEN REWRITE_TAC[CONS_11] THEN ONCE_REWRITE_TAC[ITLIST2_0_LEMMA] THEN REWRITE_TAC[VECTOR_ADD_RID] THEN REWRITE_TAC[VECTOR_ARITH `(a + u) dot (a + u) = a dot a + &2 * (u dot a) + u dot u`] THEN REWRITE_TAC[REAL_ARITH `(a + &2 * b + c = a + c) <=> (b = &0)`] THEN GEN_REWRITE_TAC (RAND_CONV o BINOP_CONV o LAND_CONV) [GSYM ETA_AX] THEN REWRITE_TAC[DOT_LADD] THEN CONJ_TAC THENL [POP_ASSUM_LIST(MP_TAC o end_itlist CONJ) THEN SPEC_TAC(`vs:(real^N)list`,`vs:(real^N)list`) THEN SPEC_TAC(`as:(real)list`,`as:(real)list`) THEN REPEAT LIST_INDUCT_TAC THEN REWRITE_TAC[LENGTH; NOT_SUC] THEN REWRITE_TAC[ALL; ITLIST2; DOT_RZERO; SUC_INJ] THEN ASM_SIMP_TAC[DOT_RADD] THEN REWRITE_TAC[REAL_ADD_RID; DOT_RMUL] THEN SIMP_TAC[orthogonal] THEN REWRITE_TAC[REAL_MUL_RZERO]; MATCH_MP_TAC MAP_EQ THEN REWRITE_TAC[REAL_ARITH `(a + b = a) <=> (b = &0)`] THEN MATCH_MP_TAC ALL_IMP THEN EXISTS_TAC `orthogonal (u:real^N)` THEN ASM_REWRITE_TAC[] THEN SIMP_TAC[orthogonal]]);;
(* ------------------------------------------------------------------------- *) (* The implication that we normally use. *) (* ------------------------------------------------------------------------- *)
let SOLOVAY_PROFORMA = 
prove (`!P vs. (!c. &0 <= c ==> !as. (LENGTH as = LENGTH vs) ==> P (CONS ((ITLIST2 (\a v s. a % v + s) as vs (vec 0)) dot (ITLIST2 (\a v s. a % v + s) as vs (vec 0)) + c) (MAP ((dot) (ITLIST2 (\a v s. a % v + s) as vs (vec 0))) vs)) vs) ==> !w:real^N. P (MAP ((dot) w) (CONS w vs)) vs`,
REPEAT GEN_TAC THEN GEN_REWRITE_TAC RAND_CONV [SOLOVAY_PROFORMA_EQ] THEN REWRITE_TAC[RIGHT_IMP_FORALL_THM; IMP_IMP] THEN REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[DOT_POS_LE]);;
(* ------------------------------------------------------------------------- *) (* Automatically set up an implication for n (+1 eliminated) quantifier. *) (* ------------------------------------------------------------------------- *) let SOLOVAY_RULE = let v_tm = `v:(real^N)list` and d_tm = `d:real list` and elv_tm = `EL:num->(real^N)list->real^N` and eld_tm = `EL:num->(real)list->real` and rn_ty = `:real^N` and rewr_rule = REWRITE_RULE [MAP; EL; HD; TL; LENGTH; FORALL_LENGTH_CLAUSES; ITLIST2; VECTOR_ADD_RID; VECTOR_ADD_LID; DOT_LZERO] and sewr_rule = PURE_ONCE_REWRITE_RULE[DOT_SYM] in fun n -> let args = map (fun i -> mk_comb(mk_comb(elv_tm,mk_small_numeral i),v_tm)) (0--(n-1)) @ map (fun i -> mk_comb(mk_comb(eld_tm,mk_small_numeral i),d_tm)) (1--n) @ [mk_comb(mk_comb(eld_tm,mk_small_numeral 0),d_tm)] in let pty = itlist (mk_fun_ty o type_of) args bool_ty in let p_tm = list_mk_abs([d_tm;v_tm],list_mk_comb(mk_var("P",pty),args)) and vs = make_args "v" [] (replicate rn_ty n) in let th1 = ISPECL [p_tm; mk_list(vs,rn_ty)] SOLOVAY_PROFORMA in let th2 = rewr_rule(CONV_RULE(TOP_DEPTH_CONV num_CONV) th1) in let th3 = sewr_rule th2 in itlist (fun v -> MATCH_MP MONO_FORALL o GEN v) vs th3;; (* ------------------------------------------------------------------------- *) (* Now instantiate it to some special cases. *) (* ------------------------------------------------------------------------- *) let MK_SOLOVAY_PROFORMA = let preths = map SOLOVAY_RULE (0--9) in fun n -> if n < 10 then el n preths else SOLOVAY_RULE n;; (* ------------------------------------------------------------------------- *) (* Apply it to a goal. *) (* ------------------------------------------------------------------------- *) let is_vector_ty ty = match ty with Tyapp("cart",[Tyapp("real",[]);_]) -> true | _ -> false;; let SOLOVAY_REDUCE_TAC (asl,w) = let avs = sort (<) (filter (is_vector_ty o type_of) (frees w)) in (REWRITE_TAC[DOT_SYM] THEN MAP_EVERY (fun v -> SPEC_TAC(v,v)) (rev avs) THEN MATCH_MP_TAC(MK_SOLOVAY_PROFORMA (length avs - 1)) THEN REWRITE_TAC[DOT_LADD; DOT_LMUL; DOT_RADD; DOT_RMUL; DOT_LZERO; DOT_RZERO] THEN REPEAT GEN_TAC) (asl,w);; (* ------------------------------------------------------------------------- *) (* Overall tactic. *) (* ------------------------------------------------------------------------- *) let SOLOVAY_VECTOR_TAC = REWRITE_TAC[dist; real_gt; real_ge; NORM_LT; NORM_LE; GSYM DOT_POS_LT] THEN REPEAT GEN_TAC THEN SOLOVAY_TAC THEN REWRITE_TAC[DOT_LZERO; DOT_RZERO] THEN REPEAT SOLOVAY_REDUCE_TAC THEN REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO; REAL_MUL_LID; REAL_MUL_RID; REAL_ADD_LID; REAL_ADD_RID] THEN REWRITE_TAC[IMP_IMP; RIGHT_IMP_FORALL_THM; GSYM CONJ_ASSOC] THEN REPEAT GEN_TAC THEN REWRITE_TAC[REAL_MUL_LNEG; REAL_MUL_RNEG; REAL_NEG_NEG; REAL_MUL_LID; REAL_MUL_RID; GSYM real_sub];; (* ------------------------------------------------------------------------- *) (* An example where REAL_RING then works. *) (* ------------------------------------------------------------------------- *)
let PYTHAGORAS = 
prove (`!A B C:real^N. orthogonal (A - B) (C - B) ==> norm(C - A) pow 2 = norm(B - A) pow 2 + norm(C - B) pow 2`,
SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_RING);;
(*** Actually in this case we can fairly easily do things manually, though we do need to explicitly use symmetry of the dot product. let PYTHAGORAS = prove (`!A B C:real^N. orthogonal (A - B) (C - B) ==> norm(C - A) pow 2 = norm(B - A) pow 2 + norm(C - B) pow 2`, REWRITE_TAC[NORM_POW_2; orthogonal; DOT_LSUB; DOT_RSUB; DOT_SYM] THEN CONV_TAC REAL_RING);; ***) (* ------------------------------------------------------------------------- *) (* Examples. *) (* ------------------------------------------------------------------------- *) needs "Examples/sos.ml";;
let EXAMPLE_1 = 
prove (`!x y:real^N. x dot y <= norm x * norm y`,
SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);;
let EXAMPLE_2 = 
prove (`!x y:real^N. a % (x + y) = a % x + a % y`,
SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);;
(*** Takes a few minutes but does work let EXAMPLE_3 = prove (`!x y:real^N. norm (x + y) <= norm x + norm y`, SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);; ****)
let EXAMPLE_4 = 
prove (`!x y z. x dot (y + z) = (x dot y) + (x dot z)`,
SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);;
let EXAMPLE_5 = 
prove (`!x y. (x dot x = &0) ==> (x dot y = &0)`,
SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);;
(* ------------------------------------------------------------------------- *) (* This is NORM_INCREASES_ONLINE. *) (* ------------------------------------------------------------------------- *) g `!a d:real^N. ~(d = vec 0) ==> norm (a + d) > norm a \/ norm (a - d) > norm a`;; time e SOLOVAY_VECTOR_TAC;; time e (CONV_TAC REAL_SOS);; (* ------------------------------------------------------------------------- *) (* DIST_INCREASES_ONLINE *) (* ------------------------------------------------------------------------- *) g `!b a d:real^N. ~(d = vec 0) ==> dist(a,b + d) > dist(a,b) \/ dist(a,b - d) > dist(a,b)`;; time e SOLOVAY_VECTOR_TAC;; time e (CONV_TAC REAL_SOS);; (* ------------------------------------------------------------------------- *) (* This one doesn't seem to work easily, but I think it does eventually. *) (* ------------------------------------------------------------------------- *) (**** let EXAMPLE_6 = prove (`!a x. norm(a % x) = abs(a) * norm x`;; SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);; ****)
let EXAMPLE_7 = 
prove (`!x. abs(norm x) = norm x`,
SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);;
(*** But this is (at least) really slow let EXAMPLE_8 = prove (`!x y. abs(norm(x) - norm(y)) <= abs(norm(x - y))`, SOLOVAY_VECTOR_TAC THEN CONV_TAC REAL_SOS);; ****) (* ------------------------------------------------------------------------- *) (* One from separating hyperplanes with a richer structure. *) (* ------------------------------------------------------------------------- *) needs "Rqe/make.ml";;
let EXAMPLE_9 = 
prove (`!x:real^N y. x dot y > &0 ==> ?u. &0 < u /\ norm(u % y - x) < norm x`,
SOLOVAY_VECTOR_TAC THEN W(fun (asl,w) -> MAP_EVERY (fun v -> SPEC_TAC(v,v)) (frees w)) THEN CONV_TAC REAL_QELIM_CONV);;
(* ------------------------------------------------------------------------- *) (* Even richer set of quantifier alternations. *) (* ------------------------------------------------------------------------- *)
let EXAMPLE_10 = 
prove (`!x:real^N y. x dot y > &0 ==> ?u. &0 < u /\ !v. &0 < v /\ v <= u ==> norm(v % y - x) < norm x`,
SOLOVAY_VECTOR_TAC THEN W(fun (asl,w) -> MAP_EVERY (fun v -> SPEC_TAC(v,v)) (frees w)) THEN CONV_TAC REAL_QELIM_CONV);;