(* ========================================================================= *)
(* The Cayley-Hamilton theorem (for real matrices).                          *)
(* ========================================================================= *)

needs "Multivariate/complexes.ml";;

(* ------------------------------------------------------------------------- *)
(* Powers of a square matrix (mpow) and sums of matrices (msum).             *)
(* ------------------------------------------------------------------------- *)

parse_as_infix("mpow",(24,"left"));;


let mpow = define
  `(!A:real^N^N. A mpow 0 = (mat 1 :real^N^N)) /\
   (!A:real^N^N n. A mpow (SUC n) = A ** A mpow n)`;;


let msum = new_definition
 `msum = iterate (matrix_add)`;;


let NEUTRAL_MATRIX_ADD = 
prove
 (`neutral((+):real^N^M->real^N^M->real^N^M) = mat 0`,

  REWRITE_TAC[neutral] THEN MATCH_MP_TAC SELECT_UNIQUE THEN
  MESON_TAC[MATRIX_ADD_RID; MATRIX_ADD_LID]);;


let MONOIDAL_MATRIX_ADD = 
prove
 (`monoidal((+):real^N^M->real^N^M->real^N^M)`,

  REWRITE_TAC[monoidal; NEUTRAL_MATRIX_ADD] THEN
  REWRITE_TAC[MATRIX_ADD_LID; MATRIX_ADD_ASSOC] THEN
  REWRITE_TAC[MATRIX_ADD_SYM]);;


let MSUM_CLAUSES = 
prove
 (`(!(f:A->real^N^M). msum {} f = mat 0) /\
   (!x (f:A->real^N^M) s.
        FINITE(s)
        ==> (msum (x INSERT s) f =
             if x IN s then msum s f else f(x) + msum s f))`,

  REWRITE_TAC[msum; GSYM NEUTRAL_MATRIX_ADD] THEN
  ONCE_REWRITE_TAC[SWAP_FORALL_THM] THEN
  MATCH_MP_TAC ITERATE_CLAUSES THEN REWRITE_TAC[MONOIDAL_MATRIX_ADD]);;


let MSUM_MATRIX_RMUL = 
prove
 (`!(f:A->real^N^M) (A:real^P^N) s.
        FINITE s ==> msum s (\i. f(i) ** A) = msum s f ** A`,

  GEN_TAC THEN GEN_TAC THEN MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
  SIMP_TAC[MSUM_CLAUSES; MATRIX_MUL_LZERO; MATRIX_ADD_RDISTRIB]);;


let MSUM_MATRIX_LMUL = 
prove
 (`!(f:A->real^P^N) (A:real^N^M) s.
        FINITE s ==> msum s (\i. A ** f(i)) = A ** msum s f`,

  GEN_TAC THEN GEN_TAC THEN MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
  SIMP_TAC[MSUM_CLAUSES; MATRIX_MUL_RZERO; MATRIX_ADD_LDISTRIB]);;


let MSUM_ADD = 
prove
 (`!f g s. FINITE s ==> (msum s (\x. f(x) + g(x)) = msum s f + msum s g)`,

  SIMP_TAC[msum; ITERATE_OP; MONOIDAL_MATRIX_ADD]);;


let MSUM_CMUL = 
prove
 (`!(f:A->real^N^M) c s.
        FINITE s ==>  msum s (\i. c %% f(i)) = c %% msum s f`,

  GEN_TAC THEN GEN_TAC THEN MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
  SIMP_TAC[MSUM_CLAUSES; MATRIX_CMUL_ADD_LDISTRIB; MATRIX_CMUL_RZERO]);;


let MSUM_NEG = 
prove
 (`!(f:A->real^N^M) s.
        FINITE s ==>  msum s (\i. --(f(i))) = --(msum s f)`,

  ONCE_REWRITE_TAC[MATRIX_NEG_MINUS1] THEN
  REWRITE_TAC[MSUM_CMUL]);;


let MSUM_SUB = 
prove
 (`!f g s. FINITE s ==> (msum s (\x. f(x) - g(x)) = msum s f - msum s g)`,

  REWRITE_TAC[MATRIX_SUB; MATRIX_NEG_MINUS1] THEN
  SIMP_TAC[MSUM_ADD; MSUM_CMUL]);;


