(* ========================================================================= *)
(* The law of cosines, of sines, and sum of angles of a triangle.            *)
(* ========================================================================= *)

needs "Multivariate/transcendentals.ml";;

prioritize_vector();;

(* ------------------------------------------------------------------------- *)
(* Angle between vectors (always 0 <= angle <= pi).                          *)
(* ------------------------------------------------------------------------- *)

let vangle = new_definition
 `vangle x y = if x = vec 0 \/ y = vec 0 then pi / &2
               else acs((x dot y) / (norm x * norm y))`;;
(* ------------------------------------------------------------------------- *) (* Traditional geometric notion of angle (but always 0 <= theta <= pi). *) (* ------------------------------------------------------------------------- *)
let angle = new_definition
 `angle(a,b,c) = vangle (a - b) (c - b)`;;
(* ------------------------------------------------------------------------- *) (* Lemmas (more than we need for this result). *) (* ------------------------------------------------------------------------- *)
let VANGLE = 
prove (`!x y:real^N. x dot y = norm(x) * norm(y) * cos(vangle x y)`,
REPEAT GEN_TAC THEN REWRITE_TAC[vangle] THEN ASM_CASES_TAC `x:real^N = vec 0` THEN ASM_REWRITE_TAC[DOT_LZERO; NORM_0; REAL_MUL_LZERO] THEN ASM_CASES_TAC `y:real^N = vec 0` THEN ASM_REWRITE_TAC[DOT_RZERO; NORM_0; REAL_MUL_LZERO; REAL_MUL_RZERO] THEN ONCE_REWRITE_TAC[AC REAL_MUL_AC `a * b * c:real = c * a * b`] THEN ASM_SIMP_TAC[GSYM REAL_EQ_LDIV_EQ; REAL_LT_MUL; NORM_POS_LT] THEN MATCH_MP_TAC(GSYM COS_ACS) THEN ASM_SIMP_TAC[REAL_LE_RDIV_EQ; REAL_LE_LDIV_EQ; NORM_POS_LT; REAL_LT_MUL] THEN MP_TAC(SPECL [`x:real^N`; `y:real^N`] NORM_CAUCHY_SCHWARZ_ABS) THEN REAL_ARITH_TAC);;
let VANGLE_RANGE = 
prove (`!x y:real^N. &0 <= vangle x y /\ vangle x y <= pi`,
REPEAT GEN_TAC THEN REWRITE_TAC[vangle] THEN COND_CASES_TAC THENL [MP_TAC PI_POS THEN REAL_ARITH_TAC; ALL_TAC] THEN RULE_ASSUM_TAC(REWRITE_RULE[DE_MORGAN_THM]) THEN MATCH_MP_TAC ACS_BOUNDS THEN ASM_SIMP_TAC[REAL_LE_RDIV_EQ; REAL_LE_LDIV_EQ; REAL_LT_MUL; NORM_POS_LT] THEN MATCH_MP_TAC(REAL_ARITH `abs(x) <= a ==> -- &1 * a <= x /\ x <= &1 * a`) THEN REWRITE_TAC[NORM_CAUCHY_SCHWARZ_ABS]);;
let ORTHOGONAL_VANGLE = 
prove (`!x y:real^N. orthogonal x y <=> vangle x y = pi / &2`,
REPEAT STRIP_TAC THEN REWRITE_TAC[orthogonal; vangle] THEN ASM_CASES_TAC `x:real^N = vec 0` THEN ASM_REWRITE_TAC[DOT_LZERO] THEN ASM_CASES_TAC `y:real^N = vec 0` THEN ASM_REWRITE_TAC[DOT_RZERO] THEN EQ_TAC THENL [SIMP_TAC[real_div; REAL_MUL_LZERO] THEN DISCH_TAC THEN REWRITE_TAC[GSYM real_div; GSYM COS_PI2] THEN MATCH_MP_TAC ACS_COS THEN MP_TAC PI_POS THEN REAL_ARITH_TAC; MP_TAC(SPECL [`x:real^N`; `y:real^N`] NORM_CAUCHY_SCHWARZ_ABS) THEN GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [GSYM REAL_MUL_LID] THEN REWRITE_TAC[GSYM REAL_BOUNDS_LE] THEN ONCE_REWRITE_TAC[GSYM REAL_MUL_LNEG] THEN ASM_SIMP_TAC[GSYM REAL_LE_RDIV_EQ; GSYM REAL_LE_LDIV_EQ; REAL_LT_MUL; NORM_POS_LT] THEN STRIP_TAC THEN DISCH_THEN(MP_TAC o AP_TERM `cos`) THEN ASM_SIMP_TAC[COS_ACS; COS_PI2] THEN REWRITE_TAC[real_div; REAL_ENTIRE; REAL_INV_EQ_0] THEN ASM_REWRITE_TAC[NORM_EQ_0]]);;
let VANGLE_EQ_PI = 
prove (`!x y:real^N. vangle x y = pi ==> norm(x) % y + norm(y) % x = vec 0`,
REPEAT STRIP_TAC THEN MP_TAC(SPECL [`x:real^N`; `y:real^N`] VANGLE) THEN ASM_REWRITE_TAC[COS_PI] THEN STRIP_TAC THEN MP_TAC(ISPECL [`x:real^N`; `--y:real^N`] NORM_CAUCHY_SCHWARZ_EQ) THEN REWRITE_TAC[NORM_NEG; DOT_RNEG; VECTOR_MUL_RNEG] THEN ASM_REWRITE_TAC[REAL_MUL_RNEG; REAL_NEG_NEG; REAL_MUL_RID] THEN VECTOR_ARITH_TAC);;
let ANGLE_EQ_PI = 
prove (`!A B C:real^N. angle(A,B,C) = pi ==> dist(A,C) = dist(A,B) + dist(B,C)`,
REPEAT GEN_TAC THEN REWRITE_TAC[angle] THEN DISCH_THEN(MP_TAC o MATCH_MP VANGLE_EQ_PI) THEN REWRITE_TAC[VECTOR_ARITH `a + x % (b - c) = vec 0 <=> a = x % (c - b)`] THEN GEN_REWRITE_TAC (funpow 3 LAND_CONV) [NORM_SUB] THEN REWRITE_TAC[GSYM NORM_TRIANGLE_EQ] THEN REWRITE_TAC[VECTOR_ARITH `(B - A) + (C - B):real^N = C - A`] THEN REWRITE_TAC[dist; NORM_SUB]);;
let SIN_ANGLE_POS = 
prove (`!A B C. &0 <= sin(angle(A,B,C))`,
SIMP_TAC[SIN_POS_PI_LE; angle; VANGLE_RANGE]);;
let ANGLE = 
prove (`!A B C. (A - C) dot (B - C) = dist(A,C) * dist(B,C) * cos(angle(A,C,B))`,
REWRITE_TAC[angle; dist; GSYM VANGLE]);;
let ANGLE_REFL = 
prove (`!A B. angle(A,A,B) = pi / &2 /\ angle(B,A,A) = pi / &2`,
REWRITE_TAC[angle; vangle; VECTOR_SUB_REFL]);;
let ANGLE_REFL_MID = 
prove (`!A B. ~(A = B) ==> angle(A,B,A) = &0`,
SIMP_TAC[angle; vangle; VECTOR_SUB_EQ; GSYM NORM_POW_2; GSYM REAL_POW_2; REAL_DIV_REFL; ACS_1; REAL_POW_EQ_0; ARITH; NORM_EQ_0]);;
let ANGLE_SYM = 
prove (`!A B C. angle(A,B,C) = angle(C,B,A)`,
REWRITE_TAC[angle; vangle; VECTOR_SUB_EQ; DISJ_SYM; REAL_MUL_SYM; DOT_SYM]);;
let ANGLE_RANGE = 
prove (`!A B C. &0 <= angle(A,B,C) /\ angle(A,B,C) <= pi`,
REWRITE_TAC[angle; VANGLE_RANGE]);;
(* ------------------------------------------------------------------------- *) (* The law of cosines. *) (* ------------------------------------------------------------------------- *)
let LAW_OF_COSINES = 
prove (`!A B C:real^N. dist(B,C) pow 2 = dist(A,B) pow 2 + dist(A,C) pow 2 - &2 * dist(A,B) * dist(A,C) * cos(angle(B,A,C))`,
REPEAT GEN_TAC THEN REWRITE_TAC[angle; ONCE_REWRITE_RULE[NORM_SUB] dist; GSYM VANGLE; NORM_POW_2] THEN VECTOR_ARITH_TAC);;
(* ------------------------------------------------------------------------- *) (* The law of sines. *) (* ------------------------------------------------------------------------- *)
let LAW_OF_SINES = 
prove (`!A B C:real^N. sin(angle(A,B,C)) * dist(B,C) = sin(angle(B,A,C)) * dist(A,C)`,
REPEAT GEN_TAC THEN MATCH_MP_TAC REAL_POW_EQ THEN EXISTS_TAC `2` THEN SIMP_TAC[SIN_ANGLE_POS; DIST_POS_LE; REAL_LE_MUL; ARITH] THEN REWRITE_TAC[REAL_POW_MUL; MATCH_MP (REAL_ARITH `x + y = &1 ==> x = &1 - y`) (SPEC_ALL SIN_CIRCLE)] THEN ASM_CASES_TAC `A:real^N = B` THEN ASM_REWRITE_TAC[ANGLE_REFL; COS_PI2] THEN RULE_ASSUM_TAC(ONCE_REWRITE_RULE[GSYM VECTOR_SUB_EQ]) THEN RULE_ASSUM_TAC(REWRITE_RULE[GSYM NORM_EQ_0]) THEN FIRST_X_ASSUM(MATCH_MP_TAC o MATCH_MP (REAL_RING `~(a = &0) ==> a pow 2 * x = a pow 2 * y ==> x = y`)) THEN ONCE_REWRITE_TAC[DIST_SYM] THEN REWRITE_TAC[GSYM dist] THEN GEN_REWRITE_TAC (RAND_CONV o LAND_CONV o ONCE_DEPTH_CONV) [DIST_SYM] THEN REWRITE_TAC[REAL_RING `a * (&1 - x) * b = c * (&1 - y) * d <=> a * b - a * b * x = c * d - c * d * y`] THEN REWRITE_TAC[GSYM REAL_POW_MUL; GSYM ANGLE] THEN REWRITE_TAC[REAL_POW_MUL; dist; NORM_POW_2] THEN REWRITE_TAC[DOT_LSUB; DOT_RSUB; DOT_SYM] THEN CONV_TAC REAL_RING);;
(* ------------------------------------------------------------------------- *) (* Hence the sum of the angles of a triangle. *) (* ------------------------------------------------------------------------- *)
let TRIANGLE_ANGLE_SUM_LEMMA = 
prove (`!A B C:real^N. ~(A = B) /\ ~(A = C) /\ ~(B = C) ==> cos(angle(B,A,C) + angle(A,B,C) + angle(B,C,A)) = -- &1`,
REPEAT GEN_TAC THEN ONCE_REWRITE_TAC[GSYM VECTOR_SUB_EQ] THEN REWRITE_TAC[GSYM NORM_EQ_0] THEN MP_TAC(ISPECL [`A:real^N`; `B:real^N`; `C:real^N`] LAW_OF_COSINES) THEN MP_TAC(ISPECL [`B:real^N`; `A:real^N`; `C:real^N`] LAW_OF_COSINES) THEN MP_TAC(ISPECL [`C:real^N`; `B:real^N`; `A:real^N`] LAW_OF_COSINES) THEN MP_TAC(ISPECL [`A:real^N`; `B:real^N`; `C:real^N`] LAW_OF_SINES) THEN MP_TAC(ISPECL [`B:real^N`; `A:real^N`; `C:real^N`] LAW_OF_SINES) THEN MP_TAC(ISPECL [`B:real^N`; `C:real^N`; `A:real^N`] LAW_OF_SINES) THEN REWRITE_TAC[COS_ADD; SIN_ADD; dist; NORM_SUB] THEN MAP_EVERY (fun t -> MP_TAC(SPEC t SIN_CIRCLE)) [`angle(B:real^N,A,C)`; `angle(A:real^N,B,C)`; `angle(B:real^N,C,A)`] THEN REWRITE_TAC[COS_ADD; SIN_ADD; ANGLE_SYM] THEN CONV_TAC REAL_RING);;
let COS_MINUS1_LEMMA = 
prove (`!x. cos(x) = -- &1 /\ &0 <= x /\ x < &3 * pi ==> x = pi`,
REPEAT STRIP_TAC THEN SUBGOAL_THEN `?n. integer n /\ x = n * pi` (X_CHOOSE_THEN `nn:real` (CONJUNCTS_THEN2 ASSUME_TAC SUBST_ALL_TAC)) THEN REWRITE_TAC[GSYM SIN_EQ_0] THENL [MP_TAC(SPEC `x:real` SIN_CIRCLE) THEN ASM_REWRITE_TAC[] THEN CONV_TAC REAL_RING; ALL_TAC] THEN SUBGOAL_THEN `?n. nn = &n` (X_CHOOSE_THEN `n:num` SUBST_ALL_TAC) THENL [FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [REAL_MUL_POS_LE]) THEN SIMP_TAC[PI_POS; REAL_ARITH `&0 < p ==> ~(p < &0) /\ ~(p = &0)`] THEN ASM_MESON_TAC[INTEGER_POS; REAL_LT_LE]; ALL_TAC] THEN MATCH_MP_TAC(REAL_RING `n = &1 ==> n * p = p`) THEN REWRITE_TAC[REAL_OF_NUM_EQ] THEN MATCH_MP_TAC(ARITH_RULE `n < 3 /\ ~(n = 0) /\ ~(n = 2) ==> n = 1`) THEN RULE_ASSUM_TAC(SIMP_RULE[REAL_LT_RMUL_EQ; PI_POS; REAL_OF_NUM_LT]) THEN ASM_REWRITE_TAC[] THEN CONJ_TAC THEN DISCH_THEN SUBST_ALL_TAC THEN REPEAT(POP_ASSUM MP_TAC) THEN SIMP_TAC[COS_0; REAL_MUL_LZERO; COS_NPI] THEN REAL_ARITH_TAC);;
let TRIANGLE_ANGLE_SUM = 
prove (`!A B C:real^N. ~(A = B) /\ ~(A = C) /\ ~(B = C) ==> angle(B,A,C) + angle(A,B,C) + angle(B,C,A) = pi`,
REPEAT STRIP_TAC THEN MATCH_MP_TAC COS_MINUS1_LEMMA THEN ASM_SIMP_TAC[TRIANGLE_ANGLE_SUM_LEMMA; REAL_LE_ADD; ANGLE_RANGE] THEN MATCH_MP_TAC(REAL_ARITH `&0 <= x /\ x <= p /\ &0 <= y /\ y <= p /\ &0 <= z /\ z <= p /\ ~(x = p /\ y = p /\ z = p) ==> x + y + z < &3 * p`) THEN ASM_SIMP_TAC[ANGLE_RANGE] THEN REPEAT STRIP_TAC THEN REPEAT(FIRST_X_ASSUM(MP_TAC o MATCH_MP ANGLE_EQ_PI)) THEN REPEAT(FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV [GSYM VECTOR_SUB_EQ])) THEN REWRITE_TAC[GSYM NORM_EQ_0; dist; NORM_SUB] THEN REAL_ARITH_TAC);;