let MSUM_CLAUSES_LEFT = 
prove
 (`!f m n. m <= n ==> msum(m..n) f = f(m) + msum(m+1..n) f`,

  SIMP_TAC[GSYM NUMSEG_LREC; MSUM_CLAUSES; FINITE_NUMSEG; IN_NUMSEG] THEN
  ARITH_TAC);;


let MSUM_SUPPORT = 
prove
 (`!f s. msum (support (+) f s) f = msum s f`,

  SIMP_TAC[msum; iterate; SUPPORT_SUPPORT]);;


let MSUM_EQ = 
prove
 (`!f g s. (!x. x IN s ==> (f x = g x)) ==> (msum s f = msum s g)`,

  REWRITE_TAC[msum] THEN
  MATCH_MP_TAC ITERATE_EQ THEN REWRITE_TAC[MONOIDAL_MATRIX_ADD]);;


let MSUM_RESTRICT_SET = 
prove
 (`!P s f. msum {x:A | x IN s /\ P x} f =
           msum s (\x. if P x then f(x) else mat 0)`,

  ONCE_REWRITE_TAC[GSYM MSUM_SUPPORT] THEN
  REWRITE_TAC[support; NEUTRAL_MATRIX_ADD; IN_ELIM_THM] THEN
  REWRITE_TAC[MESON[] `~((if P x then f x else a) = a) <=> P x /\ ~(f x = a)`;
              GSYM CONJ_ASSOC] THEN
  REPEAT GEN_TAC THEN MATCH_MP_TAC MSUM_EQ THEN SIMP_TAC[IN_ELIM_THM]);;


let MSUM_SING = 
prove
 (`!f x. msum {x} f = f(x)`,


let MSUM_IMAGE = 
prove
 (`!f g s. (!x y. x IN s /\ y IN s /\ (f x = f y) ==> (x = y))
           ==> (msum (IMAGE f s) g = msum s (g o f))`,

  REWRITE_TAC[msum; GSYM NEUTRAL_MATRIX_ADD] THEN
  MATCH_MP_TAC ITERATE_IMAGE THEN REWRITE_TAC[MONOIDAL_MATRIX_ADD]);;


let MSUM_OFFSET = 
prove
 (`!p f m n. msum(m+p..n+p) f = msum(m..n) (\i. f(i + p))`,


let MSUM_COMPONENT = 
prove
 (`!f:A->real^N^M i j s.
        1 <= i /\ i <= dimindex(:M) /\
        1 <= j /\ j <= dimindex(:N) /\
        FINITE s
        ==> (msum s f)$i$j = sum s (\a. f(a)$i$j)`,

  REPLICATE_TAC 3 GEN_TAC THEN
  REWRITE_TAC[IMP_CONJ; RIGHT_FORALL_IMP_THM] THEN
  REPEAT DISCH_TAC THEN
  MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
  SIMP_TAC[MSUM_CLAUSES; SUM_CLAUSES] THEN
  ASM_SIMP_TAC[MATRIX_ADD_COMPONENT; MAT_COMPONENT; COND_ID]);;


let MSUM_CLAUSES_NUMSEG = 
prove
 (`(!m. msum(m..0) f = if m = 0 then f(0) else mat 0) /\
   (!m n. msum(m..SUC n) f = if m <= SUC n then msum(m..n) f + f(SUC n)
                             else msum(m..n) f)`,


let MPOW_ADD = 
prove
 (`!A:real^N^N m n. A mpow (m + n) = A mpow m ** A mpow n`,

  GEN_TAC THEN INDUCT_TAC THEN
  ASM_REWRITE_TAC[ADD_CLAUSES; mpow; MATRIX_MUL_LID] THEN
  REWRITE_TAC[MATRIX_MUL_ASSOC]);;


let MPOW_1 = 
prove
 (`!A:real^N^N. A mpow 1 = A`,

  REWRITE_TAC[num_CONV `1`; mpow] THEN
  REWRITE_TAC[SYM(num_CONV `1`); MATRIX_MUL_RID]);;


let MPOW_SUC = 
prove
 (`!A:real^N^N n. A mpow (SUC n) = A mpow n ** A`,

  REWRITE_TAC[ADD1; MPOW_ADD; MPOW_1]);;

(* ------------------------------------------------------------------------- *)
(* The main lemma underlying the proof.                                      *)
(* ------------------------------------------------------------------------- *)


let MATRIC_POLYFUN_EQ_0 = 
prove
 (`!n A:num->real^N^M.
        (!x. msum(0..n) (\i. x pow i %% A i) = mat 0) <=>
        (!i. i IN 0..n ==> A i = mat 0)`,

  SIMP_TAC[CART_EQ; MSUM_COMPONENT; MAT_COMPONENT; LAMBDA_BETA;
           FINITE_NUMSEG; COND_ID;
           ONCE_REWRITE_RULE[REAL_MUL_SYM] MATRIX_CMUL_COMPONENT] THEN
  REWRITE_TAC[MESON[]
   `(!x i. P i ==> !j. Q j ==> R x i j) <=>
    (!i. P i ==> !j. Q j ==> !x. R x i j)`] THEN
  REWRITE_TAC[REAL_POLYFUN_EQ_0] THEN MESON_TAC[]);;


let MATRIC_POLY_LEMMA = 
prove
 (`!(A:real^N^N) B (C:real^N^N) n.
        (!x. msum (0..n) (\i. (x pow i) %% B i) ** (A - x %% mat 1) = C)
        ==> C = mat 0`,

  SIMP_TAC[GSYM MSUM_MATRIX_RMUL; FINITE_NUMSEG; MATRIX_SUB_LDISTRIB] THEN
  REWRITE_TAC[MATRIX_MUL_RMUL] THEN ONCE_REWRITE_TAC[MATRIX_MUL_LMUL] THEN
  ONCE_REWRITE_TAC[MATRIX_CMUL_ASSOC] THEN
  REWRITE_TAC[GSYM(CONJUNCT2 real_pow)] THEN
  SIMP_TAC[MSUM_SUB; FINITE_NUMSEG] THEN REPEAT STRIP_TAC THEN
  SUBGOAL_THEN
   `!x. msum(0..SUC n)
         (\i. x pow i %% (((if i = 0 then (--C:real^N^N) else mat 0) +
                           (if i <= n then B i ** (A:real^N^N) else mat 0)) -
                          (if i = 0 then mat 0 else B(i - 1) ** mat 1))) =
        mat 0`
  MP_TAC THENL
   [SIMP_TAC[MATRIX_CMUL_SUB_LDISTRIB; MSUM_SUB; FINITE_NUMSEG;
             MATRIX_CMUL_ADD_LDISTRIB; MSUM_ADD] THEN
    ONCE_REWRITE_TAC[COND_RAND] THEN REWRITE_TAC[MATRIX_CMUL_RZERO] THEN
    ONCE_REWRITE_TAC[MESON[]
     `(if p then mat 0 else x) = (if ~p then x else mat 0)`] THEN
    REWRITE_TAC[GSYM MSUM_RESTRICT_SET; IN_NUMSEG] THEN
    REWRITE_TAC[ARITH_RULE `(0 <= i /\ i <= SUC n) /\ i = 0 <=> i = 0`;
      ARITH_RULE `(0 <= i /\ i <= SUC n) /\ i <= n <=> 0 <= i /\ i <= n`;
      ARITH_RULE `(0 <= i /\ i <= SUC n) /\ ~(i = 0) <=>
                  1 <= i /\ i <= SUC n`] THEN
    REWRITE_TAC[SING_GSPEC; GSYM numseg; MSUM_SING; real_pow] THEN
    REWRITE_TAC[MATRIX_CMUL_LID] THEN
    FIRST_X_ASSUM(fun th -> GEN_REWRITE_TAC ONCE_DEPTH_CONV [GSYM th]) THEN
    REWRITE_TAC[MATRIX_NEG_SUB] THEN REWRITE_TAC[MATRIX_SUB; AC MATRIX_ADD_AC
     `(((A:real^N^N) + --B) + B) + C = (--B + B) + A + C`] THEN
    REWRITE_TAC[MATRIX_ADD_LNEG; MATRIX_ADD_LID] THEN
    REWRITE_TAC[num_CONV `1`] THEN REWRITE_TAC[ADD1; MSUM_OFFSET] THEN
    REWRITE_TAC[ADD_CLAUSES; ADD_SUB; MATRIX_ADD_RNEG];
    REWRITE_TAC[MATRIC_POLYFUN_EQ_0; IN_NUMSEG; LE_0] THEN DISCH_TAC THEN
    SUBGOAL_THEN `!i:num. B(n - i) = (mat 0:real^N^N)` MP_TAC THENL
     [MATCH_MP_TAC num_INDUCTION THEN CONJ_TAC THENL
       [FIRST_X_ASSUM(MP_TAC o SPEC `SUC n`) THEN
        REWRITE_TAC[LE_REFL; SUB_0; NOT_SUC; ARITH_RULE `~(SUC n <= n)`] THEN
        REWRITE_TAC[MATRIX_ADD_LID; SUC_SUB1; MATRIX_MUL_RID] THEN
        REWRITE_TAC[MATRIX_SUB_LZERO; MATRIX_NEG_EQ_0];
        X_GEN_TAC `m:num` THEN DISCH_TAC THEN
        DISJ_CASES_TAC(ARITH_RULE `n - SUC m = n - m \/ m < n`) THEN
        ASM_REWRITE_TAC[] THEN
        FIRST_X_ASSUM(MP_TAC o SPEC `n - m:num`) THEN
        ASM_SIMP_TAC[ARITH_RULE `m < n ==> ~(n - m = 0)`;
                     ARITH_RULE `n - m <= SUC n /\ n - m <= n`] THEN
        REWRITE_TAC[MATRIX_MUL_LZERO; MATRIX_ADD_LID; MATRIX_SUB_LZERO] THEN
        REWRITE_TAC[MATRIX_MUL_RID; MATRIX_NEG_EQ_0] THEN
        ASM_SIMP_TAC[ARITH_RULE `n - m - 1 = n - SUC m`]];
      DISCH_THEN(MP_TAC o SPEC `n:num`) THEN REWRITE_TAC[SUB_REFL] THEN
      DISCH_TAC THEN FIRST_X_ASSUM(MP_TAC o SPEC `0`) THEN
      ASM_REWRITE_TAC[LE_0; MATRIX_MUL_LZERO; MATRIX_ADD_RID] THEN
      REWRITE_TAC[MATRIX_SUB_RZERO; MATRIX_NEG_EQ_0]]]);;

(* ------------------------------------------------------------------------- *)
(* Show that cofactor and determinant are n-1 and n degree polynomials.      *)
(* ------------------------------------------------------------------------- *)


let POLYFUN_N_CONST = 
prove
 (`!c n. ?b. !x. c = sum(0..n) (\i. b i * x pow i)`,

  REPEAT STRIP_TAC THEN
  EXISTS_TAC `\i. if i = 0 then c else &0` THEN
  REWRITE_TAC[COND_RAND; COND_RATOR; REAL_MUL_LZERO] THEN
  REWRITE_TAC[GSYM SUM_RESTRICT_SET; IN_NUMSEG] THEN
  REWRITE_TAC[ARITH_RULE `(0 <= i /\ i <= n) /\ i = 0 <=> i = 0`] THEN
  REWRITE_TAC[SING_GSPEC; SUM_SING; real_pow; REAL_MUL_RID]);;


let POLYFUN_N_ADD = 
prove
 (`!f g. (?b. !x. f(x) = sum(0..n) (\i. b i * x pow i)) /\
         (?c. !x. g(x) = sum(0..n) (\i. c i * x pow i))
         ==> ?d. !x. f(x) + g(x) = sum(0..n) (\i. d i * x pow i)`,

  REPEAT STRIP_TAC THEN
  EXISTS_TAC `\i. (b:num->real) i + c i` THEN
  ASM_REWRITE_TAC[SUM_ADD_NUMSEG; REAL_ADD_RDISTRIB]);;


let POLYFUN_N_CMUL = 
prove
 (`!f c. (?b. !x. f(x) = sum(0..n) (\i. b i * x pow i))
         ==> ?b. !x. c * f(x) = sum(0..n) (\i. b i * x pow i)`,

  REPEAT STRIP_TAC THEN
  EXISTS_TAC `\i. c * (b:num->real) i` THEN
  ASM_REWRITE_TAC[SUM_LMUL; GSYM REAL_MUL_ASSOC]);;


let POLYFUN_N_SUM = 
prove
 (`!f s. FINITE s /\
         (!a. a IN s ==> ?b. !x. f x a = sum(0..n) (\i. b i * x pow i))
         ==> ?b. !x. sum s (f x) = sum(0..n) (\i. b i * x pow i)`,

  GEN_TAC THEN REWRITE_TAC[IMP_CONJ] THEN
  MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
  SIMP_TAC[SUM_CLAUSES; FORALL_IN_INSERT; NOT_IN_EMPTY; POLYFUN_N_CONST] THEN
  REPEAT GEN_TAC THEN REPEAT DISCH_TAC THEN
  MATCH_MP_TAC POLYFUN_N_ADD THEN ASM_SIMP_TAC[]);;


let POLYFUN_N_PRODUCT = 
prove
 (`!f s n. FINITE s /\
           (!a:A. a IN s ==> ?c d. !x. f x a = c + d * x) /\ CARD(s) <= n
           ==> ?b. !x. product s (f x) = sum(0..n) (\i. b i * x pow i)`,

  GEN_TAC THEN REWRITE_TAC[IMP_CONJ; RIGHT_FORALL_IMP_THM] THEN
  MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
  SIMP_TAC[PRODUCT_CLAUSES; POLYFUN_N_CONST; FORALL_IN_INSERT] THEN
  REPEAT GEN_TAC THEN DISCH_THEN(fun th ->
    DISCH_THEN(CONJUNCTS_THEN ASSUME_TAC) THEN MP_TAC th) THEN
  ASM_REWRITE_TAC[] THEN STRIP_TAC THEN
  ASM_SIMP_TAC[CARD_CLAUSES] THEN
  INDUCT_TAC THENL [ARITH_TAC; REWRITE_TAC[LE_SUC]] THEN DISCH_TAC THEN
  FIRST_X_ASSUM(MP_TAC o SPEC `n:num`) THEN ASM_REWRITE_TAC[] THEN
  DISCH_THEN(X_CHOOSE_TAC `b:num->real`) THEN
  FIRST_X_ASSUM(X_CHOOSE_THEN `c:real` (X_CHOOSE_TAC `d:real`)) THEN
  ASM_REWRITE_TAC[] THEN
  EXISTS_TAC `\i. (if i <= n then c * b i else &0) +
                  (if ~(i = 0) then d * b(i - 1) else &0)` THEN
  X_GEN_TAC `x:real` THEN
  REWRITE_TAC[REAL_ADD_RDISTRIB; SUM_ADD_NUMSEG] THEN
  REWRITE_TAC[COND_RAND; COND_RATOR; GSYM SUM_LMUL; REAL_MUL_LZERO] THEN
  REWRITE_TAC[GSYM SUM_RESTRICT_SET; IN_NUMSEG] THEN
  REWRITE_TAC[ARITH_RULE
   `((0 <= i /\ i <= SUC n) /\ i <= n <=> 0 <= i /\ i <= n) /\
    ((0 <= i /\ i <= SUC n) /\ ~(i = 0) <=> 1 <= i /\ i <= SUC n)`] THEN
  REWRITE_TAC[GSYM numseg] THEN
  REWRITE_TAC[MESON[num_CONV `1`; ADD1] `1..SUC n = 0+1..n+1`] THEN
  REWRITE_TAC[SUM_OFFSET; ADD_SUB; REAL_POW_ADD] THEN
  BINOP_TAC THEN MATCH_MP_TAC SUM_EQ_NUMSEG THEN REAL_ARITH_TAC);;


let COFACTOR_ENTRY_AS_POLYFUN = 
prove
 (`!A:real^N^N x i j.
        1 <= i /\ i <= dimindex(:N) /\
        1 <= j /\ j <= dimindex(:N)
        ==> ?c. !x. cofactor(A - x %% mat 1)$i$j =
                    sum(0..dimindex(:N)-1) (\i. c(i) * x pow i)`,

  REPEAT STRIP_TAC THEN ASM_SIMP_TAC[cofactor; LAMBDA_BETA; det] THEN
  MATCH_MP_TAC POLYFUN_N_SUM THEN
  SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG; FORALL_IN_GSPEC] THEN
  X_GEN_TAC `p:num->num` THEN DISCH_TAC THEN
  MATCH_MP_TAC POLYFUN_N_CMUL THEN
  SUBGOAL_THEN `1..dimindex(:N) = i INSERT ((1..dimindex(:N)) DELETE i)`
  SUBST1_TAC THENL
   [REWRITE_TAC[EXTENSION; IN_INSERT; IN_DELETE; IN_NUMSEG] THEN ASM_ARITH_TAC;
    SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG]] THEN
  ASM_REWRITE_TAC[IN_DELETE; IN_NUMSEG] THEN
  MATCH_MP_TAC POLYFUN_N_CMUL THEN
  MATCH_MP_TAC POLYFUN_N_PRODUCT THEN
  SIMP_TAC[CARD_DELETE; FINITE_DELETE; FINITE_NUMSEG] THEN
  ASM_REWRITE_TAC[IN_NUMSEG; IN_DELETE; CARD_NUMSEG_1; LE_REFL] THEN
  X_GEN_TAC `k:num` THEN STRIP_TAC THEN
  SUBGOAL_THEN `(p:num->num) k IN 1..dimindex(:N)` MP_TAC THENL
   [ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]; ALL_TAC] THEN
  ASM_SIMP_TAC[IN_NUMSEG; LAMBDA_BETA] THEN STRIP_TAC THEN
  ASM_CASES_TAC `(p:num->num) k = j` THEN ASM_REWRITE_TAC[] THENL
   [REPEAT(EXISTS_TAC `&0`) THEN REAL_ARITH_TAC; ALL_TAC] THEN
  ASM_SIMP_TAC[MATRIX_SUB_COMPONENT; MATRIX_CMUL_COMPONENT; MAT_COMPONENT] THEN
  REWRITE_TAC[REAL_ARITH `a - x * d:real = a + (--d) * x`] THEN MESON_TAC[]);;


let DETERMINANT_AS_POLYFUN = 
prove
 (`!A:real^N^N.
        ?c. !x. det(A - x %% mat 1) =
                sum(0..dimindex(:N)) (\i. c(i) * x pow i)`,

  GEN_TAC THEN REWRITE_TAC[det] THEN
  MATCH_MP_TAC POLYFUN_N_SUM THEN
  SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG; FORALL_IN_GSPEC] THEN
  X_GEN_TAC `p:num->num` THEN DISCH_TAC THEN
  MATCH_MP_TAC POLYFUN_N_CMUL THEN MATCH_MP_TAC POLYFUN_N_PRODUCT THEN
  SIMP_TAC[FINITE_NUMSEG; CARD_NUMSEG_1; LE_REFL; IN_NUMSEG] THEN
  X_GEN_TAC `k:num` THEN STRIP_TAC THEN
  SUBGOAL_THEN `(p:num->num) k IN 1..dimindex(:N)` MP_TAC THENL
   [ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]; ALL_TAC] THEN
  ASM_SIMP_TAC[IN_NUMSEG; LAMBDA_BETA] THEN STRIP_TAC THEN
  ASM_SIMP_TAC[MATRIX_SUB_COMPONENT; MATRIX_CMUL_COMPONENT; MAT_COMPONENT] THEN
  REWRITE_TAC[REAL_ARITH `a - x * d:real = a + (--d) * x`] THEN MESON_TAC[]);;

(* ------------------------------------------------------------------------- *)
(* Hence define characteristic polynomial coefficients.                      *)
(* ------------------------------------------------------------------------- *)


let char_poly = new_specification ["char_poly"]
  (REWRITE_RULE[SKOLEM_THM] DETERMINANT_AS_POLYFUN);;

(* ------------------------------------------------------------------------- *)
(* Now the Cayley-Hamilton proof.                                            *)
(* ------------------------------------------------------------------------- *)


let COFACTOR_AS_MATRIC_POLYNOMIAL = 
prove
 (`!A:real^N^N. ?C.
      !x. cofactor(A - x %% mat 1) =
          msum(0..dimindex(:N)-1) (\i. x pow i %% C i)`,

  GEN_TAC THEN SIMP_TAC[CART_EQ; MSUM_COMPONENT; FINITE_NUMSEG] THEN
  MP_TAC(ISPEC `A:real^N^N` COFACTOR_ENTRY_AS_POLYFUN) THEN
  REWRITE_TAC[IMP_CONJ; RIGHT_FORALL_IMP_THM] THEN
  REWRITE_TAC[IMP_IMP] THEN REWRITE_TAC[LAMBDA_SKOLEM] THEN
  DISCH_THEN(X_CHOOSE_THEN `c:(num->real)^N^N` ASSUME_TAC) THEN
  EXISTS_TAC `(\i. lambda j k. ((c:(num->real)^N^N)$j$k) i):num->real^N^N` THEN
  MAP_EVERY X_GEN_TAC [`x:real`; `i:num`] THEN STRIP_TAC THEN
  X_GEN_TAC `j:num` THEN STRIP_TAC THEN ASM_SIMP_TAC[] THEN
  MATCH_MP_TAC SUM_EQ_NUMSEG THEN REPEAT STRIP_TAC THEN
  ASM_SIMP_TAC[MATRIX_CMUL_COMPONENT; LAMBDA_BETA] THEN REAL_ARITH_TAC);;


let MATRIC_POWER_DIFFERENCE = 
prove
 (`!A:real^N^N x n.
        A mpow (SUC n) - x pow (SUC n) %% mat 1 =
        msum (0..n) (\i. x pow i %% A mpow (n - i)) ** (A - x %% mat 1)`,

  GEN_TAC THEN GEN_TAC THEN INDUCT_TAC THENL
   [REWRITE_TAC[MSUM_CLAUSES_NUMSEG; real_pow; SUB_0; mpow] THEN
    REWRITE_TAC[MATRIX_MUL_RID; MATRIX_CMUL_LID; MATRIX_MUL_LID] THEN
    REWRITE_TAC[REAL_MUL_RID];
    MATCH_MP_TAC EQ_TRANS THEN
    EXISTS_TAC
     `(A mpow SUC n - x pow SUC n %% mat 1) ** A +
      (x pow (SUC n) %% mat 1 :real^N^N) ** (A - x %% mat 1:real^N^N)` THEN
    CONJ_TAC THENL
     [GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [MPOW_SUC] THEN
      REWRITE_TAC[MATRIX_SUB_RDISTRIB; MATRIX_SUB_LDISTRIB] THEN
      REWRITE_TAC[MATRIX_SUB; MATRIX_MUL_LMUL; MATRIX_MUL_LID] THEN
      REWRITE_TAC[GSYM MATRIX_ADD_ASSOC] THEN AP_TERM_TAC THEN
      REWRITE_TAC[MATRIX_ADD_ASSOC; MATRIX_ADD_LNEG; MATRIX_ADD_LID] THEN
      REWRITE_TAC[real_pow; MATRIX_CMUL_ASSOC] THEN REWRITE_TAC[REAL_MUL_AC];

      ASM_REWRITE_TAC[MSUM_CLAUSES_NUMSEG; LE_0] THEN
      REWRITE_TAC[SUB_REFL; mpow; MATRIX_ADD_RDISTRIB] THEN
      AP_THM_TAC THEN AP_TERM_TAC THEN
      SIMP_TAC[GSYM MSUM_MATRIX_RMUL; FINITE_NUMSEG] THEN
      MATCH_MP_TAC MSUM_EQ THEN REWRITE_TAC[FINITE_NUMSEG] THEN
      X_GEN_TAC `i:num` THEN REWRITE_TAC[IN_NUMSEG] THEN STRIP_TAC THEN
      ASM_SIMP_TAC[MATRIX_MUL_LMUL] THEN AP_TERM_TAC THEN
      ASM_SIMP_TAC[ARITH_RULE `i <= n ==> SUC n - i = SUC(n - i)`] THEN
      REWRITE_TAC[MPOW_SUC; GSYM MATRIX_MUL_ASSOC] THEN AP_TERM_TAC THEN
      REWRITE_TAC[MATRIX_SUB_LDISTRIB; MATRIX_SUB_RDISTRIB] THEN
      REWRITE_TAC[MATRIX_MUL_RMUL; MATRIX_MUL_LMUL] THEN
      REWRITE_TAC[MATRIX_MUL_LID; MATRIX_MUL_RID]]]);;


let MATRIC_CHARPOLY_DIFFERENCE = 
prove
 (`!A:real^N^N. ?B.
      !x. msum(0..dimindex(:N)) (\i. char_poly A i %% A mpow i) -
          sum(0..dimindex(:N)) (\i. char_poly A i * x pow i) %% mat 1 =
          msum(0..(dimindex(:N)-1)) (\i. x pow i %% B i) ** (A - x %% mat 1)`,

  GEN_TAC THEN SPEC_TAC(`dimindex(:N)`,`n:num`) THEN
  SPEC_TAC(`char_poly(A:real^N^N)`,`c:num->real`) THEN
  GEN_TAC THEN INDUCT_TAC THEN
  SIMP_TAC[MSUM_CLAUSES_NUMSEG; SUM_CLAUSES_NUMSEG; LE_0] THENL
   [EXISTS_TAC `(\i. mat 0):num->real^N^N` THEN
    CONV_TAC NUM_REDUCE_CONV THEN REWRITE_TAC[MSUM_CLAUSES_NUMSEG] THEN
    REWRITE_TAC[real_pow; MATRIX_MUL_LMUL; MATRIX_MUL_LZERO; mpow;
                REAL_MUL_RID; MATRIX_CMUL_RZERO; MATRIX_SUB_REFL];
    FIRST_X_ASSUM(X_CHOOSE_TAC `B:num->real^N^N`) THEN
    REWRITE_TAC[MATRIX_SUB; MATRIX_NEG_ADD; MATRIX_CMUL_ADD_RDISTRIB] THEN
    ONCE_REWRITE_TAC[AC MATRIX_ADD_AC
     `(A + B) + (C + D):real^N^N = (A + C) + (B + D)`] THEN
    ASM_REWRITE_TAC[GSYM MATRIX_SUB] THEN
    REWRITE_TAC[GSYM MATRIX_CMUL_ASSOC; GSYM MATRIX_CMUL_SUB_LDISTRIB] THEN
    REWRITE_TAC[MATRIC_POWER_DIFFERENCE; SUC_SUB1] THEN
    EXISTS_TAC `(\i. (if i <= n - 1 then B i else mat 0) +
                     c(SUC n) %% A mpow (n - i)):num->real^N^N` THEN
    X_GEN_TAC `x:real` THEN REWRITE_TAC[MATRIX_CMUL_ADD_LDISTRIB] THEN
    SIMP_TAC[MSUM_ADD; FINITE_NUMSEG; MATRIX_ADD_RDISTRIB] THEN
    REWRITE_TAC[GSYM MATRIX_MUL_LMUL] THEN
    BINOP_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THENL
     [REWRITE_TAC[COND_RAND; COND_RATOR; MATRIX_CMUL_RZERO] THEN
      REWRITE_TAC[GSYM MSUM_RESTRICT_SET; IN_NUMSEG] THEN
      REWRITE_TAC[numseg; ARITH_RULE
       `(0 <= i /\ i <= n) /\ i <= n - 1 <=> 0 <= i /\ i <= n - 1`];
      SIMP_TAC[GSYM MSUM_CMUL; FINITE_NUMSEG; MATRIX_CMUL_ASSOC] THEN
      REWRITE_TAC[REAL_MUL_SYM]]]);;


let CAYLEY_HAMILTON = 
prove
 (`!A:real^N^N. msum(0..dimindex(:N)) (\i. char_poly A i %% A mpow i) = mat 0`,

  GEN_TAC THEN MATCH_MP_TAC MATRIC_POLY_LEMMA THEN MATCH_MP_TAC(MESON[]
   `!g. (!x. g x = k) /\ (?a b c. !x. f a b c x = g x)
        ==> ?a b c. !x. f a b c x = k`) THEN
  EXISTS_TAC
   `\x. (msum(0..dimindex(:N)) (\i. char_poly A i %% (A:real^N^N) mpow i) -
         sum(0..dimindex(:N)) (\i. char_poly A i * x pow i) %% mat 1) +
        sum(0..dimindex(:N)) (\i. char_poly A i * x pow i) %% mat 1` THEN
  REWRITE_TAC[] THEN CONJ_TAC THENL
   [REWRITE_TAC[MATRIX_SUB; GSYM MATRIX_ADD_ASSOC; MATRIX_ADD_LNEG] THEN
    REWRITE_TAC[MATRIX_ADD_RID];
    X_CHOOSE_THEN `B:num->real^N^N` (fun th -> ONCE_REWRITE_TAC[th])
     (ISPEC `A:real^N^N` MATRIC_CHARPOLY_DIFFERENCE) THEN
    REWRITE_TAC[GSYM char_poly; GSYM MATRIX_MUL_LEFT_COFACTOR] THEN
    REWRITE_TAC[GSYM MATRIX_ADD_RDISTRIB] THEN
    REWRITE_TAC[GSYM COFACTOR_TRANSP; TRANSP_MATRIX_SUB] THEN
    REWRITE_TAC[TRANSP_MATRIX_CMUL; TRANSP_MAT] THEN
    X_CHOOSE_THEN `C:num->real^N^N` (fun th -> ONCE_REWRITE_TAC[th])
     (ISPEC `transp A:real^N^N` COFACTOR_AS_MATRIC_POLYNOMIAL) THEN
    MAP_EVERY EXISTS_TAC
     [`A:real^N^N`; `(\i. B i + C i):num->real^N^N`; `dimindex(:N)-1`] THEN
    SIMP_TAC[GSYM MSUM_ADD; FINITE_NUMSEG; MATRIX_CMUL_ADD_LDISTRIB]]);;