colo12-c703.uibk.ac.at Git - Multivariate Analysis/.git/blob - Multivariate/determinants.ml

   1 (* ========================================================================= *)
   2 (* Determinant and trace of a square matrix.                                 *)
   3 (*                                                                           *)
   4 (*              (c) Copyright, John Harrison 1998-2008                       *)
   5 (* ========================================================================= *)
   6
   7 needs "Multivariate/vectors.ml";;
   8 needs "Library/permutations.ml";;
   9 needs "Library/floor.ml";;
  10 needs "Library/products.ml";;
  11
  12 prioritize_real();;
  13
  14 (* ------------------------------------------------------------------------- *)
  15 (* Trace of a matrix (this is relatively easy).                              *)
  16 (* ------------------------------------------------------------------------- *)
  17
  18 let trace = new_definition
  19   `(trace:real^N^N->real) A = sum(1..dimindex(:N)) (\i. A$i$i)`;;
  20
  21 let TRACE_0 = prove
  22  (`trace(mat 0) = &0`,
  23   SIMP_TAC[trace; mat; LAMBDA_BETA; SUM_0]);;
  24
  25 let TRACE_I = prove
  26  (`trace(mat 1 :real^N^N) = &(dimindex(:N))`,
  27   SIMP_TAC[trace; mat; LAMBDA_BETA; SUM_CONST_NUMSEG; REAL_MUL_RID] THEN
  28   AP_TERM_TAC THEN ARITH_TAC);;
  29
  30 let TRACE_ADD = prove
  31  (`!A B:real^N^N. trace(A + B) = trace(A) + trace(B)`,
  32   SIMP_TAC[trace; matrix_add; SUM_ADD_NUMSEG; LAMBDA_BETA]);;
  33
  34 let TRACE_SUB = prove
  35  (`!A B:real^N^N. trace(A - B) = trace(A) - trace(B)`,
  36   SIMP_TAC[trace; matrix_sub; SUM_SUB_NUMSEG; LAMBDA_BETA]);;
  37
  38 let TRACE_MUL_SYM = prove
  39  (`!A B:real^N^N. trace(A ** B) = trace(B ** A)`,
  40   REPEAT GEN_TAC THEN SIMP_TAC[trace; matrix_mul; LAMBDA_BETA] THEN
  41   GEN_REWRITE_TAC RAND_CONV [SUM_SWAP_NUMSEG] THEN REWRITE_TAC[REAL_MUL_SYM]);;
  42
  43 let TRACE_TRANSP = prove
  44  (`!A:real^N^N. trace(transp A) = trace A`,
  45   SIMP_TAC[trace; transp; LAMBDA_BETA]);;
  46
  47 let TRACE_CONJUGATE = prove
  48  (`!A:real^N^N U:real^N^N.
  49         invertible U ==> trace(matrix_inv U ** A ** U) = trace A`,
  50   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[TRACE_MUL_SYM] THEN
  51   ASM_SIMP_TAC[GSYM MATRIX_MUL_ASSOC; MATRIX_INV; MATRIX_MUL_RID]);;
  52
  53 (* ------------------------------------------------------------------------- *)
  54 (* Definition of determinant.                                                *)
  55 (* ------------------------------------------------------------------------- *)
  56
  57 let det = new_definition
  58  `det(A:real^N^N) =
  59         sum { p | p permutes 1..dimindex(:N) }
  60             (\p. sign(p) * product (1..dimindex(:N)) (\i. A$i$(p i)))`;;
  61
  62 (* ------------------------------------------------------------------------- *)
  63 (* A few general lemmas we need below.                                       *)
  64 (* ------------------------------------------------------------------------- *)
  65
  66 let IN_DIMINDEX_SWAP = prove
  67  (`!m n j. 1 <= m /\ m <= dimindex(:N) /\
  68              1 <= n /\ n <= dimindex(:N) /\
  69              1 <= j /\ j <= dimindex(:N)
  70            ==> 1 <= swap(m,n) j /\ swap(m,n) j <= dimindex(:N)`,
  71   REWRITE_TAC[swap] THEN ARITH_TAC);;
  72
  73 let LAMBDA_BETA_PERM = prove
  74  (`!p i. p permutes 1..dimindex(:N) /\ 1 <= i /\ i <= dimindex(:N)
  75          ==> ((lambda) g :A^N) $ p(i) = g(p i)`,
  76   ASM_MESON_TAC[LAMBDA_BETA; PERMUTES_IN_IMAGE; IN_NUMSEG]);;
  77
  78 let PRODUCT_PERMUTE = prove
  79  (`!f p s. p permutes s ==> product s f = product s (f o p)`,
  80   REWRITE_TAC[product] THEN MATCH_MP_TAC ITERATE_PERMUTE THEN
  81   REWRITE_TAC[MONOIDAL_REAL_MUL]);;
  82
  83 let PRODUCT_PERMUTE_NUMSEG = prove
  84  (`!f p m n. p permutes m..n ==> product(m..n) f = product(m..n) (f o p)`,
  85   MESON_TAC[PRODUCT_PERMUTE; FINITE_NUMSEG]);;
  86
  87 let REAL_MUL_SUM = prove
  88  (`!s t f g.
  89         FINITE s /\ FINITE t
  90         ==> sum s f * sum t g = sum s (\i. sum t (\j. f(i) * g(j)))`,
  91   SIMP_TAC[SUM_LMUL] THEN
  92   ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN SIMP_TAC[SUM_LMUL]);;
  93
  94 let REAL_MUL_SUM_NUMSEG = prove
  95  (`!m n p q. sum(m..n) f * sum(p..q) g =
  96              sum(m..n) (\i. sum(p..q) (\j. f(i) * g(j)))`,
  97   SIMP_TAC[REAL_MUL_SUM; FINITE_NUMSEG]);;
  98
  99 (* ------------------------------------------------------------------------- *)
 100 (* Basic determinant properties.                                             *)
 101 (* ------------------------------------------------------------------------- *)
 102
 103 let DET_CMUL = prove
 104  (`!A:real^N^N c. det(c %% A) = c pow dimindex(:N) * det A`,
 105   REPEAT GEN_TAC THEN
 106   SIMP_TAC[det; MATRIX_CMUL_COMPONENT; PRODUCT_MUL; FINITE_NUMSEG] THEN
 107   SIMP_TAC[PRODUCT_CONST_NUMSEG_1; GSYM SUM_LMUL] THEN
 108   REWRITE_TAC[REAL_MUL_AC]);;
 109
 110 let DET_NEG = prove
 111  (`!A:real^N^N. det(--A) = --(&1) pow dimindex(:N) * det A`,
 112   REWRITE_TAC[MATRIX_NEG_MINUS1; DET_CMUL]);;
 113
 114 let DET_TRANSP = prove
 115  (`!A:real^N^N. det(transp A) = det A`,
 116   GEN_TAC THEN REWRITE_TAC[det] THEN
 117   GEN_REWRITE_TAC LAND_CONV [SUM_PERMUTATIONS_INVERSE] THEN
 118   MATCH_MP_TAC SUM_EQ THEN
 119   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN X_GEN_TAC `p:num->num` THEN
 120   REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN BINOP_TAC THENL
 121    [ASM_MESON_TAC[SIGN_INVERSE; PERMUTATION_PERMUTES; FINITE_NUMSEG];
 122     ALL_TAC] THEN
 123   FIRST_ASSUM(fun th -> GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV)
 124     [GSYM(MATCH_MP PERMUTES_IMAGE th)]) THEN
 125   MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC
 126    `product(1..dimindex(:N))
 127        ((\i. (transp A:real^N^N)$i$inverse p(i)) o p)` THEN
 128   CONJ_TAC THENL
 129    [MATCH_MP_TAC PRODUCT_IMAGE THEN
 130     ASM_MESON_TAC[FINITE_NUMSEG; PERMUTES_INJECTIVE; PERMUTES_INVERSE];
 131     MATCH_MP_TAC PRODUCT_EQ THEN REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG] THEN
 132     SIMP_TAC[transp; LAMBDA_BETA; o_THM] THEN
 133     FIRST_ASSUM(MP_TAC o MATCH_MP PERMUTES_INVERSES_o) THEN
 134     SIMP_TAC[FUN_EQ_THM; I_THM; o_THM] THEN STRIP_TAC THEN
 135     ASM_SIMP_TAC[PERMUTES_IN_NUMSEG; LAMBDA_BETA_PERM; LAMBDA_BETA]]);;
 136
 137 let DET_LOWERTRIANGULAR = prove
 138  (`!A:real^N^N.
 139         (!i j. 1 <= i /\ i <= dimindex(:N) /\
 140                1 <= j /\ j <= dimindex(:N) /\ i < j ==> A$i$j = &0)
 141         ==> det(A) = product(1..dimindex(:N)) (\i. A$i$i)`,
 142   REPEAT STRIP_TAC THEN REWRITE_TAC[det] THEN MATCH_MP_TAC EQ_TRANS THEN
 143   EXISTS_TAC `sum {I}
 144      (\p. sign p * product(1..dimindex(:N)) (\i. (A:real^N^N)$i$p(i)))` THEN
 145   CONJ_TAC THENL
 146    [ALL_TAC; REWRITE_TAC[SUM_SING; SIGN_I; REAL_MUL_LID; I_THM]] THEN
 147   MATCH_MP_TAC SUM_SUPERSET THEN
 148   SIMP_TAC[IN_SING; FINITE_RULES; SUBSET; IN_ELIM_THM; PERMUTES_I] THEN
 149   X_GEN_TAC `p:num->num` THEN STRIP_TAC THEN
 150   ASM_REWRITE_TAC[PRODUCT_EQ_0_NUMSEG; REAL_ENTIRE; SIGN_NZ] THEN
 151   MP_TAC(SPECL [`p:num->num`; `1..dimindex(:N)`] PERMUTES_NUMSET_LE) THEN
 152   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; NOT_LT]);;
 153
 154 let DET_UPPERTRIANGULAR = prove
 155  (`!A:real^N^N.
 156         (!i j. 1 <= i /\ i <= dimindex(:N) /\
 157                1 <= j /\ j <= dimindex(:N) /\ j < i ==> A$i$j = &0)
 158         ==> det(A) = product(1..dimindex(:N)) (\i. A$i$i)`,
 159   REPEAT STRIP_TAC THEN REWRITE_TAC[det] THEN MATCH_MP_TAC EQ_TRANS THEN
 160   EXISTS_TAC `sum {I}
 161      (\p. sign p * product(1..dimindex(:N)) (\i. (A:real^N^N)$i$p(i)))` THEN
 162   CONJ_TAC THENL
 163    [ALL_TAC; REWRITE_TAC[SUM_SING; SIGN_I; REAL_MUL_LID; I_THM]] THEN
 164   MATCH_MP_TAC SUM_SUPERSET THEN
 165   SIMP_TAC[IN_SING; FINITE_RULES; SUBSET; IN_ELIM_THM; PERMUTES_I] THEN
 166   X_GEN_TAC `p:num->num` THEN STRIP_TAC THEN
 167   ASM_REWRITE_TAC[PRODUCT_EQ_0_NUMSEG; REAL_ENTIRE; SIGN_NZ] THEN
 168   MP_TAC(SPECL [`p:num->num`; `1..dimindex(:N)`] PERMUTES_NUMSET_GE) THEN
 169   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; NOT_LT]);;
 170
 171 let DET_DIAGONAL = prove
 172  (`!A:real^N^N.
 173         (!i j. 1 <= i /\ i <= dimindex(:N) /\
 174                1 <= j /\ j <= dimindex(:N) /\ ~(i = j) ==> A$i$j = &0)
 175         ==> det(A) = product(1..dimindex(:N)) (\i. A$i$i)`,
 176   REPEAT STRIP_TAC THEN MATCH_MP_TAC DET_LOWERTRIANGULAR THEN
 177   REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[] THEN
 178   ASM_MESON_TAC[LT_REFL]);;
 179
 180 let DET_I = prove
 181  (`det(mat 1 :real^N^N) = &1`,
 182   MATCH_MP_TAC EQ_TRANS THEN
 183   EXISTS_TAC `product(1..dimindex(:N)) (\i. (mat 1:real^N^N)$i$i)` THEN
 184   CONJ_TAC THENL
 185    [MATCH_MP_TAC DET_LOWERTRIANGULAR;
 186     MATCH_MP_TAC PRODUCT_EQ_1_NUMSEG] THEN
 187   SIMP_TAC[mat; LAMBDA_BETA] THEN MESON_TAC[LT_REFL]);;
 188
 189 let DET_0 = prove
 190  (`det(mat 0 :real^N^N) = &0`,
 191   MATCH_MP_TAC EQ_TRANS THEN
 192   EXISTS_TAC `product(1..dimindex(:N)) (\i. (mat 0:real^N^N)$i$i)` THEN
 193   CONJ_TAC THENL
 194    [MATCH_MP_TAC DET_LOWERTRIANGULAR;
 195     REWRITE_TAC[PRODUCT_EQ_0_NUMSEG] THEN EXISTS_TAC `1`] THEN
 196   SIMP_TAC[mat; LAMBDA_BETA; COND_ID; DIMINDEX_GE_1; LE_REFL]);;
 197
 198 let DET_PERMUTE_ROWS = prove
 199  (`!A:real^N^N p.
 200         p permutes 1..dimindex(:N)
 201         ==> det(lambda i. A$p(i)) = sign(p) * det(A)`,
 202   REWRITE_TAC[det] THEN SIMP_TAC[LAMBDA_BETA] THEN REPEAT STRIP_TAC THEN
 203   SIMP_TAC[GSYM SUM_LMUL; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 204   FIRST_ASSUM(fun th -> GEN_REWRITE_TAC LAND_CONV
 205     [MATCH_MP SUM_PERMUTATIONS_COMPOSE_R th]) THEN
 206   MATCH_MP_TAC SUM_EQ THEN
 207   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN X_GEN_TAC `q:num->num` THEN
 208   REWRITE_TAC[IN_ELIM_THM; REAL_MUL_ASSOC] THEN DISCH_TAC THEN BINOP_TAC THENL
 209    [ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
 210     ASM_MESON_TAC[SIGN_COMPOSE; PERMUTATION_PERMUTES; FINITE_NUMSEG];
 211     ALL_TAC] THEN
 212   MP_TAC(MATCH_MP PERMUTES_INVERSE (ASSUME `p permutes 1..dimindex(:N)`)) THEN
 213   DISCH_THEN(fun th -> GEN_REWRITE_TAC LAND_CONV
 214     [MATCH_MP PRODUCT_PERMUTE_NUMSEG th]) THEN
 215   MATCH_MP_TAC PRODUCT_EQ THEN REWRITE_TAC[IN_NUMSEG; FINITE_NUMSEG] THEN
 216   REPEAT STRIP_TAC THEN REWRITE_TAC[o_THM] THEN
 217   ASM_MESON_TAC[PERMUTES_INVERSES]);;
 218
 219 let DET_PERMUTE_COLUMNS = prove
 220  (`!A:real^N^N p.
 221         p permutes 1..dimindex(:N)
 222         ==> det((lambda i j. A$i$p(j)):real^N^N) = sign(p) * det(A)`,
 223   REPEAT STRIP_TAC THEN
 224   GEN_REWRITE_TAC (funpow 2 RAND_CONV) [GSYM DET_TRANSP] THEN
 225   FIRST_ASSUM(fun th -> ONCE_REWRITE_TAC
 226    [GSYM(MATCH_MP DET_PERMUTE_ROWS th)]) THEN
 227   GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN AP_TERM_TAC THEN
 228   ASM_SIMP_TAC[CART_EQ; transp; LAMBDA_BETA; LAMBDA_BETA_PERM]);;
 229
 230 let DET_IDENTICAL_ROWS = prove
 231  (`!A:real^N^N i j. 1 <= i /\ i <= dimindex(:N) /\
 232                     1 <= j /\ j <= dimindex(:N) /\ ~(i = j) /\
 233                     row i A = row j A
 234                     ==> det A = &0`,
 235   REPEAT STRIP_TAC THEN
 236   MP_TAC(SPECL [`A:real^N^N`; `swap(i:num,j:num)`] DET_PERMUTE_ROWS) THEN
 237   ASM_SIMP_TAC[PERMUTES_SWAP; IN_NUMSEG; SIGN_SWAP] THEN
 238   MATCH_MP_TAC(REAL_ARITH `a = b ==> b = -- &1 * a ==> a = &0`) THEN
 239   AP_TERM_TAC THEN FIRST_X_ASSUM(MP_TAC o SYM) THEN
 240   SIMP_TAC[row; CART_EQ; LAMBDA_BETA] THEN
 241   REWRITE_TAC[swap] THEN ASM_MESON_TAC[]);;
 242
 243 let DET_IDENTICAL_COLUMNS = prove
 244  (`!A:real^N^N i j. 1 <= i /\ i <= dimindex(:N) /\
 245                     1 <= j /\ j <= dimindex(:N) /\ ~(i = j) /\
 246                     column i A = column j A
 247                     ==> det A = &0`,
 248   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[GSYM DET_TRANSP] THEN
 249   MATCH_MP_TAC DET_IDENTICAL_ROWS THEN ASM_MESON_TAC[ROW_TRANSP]);;
 250
 251 let DET_ZERO_ROW = prove
 252  (`!A:real^N^N i.
 253        1 <= i /\ i <= dimindex(:N) /\ row i A = vec 0  ==> det A = &0`,
 254   SIMP_TAC[det; row; CART_EQ; LAMBDA_BETA; VEC_COMPONENT] THEN
 255   REPEAT STRIP_TAC THEN MATCH_MP_TAC SUM_EQ_0 THEN
 256   REWRITE_TAC[IN_ELIM_THM; REAL_ENTIRE; SIGN_NZ] THEN REPEAT STRIP_TAC THEN
 257   SIMP_TAC[PRODUCT_EQ_0; FINITE_NUMSEG; IN_NUMSEG] THEN
 258   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]);;
 259
 260 let DET_ZERO_COLUMN = prove
 261  (`!A:real^N^N i.
 262        1 <= i /\ i <= dimindex(:N) /\ column i A = vec 0  ==> det A = &0`,
 263   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[GSYM DET_TRANSP] THEN
 264   MATCH_MP_TAC DET_ZERO_ROW THEN ASM_MESON_TAC[ROW_TRANSP]);;
 265
 266 let DET_ROW_ADD = prove
 267  (`!a b c k.
 268          1 <= k /\ k <= dimindex(:N)
 269          ==> det((lambda i. if i = k then a + b else c i):real^N^N) =
 270              det((lambda i. if i = k then a else c i):real^N^N) +
 271              det((lambda i. if i = k then b else c i):real^N^N)`,
 272   SIMP_TAC[det; LAMBDA_BETA; GSYM SUM_ADD;
 273            FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 274   REPEAT STRIP_TAC THEN MATCH_MP_TAC SUM_EQ THEN
 275   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 276   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN
 277   DISCH_TAC THEN REWRITE_TAC[GSYM REAL_ADD_LDISTRIB] THEN AP_TERM_TAC THEN
 278   SUBGOAL_THEN `1..dimindex(:N) = k INSERT ((1..dimindex(:N)) DELETE k)`
 279   SUBST1_TAC THENL [ASM_MESON_TAC[INSERT_DELETE; IN_NUMSEG]; ALL_TAC] THEN
 280   SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE] THEN
 281   MATCH_MP_TAC(REAL_RING
 282    `c = a + b /\ y = x:real /\ z = x ==> c * x = a * y + b * z`) THEN
 283   REWRITE_TAC[VECTOR_ADD_COMPONENT] THEN
 284   CONJ_TAC THEN MATCH_MP_TAC PRODUCT_EQ THEN
 285   SIMP_TAC[IN_DELETE; FINITE_DELETE; FINITE_NUMSEG]);;
 286
 287 let DET_ROW_MUL = prove
 288  (`!a b c k.
 289         1 <= k /\ k <= dimindex(:N)
 290         ==> det((lambda i. if i = k then c % a else b i):real^N^N) =
 291             c * det((lambda i. if i = k then a else b i):real^N^N)`,
 292   SIMP_TAC[det; LAMBDA_BETA; GSYM SUM_LMUL;
 293            FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 294   REPEAT STRIP_TAC THEN MATCH_MP_TAC SUM_EQ THEN
 295   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 296   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 297   SUBGOAL_THEN `1..dimindex(:N) = k INSERT ((1..dimindex(:N)) DELETE k)`
 298   SUBST1_TAC THENL [ASM_MESON_TAC[INSERT_DELETE; IN_NUMSEG]; ALL_TAC] THEN
 299   SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE] THEN
 300   MATCH_MP_TAC(REAL_RING
 301    `cp = c * p /\ p1 = p2:real ==> s * cp * p1 = c * s * p * p2`) THEN
 302   REWRITE_TAC[VECTOR_MUL_COMPONENT] THEN MATCH_MP_TAC PRODUCT_EQ THEN
 303   SIMP_TAC[IN_DELETE; FINITE_DELETE; FINITE_NUMSEG]);;
 304
 305 let DET_ROW_OPERATION = prove
 306  (`!A:real^N^N i.
 307         1 <= i /\ i <= dimindex(:N) /\
 308         1 <= j /\ j <= dimindex(:N) /\ ~(i = j)
 309         ==> det(lambda k. if k = i then row i A + c % row j A else row k A) =
 310             det A`,
 311   REPEAT STRIP_TAC THEN ASM_SIMP_TAC[DET_ROW_ADD; DET_ROW_MUL] THEN
 312   MATCH_MP_TAC(REAL_RING `a = b /\ d = &0 ==> a + c * d = b`) THEN
 313   CONJ_TAC THENL
 314    [AP_TERM_TAC THEN ASM_SIMP_TAC[LAMBDA_BETA; CART_EQ] THEN
 315     REPEAT STRIP_TAC THEN COND_CASES_TAC THEN
 316     ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ];
 317     MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 318     MAP_EVERY EXISTS_TAC [`i:num`; `j:num`] THEN
 319     ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ]]);;
 320
 321 let DET_ROW_SPAN = prove
 322  (`!A:real^N^N i x.
 323         1 <= i /\ i <= dimindex(:N) /\
 324         x IN span {row j A | 1 <= j /\ j <= dimindex(:N) /\ ~(j = i)}
 325         ==> det(lambda k. if k = i then row i A + x else row k A) =
 326             det A`,
 327   GEN_TAC THEN GEN_TAC THEN
 328   REWRITE_TAC[IMP_CONJ; RIGHT_FORALL_IMP_THM] THEN REPEAT DISCH_TAC THEN
 329   MATCH_MP_TAC SPAN_INDUCT_ALT THEN CONJ_TAC THENL
 330    [AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA; VECTOR_ADD_RID] THEN
 331     REPEAT STRIP_TAC THEN COND_CASES_TAC THEN ASM_SIMP_TAC[row; LAMBDA_BETA];
 332     ALL_TAC] THEN
 333   REPEAT GEN_TAC THEN REWRITE_TAC[IN_ELIM_THM] THEN
 334   DISCH_THEN(CONJUNCTS_THEN2 (X_CHOOSE_TAC `j:num`) (SUBST_ALL_TAC o SYM)) THEN
 335   ONCE_REWRITE_TAC[VECTOR_ARITH
 336      `a + c % x + y:real^N = (a + y) + c % x`] THEN
 337   ABBREV_TAC `z = row i (A:real^N^N) + y` THEN
 338   ASM_SIMP_TAC[DET_ROW_MUL; DET_ROW_ADD] THEN
 339   MATCH_MP_TAC(REAL_RING `d = &0 ==> a + c * d = a`) THEN
 340   MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 341   MAP_EVERY EXISTS_TAC [`i:num`; `j:num`] THEN
 342   ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ]);;
 343
 344 (* ------------------------------------------------------------------------- *)
 345 (* May as well do this, though it's a bit unsatisfactory since it ignores    *)
 346 (* exact duplicates by considering the rows/columns as a set.                *)
 347 (* ------------------------------------------------------------------------- *)
 348
 349 let DET_DEPENDENT_ROWS = prove
 350  (`!A:real^N^N. dependent(rows A) ==> det A = &0`,
 351   GEN_TAC THEN
 352   REWRITE_TAC[dependent; rows; IN_ELIM_THM; LEFT_AND_EXISTS_THM] THEN
 353   REWRITE_TAC[LEFT_IMP_EXISTS_THM] THEN GEN_TAC THEN X_GEN_TAC `i:num` THEN
 354   STRIP_TAC THEN FIRST_X_ASSUM SUBST_ALL_TAC THEN
 355   ASM_CASES_TAC
 356    `?i j. 1 <= i /\ i <= dimindex(:N) /\
 357           1 <= j /\ j <= dimindex(:N) /\ ~(i = j) /\
 358           row i (A:real^N^N) = row j A`
 359   THENL [ASM_MESON_TAC[DET_IDENTICAL_ROWS]; ALL_TAC] THEN
 360   MP_TAC(SPECL [`A:real^N^N`; `i:num`; `--(row i (A:real^N^N))`]
 361     DET_ROW_SPAN) THEN
 362   ANTS_TAC THENL
 363    [ASM_REWRITE_TAC[] THEN MATCH_MP_TAC SPAN_NEG THEN
 364     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [IN]) THEN
 365     MATCH_MP_TAC(TAUT `a = b ==> a ==> b`) THEN
 366     REWRITE_TAC[IN] THEN AP_THM_TAC THEN AP_TERM_TAC THEN
 367     REWRITE_TAC[EXTENSION; IN_DELETE; IN_ELIM_THM] THEN ASM_MESON_TAC[];
 368     DISCH_THEN(SUBST1_TAC o SYM) THEN MATCH_MP_TAC DET_ZERO_ROW THEN
 369     EXISTS_TAC `i:num` THEN
 370     ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ; VECTOR_ADD_COMPONENT;
 371                  VECTOR_NEG_COMPONENT; VEC_COMPONENT] THEN
 372     REAL_ARITH_TAC]);;
 373
 374 let DET_DEPENDENT_COLUMNS = prove
 375  (`!A:real^N^N. dependent(columns A) ==> det A = &0`,
 376   MESON_TAC[DET_DEPENDENT_ROWS; ROWS_TRANSP; DET_TRANSP]);;
 377
 378 (* ------------------------------------------------------------------------- *)
 379 (* Multilinearity and the multiplication formula.                            *)
 380 (* ------------------------------------------------------------------------- *)
 381
 382 let DET_LINEAR_ROW_VSUM = prove
 383  (`!a c s k.
 384          FINITE s /\ 1 <= k /\ k <= dimindex(:N)
 385          ==> det((lambda i. if i = k then vsum s a else c i):real^N^N) =
 386              sum s
 387                (\j. det((lambda i. if i = k then a(j) else c i):real^N^N))`,
 388   GEN_TAC THEN GEN_TAC THEN ONCE_REWRITE_TAC[IMP_CONJ] THEN
 389   REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN
 390   MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
 391   SIMP_TAC[VSUM_CLAUSES; SUM_CLAUSES; DET_ROW_ADD] THEN
 392   REPEAT STRIP_TAC THEN MATCH_MP_TAC DET_ZERO_ROW THEN EXISTS_TAC `k:num` THEN
 393   ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ; VEC_COMPONENT]);;
 394
 395 let BOUNDED_FUNCTIONS_BIJECTIONS_1 = prove
 396  (`!p. p IN {(y,g) | y IN s /\
 397                      g IN {f | (!i. 1 <= i /\ i <= k ==> f i IN s) /\
 398                                (!i. ~(1 <= i /\ i <= k) ==> f i = i)}}
 399        ==> (\(y,g) i. if i = SUC k then y else g(i)) p IN
 400              {f | (!i. 1 <= i /\ i <= SUC k ==> f i IN s) /\
 401                   (!i. ~(1 <= i /\ i <= SUC k) ==> f i = i)} /\
 402            (\h. h(SUC k),(\i. if i = SUC k then i else h(i)))
 403             ((\(y,g) i. if i = SUC k then y else g(i)) p) = p`,
 404   REWRITE_TAC[FORALL_PAIR_THM; IN_ELIM_PAIR_THM] THEN
 405   CONV_TAC(REDEPTH_CONV GEN_BETA_CONV) THEN REWRITE_TAC[IN_ELIM_THM] THEN
 406   MAP_EVERY X_GEN_TAC [`y:num`; `h:num->num`] THEN REPEAT STRIP_TAC THENL
 407    [ASM_MESON_TAC[LE];
 408     ASM_MESON_TAC[LE; ARITH_RULE `~(1 <= i /\ i <= SUC k) ==> ~(i = SUC k)`];
 409     REWRITE_TAC[PAIR_EQ; FUN_EQ_THM] THEN
 410     ASM_MESON_TAC[ARITH_RULE `~(SUC k <= k)`]]);;
 411
 412 let BOUNDED_FUNCTIONS_BIJECTIONS_2 = prove
 413  (`!h. h IN {f | (!i. 1 <= i /\ i <= SUC k ==> f i IN s) /\
 414                  (!i. ~(1 <= i /\ i <= SUC k) ==> f i = i)}
 415        ==> (\h. h(SUC k),(\i. if i = SUC k then i else h(i))) h IN
 416            {(y,g) | y IN s /\
 417                      g IN {f | (!i. 1 <= i /\ i <= k ==> f i IN s) /\
 418                                (!i. ~(1 <= i /\ i <= k) ==> f i = i)}} /\
 419            (\(y,g) i. if i = SUC k then y else g(i))
 420               ((\h. h(SUC k),(\i. if i = SUC k then i else h(i))) h) = h`,
 421   REWRITE_TAC[IN_ELIM_PAIR_THM] THEN
 422   CONV_TAC(REDEPTH_CONV GEN_BETA_CONV) THEN REWRITE_TAC[IN_ELIM_THM] THEN
 423   X_GEN_TAC `h:num->num` THEN REPEAT STRIP_TAC THENL
 424    [FIRST_X_ASSUM MATCH_MP_TAC THEN ARITH_TAC;
 425     ASM_MESON_TAC[ARITH_RULE `i <= k ==> i <= SUC k /\ ~(i = SUC k)`];
 426     ASM_MESON_TAC[ARITH_RULE `i <= SUC k /\ ~(i = SUC k) ==> i <= k`];
 427     REWRITE_TAC[FUN_EQ_THM] THEN ASM_MESON_TAC[LE_REFL]]);;
 428
 429 let FINITE_BOUNDED_FUNCTIONS = prove
 430  (`!s k. FINITE s
 431          ==> FINITE {f | (!i. 1 <= i /\ i <= k ==> f(i) IN s) /\
 432                          (!i. ~(1 <= i /\ i <= k) ==> f(i) = i)}`,
 433   REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN GEN_TAC THEN DISCH_TAC THEN
 434   INDUCT_TAC THENL
 435    [REWRITE_TAC[ARITH_RULE `~(1 <= i /\ i <= 0)`] THEN
 436     SIMP_TAC[GSYM FUN_EQ_THM; SET_RULE `{x | x = y} = {y}`; FINITE_RULES];
 437     ALL_TAC] THEN
 438   UNDISCH_TAC `FINITE(s:num->bool)` THEN POP_ASSUM MP_TAC THEN
 439   REWRITE_TAC[TAUT `a ==> b ==> c <=> b /\ a ==> c`] THEN
 440   DISCH_THEN(MP_TAC o MATCH_MP FINITE_PRODUCT) THEN
 441   DISCH_THEN(MP_TAC o ISPEC `\(y:num,g) i. if i = SUC k then y else g(i)` o
 442                       MATCH_MP FINITE_IMAGE) THEN
 443   MATCH_MP_TAC(TAUT `a = b ==> a ==> b`) THEN AP_TERM_TAC THEN
 444   REWRITE_TAC[EXTENSION; IN_IMAGE] THEN
 445   X_GEN_TAC `h:num->num` THEN EQ_TAC THENL
 446    [STRIP_TAC THEN ASM_SIMP_TAC[BOUNDED_FUNCTIONS_BIJECTIONS_1]; ALL_TAC] THEN
 447   DISCH_TAC THEN EXISTS_TAC
 448     `(\h. h(SUC k),(\i. if i = SUC k then i else h(i))) h` THEN
 449   PURE_ONCE_REWRITE_TAC[CONJ_SYM] THEN CONV_TAC (RAND_CONV SYM_CONV) THEN
 450   MATCH_MP_TAC BOUNDED_FUNCTIONS_BIJECTIONS_2 THEN ASM_REWRITE_TAC[]);;
 451
 452 let DET_LINEAR_ROWS_VSUM_LEMMA = prove
 453  (`!s k a c.
 454          FINITE s /\ k <= dimindex(:N)
 455          ==> det((lambda i. if i <= k then vsum s (a i) else c i):real^N^N) =
 456              sum {f | (!i. 1 <= i /\ i <= k ==> f(i) IN s) /\
 457                       !i. ~(1 <= i /\ i <= k) ==> f(i) = i}
 458                  (\f. det((lambda i. if i <= k then a i (f i) else c i)
 459                           :real^N^N))`,
 460   let lemma = prove
 461    (`(lambda i. if i <= 0 then x(i) else y(i)) = (lambda i. y i)`,
 462     SIMP_TAC[CART_EQ; ARITH; LAMBDA_BETA; ARITH_RULE
 463                  `1 <= k ==> ~(k <= 0)`]) in
 464   ONCE_REWRITE_TAC[IMP_CONJ] THEN
 465   REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN GEN_TAC THEN DISCH_TAC THEN
 466   INDUCT_TAC THENL
 467    [REWRITE_TAC[lemma; LE_0] THEN GEN_TAC THEN
 468     REWRITE_TAC[ARITH_RULE `~(1 <= i /\ i <= 0)`] THEN
 469     REWRITE_TAC[GSYM FUN_EQ_THM; SET_RULE `{x | x = y} = {y}`] THEN
 470     REWRITE_TAC[SUM_SING];
 471     ALL_TAC] THEN
 472   DISCH_TAC THEN FIRST_X_ASSUM(MP_TAC o check (is_imp o concl)) THEN
 473   ASM_SIMP_TAC[ARITH_RULE `SUC k <= n ==> k <= n`] THEN REPEAT STRIP_TAC THEN
 474   GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [LE] THEN
 475   REWRITE_TAC[TAUT
 476    `(if a \/ b then c else d) = (if a then c else if b then c else d)`] THEN
 477   ASM_SIMP_TAC[DET_LINEAR_ROW_VSUM; ARITH_RULE `1 <= SUC k`] THEN
 478   ONCE_REWRITE_TAC[TAUT
 479     `(if a then b else if c then d else e) =
 480      (if c then (if a then b else d) else (if a then b else e))`] THEN
 481   ASM_SIMP_TAC[ARITH_RULE `i <= k ==> ~(i = SUC k)`] THEN
 482   ASM_SIMP_TAC[SUM_SUM_PRODUCT; FINITE_BOUNDED_FUNCTIONS] THEN
 483   MATCH_MP_TAC SUM_EQ_GENERAL_INVERSES THEN
 484   EXISTS_TAC `\(y:num,g) i. if i = SUC k then y else g(i)` THEN
 485   EXISTS_TAC `\h. h(SUC k),(\i. if i = SUC k then i else h(i))` THEN
 486   CONJ_TAC THENL [ACCEPT_TAC BOUNDED_FUNCTIONS_BIJECTIONS_2; ALL_TAC] THEN
 487   X_GEN_TAC `p:num#(num->num)` THEN
 488   DISCH_THEN(STRIP_ASSUME_TAC o MATCH_MP BOUNDED_FUNCTIONS_BIJECTIONS_1) THEN
 489   ASM_REWRITE_TAC[] THEN
 490   SPEC_TAC(`p:num#(num->num)`,`q:num#(num->num)`) THEN
 491   REWRITE_TAC[FORALL_PAIR_THM] THEN
 492   CONV_TAC(ONCE_DEPTH_CONV GEN_BETA_CONV) THEN
 493   MAP_EVERY X_GEN_TAC [`y:num`; `g:num->num`] THEN AP_TERM_TAC THEN
 494   SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 495   REPEAT STRIP_TAC THEN REPEAT(COND_CASES_TAC THEN ASM_REWRITE_TAC[]) THEN
 496   ASM_MESON_TAC[LE; ARITH_RULE `~(SUC k <= k)`]);;
 497
 498 let DET_LINEAR_ROWS_VSUM = prove
 499  (`!s a.
 500          FINITE s
 501          ==> det((lambda i. vsum s (a i)):real^N^N) =
 502              sum {f | (!i. 1 <= i /\ i <= dimindex(:N) ==> f(i) IN s) /\
 503                       !i. ~(1 <= i /\ i <= dimindex(:N)) ==> f(i) = i}
 504                  (\f. det((lambda i. a i (f i)):real^N^N))`,
 505   let lemma = prove
 506    (`(lambda i. if i <= dimindex(:N) then x(i) else y(i)):real^N^N =
 507      (lambda i. x(i))`,
 508     SIMP_TAC[CART_EQ; LAMBDA_BETA]) in
 509   REPEAT STRIP_TAC THEN
 510   MP_TAC(SPECL [`s:num->bool`; `dimindex(:N)`] DET_LINEAR_ROWS_VSUM_LEMMA) THEN
 511   ASM_REWRITE_TAC[LE_REFL; lemma] THEN SIMP_TAC[]);;
 512
 513 let MATRIX_MUL_VSUM_ALT = prove
 514  (`!A:real^N^N B:real^N^N. A ** B =
 515                   lambda i. vsum (1..dimindex(:N)) (\k. A$i$k % B$k)`,
 516   SIMP_TAC[matrix_mul; CART_EQ; LAMBDA_BETA; VECTOR_MUL_COMPONENT;
 517            VSUM_COMPONENT]);;
 518
 519 let DET_ROWS_MUL = prove
 520  (`!a c. det((lambda i. c(i) % a(i)):real^N^N) =
 521          product(1..dimindex(:N)) (\i. c(i)) *
 522          det((lambda i. a(i)):real^N^N)`,
 523   REPEAT GEN_TAC THEN SIMP_TAC[det; LAMBDA_BETA] THEN
 524   SIMP_TAC[GSYM SUM_LMUL; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 525   MATCH_MP_TAC SUM_EQ THEN SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 526   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 527   MATCH_MP_TAC(REAL_RING `b = c * d ==> s * b = c * s * d`) THEN
 528   SIMP_TAC[GSYM PRODUCT_MUL_NUMSEG] THEN
 529   MATCH_MP_TAC PRODUCT_EQ_NUMSEG THEN
 530   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; VECTOR_MUL_COMPONENT]);;
 531
 532 let DET_MUL = prove
 533  (`!A B:real^N^N. det(A ** B) = det(A) * det(B)`,
 534   REPEAT GEN_TAC THEN REWRITE_TAC[MATRIX_MUL_VSUM_ALT] THEN
 535   SIMP_TAC[DET_LINEAR_ROWS_VSUM; FINITE_NUMSEG] THEN
 536   MATCH_MP_TAC EQ_TRANS THEN
 537   EXISTS_TAC `sum {p | p permutes 1..dimindex(:N)}
 538             (\f. det (lambda i. (A:real^N^N)$i$f i % (B:real^N^N)$f i))` THEN
 539   CONJ_TAC THENL
 540    [REWRITE_TAC[DET_ROWS_MUL] THEN
 541     MATCH_MP_TAC SUM_SUPERSET THEN
 542     SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 543     REWRITE_TAC[SUBSET; IN_ELIM_THM] THEN CONJ_TAC THENL
 544      [MESON_TAC[permutes; IN_NUMSEG]; ALL_TAC] THEN
 545     X_GEN_TAC `f:num->num` THEN REWRITE_TAC[permutes; IN_NUMSEG] THEN
 546     DISCH_THEN(CONJUNCTS_THEN2 STRIP_ASSUME_TAC MP_TAC) THEN
 547     ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
 548     REWRITE_TAC[REAL_ENTIRE] THEN DISJ2_TAC THEN
 549     MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 550     MP_TAC(ISPECL [`1..dimindex(:N)`; `f:num->num`]
 551        SURJECTIVE_IFF_INJECTIVE) THEN
 552     ASM_REWRITE_TAC[SUBSET; IN_NUMSEG; FINITE_NUMSEG; FORALL_IN_IMAGE] THEN
 553     MATCH_MP_TAC(TAUT `(~b ==> c) /\ (b ==> ~a) ==> (a <=> b) ==> c`) THEN
 554     CONJ_TAC THENL
 555      [REWRITE_TAC[NOT_FORALL_THM] THEN
 556       REPEAT(MATCH_MP_TAC MONO_EXISTS THEN GEN_TAC) THEN
 557       SIMP_TAC[CART_EQ; LAMBDA_BETA; row; NOT_IMP];
 558       ALL_TAC] THEN
 559     DISCH_TAC THEN
 560     SUBGOAL_THEN `!x y. (f:num->num)(x) = f(y) ==> x = y` ASSUME_TAC THENL
 561      [REPEAT GEN_TAC THEN
 562       ASM_CASES_TAC `1 <= x /\ x <= dimindex(:N)` THEN
 563       ASM_CASES_TAC `1 <= y /\ y <= dimindex(:N)` THEN
 564       ASM_MESON_TAC[];
 565       ALL_TAC] THEN
 566     ASM_MESON_TAC[];
 567     ALL_TAC] THEN
 568   SIMP_TAC[det; REAL_MUL_SUM; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 569   MATCH_MP_TAC SUM_EQ THEN SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 570   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 571   FIRST_ASSUM(fun th -> GEN_REWRITE_TAC RAND_CONV
 572     [MATCH_MP SUM_PERMUTATIONS_COMPOSE_R (MATCH_MP PERMUTES_INVERSE th)]) THEN
 573   MATCH_MP_TAC SUM_EQ THEN SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 574   X_GEN_TAC `q:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 575   REWRITE_TAC[o_THM] THEN ONCE_REWRITE_TAC[AC REAL_MUL_AC
 576    `(p * x) * (q * y) = (p * q) * (x * y)`] THEN
 577   BINOP_TAC THENL
 578    [SUBGOAL_THEN `sign(q o inverse p) = sign(p:num->num) * sign(q:num->num)`
 579      (fun t -> SIMP_TAC[REAL_MUL_ASSOC; SIGN_IDEMPOTENT; REAL_MUL_LID; t]) THEN
 580     ASM_MESON_TAC[SIGN_COMPOSE; PERMUTES_INVERSE; PERMUTATION_PERMUTES;
 581                   FINITE_NUMSEG; SIGN_INVERSE; REAL_MUL_SYM];
 582     ALL_TAC] THEN
 583   GEN_REWRITE_TAC (RAND_CONV o RAND_CONV)
 584    [MATCH_MP PRODUCT_PERMUTE_NUMSEG (ASSUME `p permutes 1..dimindex(:N)`)] THEN
 585   SIMP_TAC[GSYM PRODUCT_MUL; FINITE_NUMSEG] THEN
 586   MATCH_MP_TAC PRODUCT_EQ_NUMSEG THEN
 587   ASM_SIMP_TAC[LAMBDA_BETA; LAMBDA_BETA_PERM; o_THM] THEN
 588   X_GEN_TAC `i:num` THEN STRIP_TAC THEN MATCH_MP_TAC EQ_TRANS THEN
 589   EXISTS_TAC `(A:real^N^N)$i$p(i) * (B:real^N^N)$p(i)$q(i)` THEN CONJ_TAC THENL
 590    [ASM_MESON_TAC[VECTOR_MUL_COMPONENT; PERMUTES_IN_IMAGE; IN_NUMSEG];
 591     ASM_MESON_TAC[PERMUTES_INVERSES]]);;
 592
 593 (* ------------------------------------------------------------------------- *)
 594 (* Relation to invertibility.                                                *)
 595 (* ------------------------------------------------------------------------- *)
 596
 597 let INVERTIBLE_DET_NZ = prove
 598  (`!A:real^N^N. invertible(A) <=> ~(det A = &0)`,
 599   GEN_TAC THEN EQ_TAC THENL
 600    [REWRITE_TAC[INVERTIBLE_RIGHT_INVERSE; LEFT_IMP_EXISTS_THM] THEN
 601     GEN_TAC THEN DISCH_THEN(MP_TAC o AP_TERM `det:real^N^N->real`) THEN
 602     REWRITE_TAC[DET_MUL; DET_I] THEN CONV_TAC REAL_RING;
 603     ALL_TAC] THEN
 604   ONCE_REWRITE_TAC[GSYM CONTRAPOS_THM] THEN
 605   REWRITE_TAC[INVERTIBLE_RIGHT_INVERSE] THEN
 606   REWRITE_TAC[MATRIX_RIGHT_INVERTIBLE_INDEPENDENT_ROWS] THEN
 607   REWRITE_TAC[NOT_FORALL_THM; NOT_IMP] THEN
 608   REWRITE_TAC[RIGHT_AND_EXISTS_THM; LEFT_IMP_EXISTS_THM] THEN
 609   MAP_EVERY X_GEN_TAC [`c:num->real`; `i:num`] THEN STRIP_TAC THEN
 610   MP_TAC(SPECL [`A:real^N^N`; `i:num`; `--(row i (A:real^N^N))`]
 611     DET_ROW_SPAN) THEN
 612   ANTS_TAC THENL
 613    [ASM_REWRITE_TAC[] THEN
 614     SUBGOAL_THEN
 615       `--(row i (A:real^N^N)) =
 616        vsum ((1..dimindex(:N)) DELETE i) (\j. inv(c i) % c j % row j A)`
 617     SUBST1_TAC THENL
 618      [ASM_SIMP_TAC[VSUM_DELETE_CASES; FINITE_NUMSEG; IN_NUMSEG; VSUM_LMUL] THEN
 619       ASM_SIMP_TAC[VECTOR_MUL_ASSOC; REAL_MUL_LINV] THEN VECTOR_ARITH_TAC;
 620       ALL_TAC] THEN
 621     MATCH_MP_TAC SPAN_VSUM THEN
 622     REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG; FINITE_DELETE; IN_DELETE] THEN
 623     X_GEN_TAC `j:num` THEN STRIP_TAC THEN REPEAT(MATCH_MP_TAC SPAN_MUL) THEN
 624     MATCH_MP_TAC(CONJUNCT1 SPAN_CLAUSES) THEN
 625     REWRITE_TAC[IN_ELIM_THM] THEN ASM_MESON_TAC[];
 626     ALL_TAC] THEN
 627   DISCH_THEN(SUBST1_TAC o SYM) THEN MATCH_MP_TAC DET_ZERO_ROW THEN
 628   EXISTS_TAC `i:num` THEN
 629   ASM_SIMP_TAC[row; CART_EQ; LAMBDA_BETA; VEC_COMPONENT;
 630                VECTOR_ARITH `x + --x:real^N = vec 0`]);;
 631
 632 let DET_EQ_0 = prove
 633  (`!A:real^N^N. det(A) = &0 <=> ~invertible(A)`,
 634   REWRITE_TAC[INVERTIBLE_DET_NZ]);;
 635
 636 let MATRIX_MUL_LINV = prove
 637  (`!A:real^N^N. ~(det A = &0) ==> matrix_inv A ** A = mat 1`,
 638   SIMP_TAC[MATRIX_INV; DET_EQ_0]);;
 639
 640 let MATRIX_MUL_RINV = prove
 641  (`!A:real^N^N. ~(det A = &0) ==> A ** matrix_inv A = mat 1`,
 642   SIMP_TAC[MATRIX_INV; DET_EQ_0]);;
 643
 644 let DET_MATRIX_EQ_0 = prove
 645  (`!f:real^N->real^N.
 646         linear f
 647         ==> (det(matrix f) = &0 <=>
 648              ~(?g. linear g /\ f o g = I /\ g o f = I))`,
 649   SIMP_TAC[DET_EQ_0; MATRIX_INVERTIBLE]);;
 650
 651 let DET_MATRIX_EQ_0_LEFT = prove
 652  (`!f:real^N->real^N.
 653         linear f
 654         ==> (det(matrix f) = &0 <=>
 655              ~(?g. linear g /\ g o f = I))`,
 656    SIMP_TAC[DET_MATRIX_EQ_0] THEN MESON_TAC[LINEAR_INVERSE_LEFT]);;
 657
 658 let DET_MATRIX_EQ_0_RIGHT = prove
 659  (`!f:real^N->real^N.
 660         linear f
 661         ==> (det(matrix f) = &0 <=>
 662              ~(?g. linear g /\ f o g = I))`,
 663    SIMP_TAC[DET_MATRIX_EQ_0] THEN MESON_TAC[LINEAR_INVERSE_LEFT]);;
 664
 665 let DET_EQ_0_RANK = prove
 666  (`!A:real^N^N. det A = &0 <=> rank A < dimindex(:N)`,
 667   REWRITE_TAC[DET_EQ_0; INVERTIBLE_LEFT_INVERSE; GSYM FULL_RANK_INJECTIVE;
 668               MATRIX_LEFT_INVERTIBLE_INJECTIVE] THEN
 669   GEN_TAC THEN MP_TAC(ISPEC `A:real^N^N` RANK_BOUND) THEN
 670   ARITH_TAC);;
 671
 672 let HOMOGENEOUS_LINEAR_EQUATIONS_DET = prove
 673  (`!A:real^N^N. (?x. ~(x = vec 0) /\ A ** x = vec 0) <=> det A = &0`,
 674   GEN_TAC THEN
 675   REWRITE_TAC[MATRIX_NONFULL_LINEAR_EQUATIONS_EQ; DET_EQ_0_RANK] THEN
 676   MATCH_MP_TAC(ARITH_RULE `r <= MIN N N ==> (~(r = N) <=> r < N)`) THEN
 677   REWRITE_TAC[RANK_BOUND]);;
 678
 679 let INVERTIBLE_MATRIX_MUL = prove
 680  (`!A:real^N^N B:real^N^N.
 681         invertible(A ** B) <=> invertible A /\ invertible B`,
 682   REWRITE_TAC[INVERTIBLE_DET_NZ; DET_MUL; DE_MORGAN_THM; REAL_ENTIRE]);;
 683
 684 let MATRIX_INV_MUL = prove
 685  (`!A:real^N^N B:real^N^N.
 686         invertible A /\ invertible B
 687         ==> matrix_inv(A ** B) = matrix_inv B ** matrix_inv A`,
 688   REPEAT STRIP_TAC THEN MATCH_MP_TAC MATRIX_INV_UNIQUE THEN
 689   ONCE_REWRITE_TAC[MATRIX_MUL_ASSOC] THEN
 690   GEN_REWRITE_TAC (BINOP_CONV o LAND_CONV o LAND_CONV)
 691    [GSYM MATRIX_MUL_ASSOC] THEN
 692   ASM_SIMP_TAC[MATRIX_MUL_LINV; DET_EQ_0; MATRIX_MUL_RID; MATRIX_MUL_RINV]);;
 693
 694 (* ------------------------------------------------------------------------- *)
 695 (* Cramer's rule.                                                            *)
 696 (* ------------------------------------------------------------------------- *)
 697
 698 let CRAMER_LEMMA_TRANSP = prove
 699  (`!A:real^N^N x:real^N.
 700         1 <= k /\ k <= dimindex(:N)
 701         ==> det((lambda i. if i = k
 702                            then vsum(1..dimindex(:N)) (\i. x$i % row i A)
 703                            else row i A):real^N^N) =
 704             x$k * det A`,
 705   REPEAT STRIP_TAC THEN
 706   SUBGOAL_THEN `1..dimindex(:N) = k INSERT ((1..dimindex(:N)) DELETE k)`
 707   SUBST1_TAC THENL [ASM_MESON_TAC[INSERT_DELETE; IN_NUMSEG]; ALL_TAC] THEN
 708   SIMP_TAC[VSUM_CLAUSES; FINITE_NUMSEG; FINITE_DELETE; IN_DELETE] THEN
 709   REWRITE_TAC[VECTOR_ARITH
 710    `(x:real^N)$k % row k (A:real^N^N) + s =
 711     (x$k - &1) % row k A + row k A + s`] THEN
 712   W(MP_TAC o PART_MATCH (lhs o rand) DET_ROW_ADD o lhand o snd) THEN
 713   ASM_SIMP_TAC[DET_ROW_MUL] THEN DISCH_THEN(K ALL_TAC) THEN
 714   MATCH_MP_TAC(REAL_RING `d = d' /\ e = d' ==> (c - &1) * d + e = c * d'`) THEN
 715   CONJ_TAC THENL
 716    [AP_TERM_TAC THEN ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 717     REPEAT STRIP_TAC THEN COND_CASES_TAC THEN ASM_SIMP_TAC[LAMBDA_BETA; row];
 718     MATCH_MP_TAC DET_ROW_SPAN THEN ASM_REWRITE_TAC[] THEN
 719     ASM_REWRITE_TAC[] THEN MATCH_MP_TAC SPAN_VSUM THEN
 720     REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG; FINITE_DELETE; IN_DELETE] THEN
 721     REPEAT STRIP_TAC THEN MATCH_MP_TAC SPAN_MUL THEN
 722     MATCH_MP_TAC(CONJUNCT1 SPAN_CLAUSES) THEN
 723     REWRITE_TAC[IN_ELIM_THM] THEN ASM_MESON_TAC[]]);;
 724
 725 let CRAMER_LEMMA = prove
 726  (`!A:real^N^N x:real^N.
 727         1 <= k /\ k <= dimindex(:N)
 728         ==> det((lambda i j. if j = k then (A**x)$i else A$i$j):real^N^N) =
 729             x$k * det(A)`,
 730   REPEAT GEN_TAC THEN DISCH_TAC THEN REWRITE_TAC[MATRIX_MUL_VSUM] THEN
 731   FIRST_ASSUM(MP_TAC o SYM o SPECL [`transp(A:real^N^N)`; `x:real^N`] o
 732               MATCH_MP CRAMER_LEMMA_TRANSP) THEN
 733   REWRITE_TAC[DET_TRANSP] THEN DISCH_THEN SUBST1_TAC THEN
 734   GEN_REWRITE_TAC LAND_CONV [GSYM DET_TRANSP] THEN AP_TERM_TAC THEN
 735   ASM_SIMP_TAC[CART_EQ; transp; LAMBDA_BETA; MATRIX_MUL_VSUM; row; column;
 736         COND_COMPONENT; VECTOR_MUL_COMPONENT; VSUM_COMPONENT]);;
 737
 738 let CRAMER = prove
 739  (`!A:real^N^N x b.
 740         ~(det(A) = &0)
 741         ==> (A ** x = b <=>
 742              x = lambda k.
 743                    det((lambda i j. if j = k then b$i else A$i$j):real^N^N) /
 744                    det(A))`,
 745   GEN_TAC THEN REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN DISCH_TAC THEN
 746   ONCE_REWRITE_TAC[SWAP_FORALL_THM] THEN GEN_TAC THEN MATCH_MP_TAC(MESON[]
 747    `(?x. p(x)) /\ (!x. p(x) ==> x = a) ==> !x. p(x) <=> x = a`) THEN
 748   CONJ_TAC THENL
 749    [MP_TAC(SPEC `A:real^N^N` INVERTIBLE_DET_NZ) THEN
 750     ASM_MESON_TAC[invertible; MATRIX_VECTOR_MUL_ASSOC; MATRIX_VECTOR_MUL_LID];
 751     GEN_TAC THEN DISCH_THEN(SUBST1_TAC o SYM) THEN
 752     ASM_SIMP_TAC[CART_EQ; CRAMER_LEMMA; LAMBDA_BETA; REAL_FIELD
 753     `~(z = &0) ==> (x = y / z <=> x * z = y)`]]);;
 754
 755 (* ------------------------------------------------------------------------- *)
 756 (* Variants of Cramer's rule for matrix-matrix multiplication.               *)
 757 (* ------------------------------------------------------------------------- *)
 758
 759 let CRAMER_MATRIX_LEFT = prove
 760  (`!A:real^N^N X:real^N^N B:real^N^N.
 761         ~(det A = &0)
 762         ==> (X ** A = B <=>
 763              X = lambda k l.
 764                    det((lambda i j. if j = l then B$k$i else A$j$i):real^N^N) /
 765                    det A)`,
 766   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[CART_EQ] THEN
 767   ASM_SIMP_TAC[MATRIX_MUL_COMPONENT; CRAMER; DET_TRANSP] THEN
 768   SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 769   REPLICATE_TAC 2 (AP_TERM_TAC THEN ABS_TAC THEN AP_TERM_TAC) THEN
 770   AP_TERM_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 771   SIMP_TAC[CART_EQ; LAMBDA_BETA; transp]);;
 772
 773 let CRAMER_MATRIX_RIGHT = prove
 774  (`!A:real^N^N X:real^N^N B:real^N^N.
 775         ~(det A = &0)
 776         ==> (A ** X = B <=>
 777              X = lambda k l.
 778                    det((lambda i j. if j = k then B$i$l else A$i$j):real^N^N) /
 779                    det A)`,
 780   REPEAT STRIP_TAC THEN
 781   GEN_REWRITE_TAC LAND_CONV [GSYM TRANSP_EQ] THEN
 782   REWRITE_TAC[MATRIX_TRANSP_MUL] THEN
 783   ASM_SIMP_TAC[CRAMER_MATRIX_LEFT; DET_TRANSP] THEN
 784   GEN_REWRITE_TAC LAND_CONV [GSYM TRANSP_EQ] THEN
 785   REWRITE_TAC[TRANSP_TRANSP] THEN AP_TERM_TAC THEN
 786   SIMP_TAC[CART_EQ; LAMBDA_BETA; transp] THEN
 787   REPEAT(GEN_TAC THEN STRIP_TAC) THEN
 788   AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 789   SIMP_TAC[CART_EQ; LAMBDA_BETA; transp]);;
 790
 791 let CRAMER_MATRIX_RIGHT_INVERSE = prove
 792  (`!A:real^N^N A':real^N^N.
 793         A ** A' = mat 1 <=>
 794         ~(det A = &0) /\
 795         A' = lambda k l.
 796                 det((lambda i j. if j = k then if i = l then &1 else &0
 797                                  else A$i$j):real^N^N) /
 798                 det A`,
 799   REPEAT GEN_TAC THEN ASM_CASES_TAC `det(A:real^N^N) = &0` THENL
 800    [ASM_REWRITE_TAC[] THEN
 801     DISCH_THEN(MP_TAC o AP_TERM `det:real^N^N->real`) THEN
 802     ASM_REWRITE_TAC[DET_MUL; DET_I] THEN REAL_ARITH_TAC;
 803     ASM_SIMP_TAC[CRAMER_MATRIX_RIGHT] THEN AP_TERM_TAC THEN
 804     SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 805     REPEAT(GEN_TAC THEN STRIP_TAC) THEN
 806     AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 807     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA; mat]]);;
 808
 809 let CRAMER_MATRIX_LEFT_INVERSE = prove
 810  (`!A:real^N^N A':real^N^N.
 811         A' ** A = mat 1 <=>
 812         ~(det A = &0) /\
 813         A' = lambda k l.
 814                 det((lambda i j. if j = l then if i = k then &1 else &0
 815                                  else A$j$i):real^N^N) /
 816                 det A`,
 817   REPEAT GEN_TAC THEN ASM_CASES_TAC `det(A:real^N^N) = &0` THENL
 818    [ASM_REWRITE_TAC[] THEN
 819     DISCH_THEN(MP_TAC o AP_TERM `det:real^N^N->real`) THEN
 820     ASM_REWRITE_TAC[DET_MUL; DET_I] THEN REAL_ARITH_TAC;
 821     ASM_SIMP_TAC[CRAMER_MATRIX_LEFT] THEN AP_TERM_TAC THEN
 822     SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 823     REPEAT(GEN_TAC THEN STRIP_TAC) THEN
 824     AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 825     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA; mat] THEN MESON_TAC[]]);;
 826
 827 (* ------------------------------------------------------------------------- *)
 828 (* Cofactors and their relationship to inverse matrices.                     *)
 829 (* ------------------------------------------------------------------------- *)
 830
 831 let cofactor = new_definition
 832   `(cofactor:real^N^N->real^N^N) A =
 833         lambda i j. det((lambda k l. if k = i /\ l = j then &1
 834                                      else if k = i \/ l = j then &0
 835                                      else A$k$l):real^N^N)`;;
 836
 837 let COFACTOR_TRANSP = prove
 838  (`!A:real^N^N. cofactor(transp A) = transp(cofactor A)`,
 839   SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA; transp] THEN REPEAT STRIP_TAC THEN
 840   GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN
 841   AP_TERM_TAC THEN SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA; transp] THEN
 842   MESON_TAC[]);;
 843
 844 let COFACTOR_COLUMN = prove
 845  (`!A:real^N^N.
 846         cofactor A =
 847         lambda i j. det((lambda k l. if l = j then if k = i then &1 else &0
 848                                      else A$k$l):real^N^N)`,
 849   GEN_TAC THEN CONV_TAC SYM_CONV THEN
 850   SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA] THEN
 851   X_GEN_TAC `i:num` THEN STRIP_TAC THEN
 852   X_GEN_TAC `j:num` THEN STRIP_TAC THEN
 853   REWRITE_TAC[det] THEN MATCH_MP_TAC SUM_EQ THEN
 854   REWRITE_TAC[FORALL_IN_GSPEC] THEN GEN_TAC THEN
 855   DISCH_TAC THEN AP_TERM_TAC THEN
 856   ASM_CASES_TAC `(p:num->num) i = j` THENL
 857    [MATCH_MP_TAC PRODUCT_EQ THEN
 858     X_GEN_TAC `k:num` THEN SIMP_TAC[IN_NUMSEG; LAMBDA_BETA] THEN STRIP_TAC THEN
 859     SUBGOAL_THEN `(p:num->num) k IN 1..dimindex(:N)` MP_TAC THENL
 860      [ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG];
 861       SIMP_TAC[LAMBDA_BETA; IN_NUMSEG] THEN STRIP_TAC] THEN
 862     ASM_CASES_TAC `(p:num->num) k = j` THEN ASM_REWRITE_TAC[] THEN
 863     COND_CASES_TAC THEN ASM_REWRITE_TAC[] THEN ASM_MESON_TAC[];
 864     MATCH_MP_TAC(REAL_ARITH `s = &0 /\ t = &0 ==> s = t`) THEN
 865     ASM_SIMP_TAC[PRODUCT_EQ_0; FINITE_NUMSEG] THEN CONJ_TAC THEN
 866     EXISTS_TAC `inverse (p:num->num) j` THEN
 867     ASM_SIMP_TAC[IN_NUMSEG; LAMBDA_BETA] THEN
 868     (SUBGOAL_THEN `inverse(p:num->num) j IN 1..dimindex(:N)` MP_TAC THENL
 869       [ASM_MESON_TAC[PERMUTES_IN_IMAGE; PERMUTES_INVERSE; IN_NUMSEG];
 870        SIMP_TAC[LAMBDA_BETA; IN_NUMSEG] THEN STRIP_TAC] THEN
 871      SUBGOAL_THEN `(p:num->num)(inverse p j) = j` SUBST1_TAC THENL
 872       [ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG];
 873        ASM_SIMP_TAC[LAMBDA_BETA] THEN
 874         ASM_MESON_TAC[PERMUTES_INVERSE_EQ]])]);;
 875
 876 let COFACTOR_ROW = prove
 877  (`!A:real^N^N.
 878         cofactor A =
 879         lambda i j. det((lambda k l. if k = i then if l = j then &1 else &0
 880                                      else A$k$l):real^N^N)`,
 881   GEN_TAC THEN ONCE_REWRITE_TAC[GSYM TRANSP_EQ] THEN
 882   REWRITE_TAC[GSYM COFACTOR_TRANSP] THEN
 883   SIMP_TAC[COFACTOR_COLUMN; CART_EQ; LAMBDA_BETA; transp] THEN
 884   REPEAT STRIP_TAC THEN
 885   GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN
 886   AP_TERM_TAC THEN SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA; transp]);;
 887
 888 let MATRIX_RIGHT_INVERSE_COFACTOR = prove
 889  (`!A:real^N^N A':real^N^N.
 890         A ** A' = mat 1 <=>
 891         ~(det A = &0) /\ A' = inv(det A) %% transp(cofactor A)`,
 892   REPEAT GEN_TAC THEN REWRITE_TAC[CRAMER_MATRIX_RIGHT_INVERSE] THEN
 893   ASM_CASES_TAC `det(A:real^N^N) = &0` THEN ASM_REWRITE_TAC[] THEN
 894   AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA; MATRIX_CMUL_COMPONENT] THEN
 895   X_GEN_TAC `k:num` THEN STRIP_TAC THEN
 896   X_GEN_TAC `l:num` THEN STRIP_TAC THEN
 897   REWRITE_TAC[ONCE_REWRITE_RULE[REAL_MUL_SYM] real_div] THEN AP_TERM_TAC THEN
 898   ASM_SIMP_TAC[transp; COFACTOR_COLUMN; LAMBDA_BETA] THEN
 899   AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA]);;
 900
 901 let MATRIX_LEFT_INVERSE_COFACTOR = prove
 902  (`!A:real^N^N A':real^N^N.
 903         A' ** A = mat 1 <=>
 904         ~(det A = &0) /\ A' = inv(det A) %% transp(cofactor A)`,
 905   REPEAT GEN_TAC THEN
 906   ONCE_REWRITE_TAC[MATRIX_LEFT_RIGHT_INVERSE] THEN
 907   REWRITE_TAC[MATRIX_RIGHT_INVERSE_COFACTOR]);;
 908
 909 let MATRIX_INV_COFACTOR = prove
 910  (`!A. ~(det A = &0) ==> matrix_inv A = inv(det A) %% transp(cofactor A)`,
 911   GEN_TAC THEN DISCH_THEN(MP_TAC o MATCH_MP MATRIX_MUL_LINV) THEN
 912   SIMP_TAC[MATRIX_LEFT_INVERSE_COFACTOR]);;
 913
 914 let COFACTOR_MATRIX_INV = prove
 915  (`!A:real^N^N. ~(det A = &0) ==> cofactor A = det(A) %% transp(matrix_inv A)`,
 916   SIMP_TAC[MATRIX_INV_COFACTOR; TRANSP_MATRIX_CMUL; TRANSP_TRANSP] THEN
 917   SIMP_TAC[MATRIX_CMUL_ASSOC; REAL_MUL_RINV; MATRIX_CMUL_LID]);;
 918
 919 let COFACTOR_I = prove
 920  (`cofactor(mat 1:real^N^N) = mat 1`,
 921   SIMP_TAC[COFACTOR_MATRIX_INV; DET_I; REAL_OF_NUM_EQ; ARITH_EQ] THEN
 922   REWRITE_TAC[MATRIX_INV_I; MATRIX_CMUL_LID; TRANSP_MAT]);;
 923
 924 let DET_COFACTOR_EXPANSION = prove
 925  (`!A:real^N^N i.
 926         1 <= i /\ i <= dimindex(:N)
 927         ==> det A = sum (1..dimindex(:N))
 928                         (\j. A$i$j * (cofactor A)$i$j)`,
 929   REPEAT STRIP_TAC THEN ASM_SIMP_TAC[COFACTOR_COLUMN; LAMBDA_BETA; det] THEN
 930   REWRITE_TAC[GSYM SUM_LMUL] THEN
 931   W(MP_TAC o PART_MATCH (lhand o rand) SUM_SWAP o rand o snd) THEN
 932   ANTS_TAC THENL [SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG]; ALL_TAC] THEN
 933   DISCH_THEN SUBST1_TAC THEN
 934   MATCH_MP_TAC SUM_EQ THEN REWRITE_TAC[FORALL_IN_GSPEC] THEN
 935   GEN_TAC THEN DISCH_TAC THEN
 936   ONCE_REWRITE_TAC[REAL_ARITH `a * s * p:real = s * a * p`] THEN
 937   REWRITE_TAC[SUM_LMUL] THEN AP_TERM_TAC THEN MATCH_MP_TAC EQ_TRANS THEN
 938   EXISTS_TAC
 939    `sum (1..dimindex (:N))
 940         (\j. (A:real^N^N)$i$j *
 941              product
 942               (inverse p j INSERT ((1..dimindex(:N)) DELETE (inverse p j)))
 943               (\k. if k = inverse p j then if k = i then &1 else &0
 944                    else A$k$(p k)))` THEN
 945   CONJ_TAC THENL
 946    [SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_PERMUTATIONS;
 947              FINITE_NUMSEG; IN_DELETE] THEN
 948     SUBGOAL_THEN `!j. inverse (p:num->num) j = i <=> j = p i`
 949      (fun th -> REWRITE_TAC[th])
 950     THENL [ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG]; ALL_TAC] THEN
 951     REWRITE_TAC[REAL_ARITH
 952      `x * (if p then &1 else &0) * y = if p then x * y else &0`] THEN
 953     SIMP_TAC[SUM_DELTA] THEN COND_CASES_TAC THENL
 954      [ALL_TAC; ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]] THEN
 955     SUBGOAL_THEN
 956      `1..dimindex(:N) = i INSERT ((1..dimindex(:N)) DELETE i)`
 957      (fun th -> GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [th])
 958     THENL
 959      [ASM_SIMP_TAC[IN_NUMSEG; SET_RULE `s = x INSERT (s DELETE x) <=> x IN s`];
 960       SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE] THEN
 961       AP_TERM_TAC THEN MATCH_MP_TAC(MESON[PRODUCT_EQ]
 962        `s = t /\ (!x. x IN t ==> f x = g x) ==> product s f = product t g`) THEN
 963       SIMP_TAC[IN_DELETE] THEN ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG]];
 964     MATCH_MP_TAC SUM_EQ_NUMSEG THEN X_GEN_TAC `j:num` THEN STRIP_TAC THEN
 965     REWRITE_TAC[] THEN AP_TERM_TAC THEN MATCH_MP_TAC(MESON[PRODUCT_EQ]
 966      `s = t /\ (!x. x IN t ==> f x = g x) ==> product s f = product t g`) THEN
 967     CONJ_TAC THENL
 968      [REWRITE_TAC[SET_RULE `x INSERT (s DELETE x) = s <=> x IN s`] THEN
 969       ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; PERMUTES_INVERSE];
 970       X_GEN_TAC `k:num` THEN REWRITE_TAC[IN_NUMSEG] THEN STRIP_TAC THEN
 971       SUBGOAL_THEN `(p:num->num) k IN 1..dimindex(:N)` MP_TAC THENL
 972        [ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]; ALL_TAC] THEN
 973       SIMP_TAC[LAMBDA_BETA; IN_NUMSEG] THEN
 974       ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG]]]);;
 975
 976 let MATRIX_MUL_RIGHT_COFACTOR = prove
 977  (`!A:real^N^N. A ** transp(cofactor A) = det(A) %% mat 1`,
 978   GEN_TAC THEN
 979   SIMP_TAC[CART_EQ; MATRIX_CMUL_COMPONENT; mat;
 980            matrix_mul; LAMBDA_BETA; transp] THEN
 981   X_GEN_TAC `i:num` THEN STRIP_TAC THEN
 982   X_GEN_TAC `i':num` THEN STRIP_TAC THEN
 983   COND_CASES_TAC THEN
 984   ASM_SIMP_TAC[GSYM DET_COFACTOR_EXPANSION; REAL_MUL_RID] THEN
 985   MATCH_MP_TAC EQ_TRANS THEN
 986   EXISTS_TAC `det((lambda k l. if k = i' then (A:real^N^N)$i$l
 987                                else A$k$l):real^N^N)` THEN
 988   CONJ_TAC THENL
 989    [MP_TAC(GEN `A:real^N^N`
 990      (ISPECL [`A:real^N^N`; `i':num`] DET_COFACTOR_EXPANSION)) THEN
 991     ASM_REWRITE_TAC[] THEN ASM_SIMP_TAC[] THEN DISCH_THEN(K ALL_TAC) THEN
 992     MATCH_MP_TAC SUM_EQ THEN X_GEN_TAC `j:num` THEN
 993     REWRITE_TAC[IN_NUMSEG] THEN STRIP_TAC THEN
 994     ASM_SIMP_TAC[LAMBDA_BETA] THEN AP_TERM_TAC THEN
 995     ASM_SIMP_TAC[cofactor; LAMBDA_BETA] THEN AP_TERM_TAC THEN
 996     SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN ASM_MESON_TAC[];
 997     REWRITE_TAC[REAL_MUL_RZERO] THEN MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 998     MAP_EVERY EXISTS_TAC [`i:num`;` i':num`] THEN
 999     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA; row]]);;
1000
1001 let MATRIX_MUL_LEFT_COFACTOR = prove
1002  (`!A:real^N^N. transp(cofactor A) ** A = det(A) %% mat 1`,
1003   GEN_TAC THEN ONCE_REWRITE_TAC[GSYM TRANSP_EQ] THEN
1004   REWRITE_TAC[MATRIX_TRANSP_MUL] THEN
1005   ONCE_REWRITE_TAC[GSYM COFACTOR_TRANSP] THEN
1006   REWRITE_TAC[MATRIX_MUL_RIGHT_COFACTOR; TRANSP_MATRIX_CMUL] THEN
1007   REWRITE_TAC[DET_TRANSP; TRANSP_MAT]);;
1008
1009 let COFACTOR_CMUL = prove
1010  (`!A:real^N^N c. cofactor(c %% A) = c pow (dimindex(:N) - 1) %% cofactor A`,
1011   REPEAT GEN_TAC THEN
1012   SIMP_TAC[CART_EQ; cofactor; LAMBDA_BETA; MATRIX_CMUL_COMPONENT] THEN
1013   X_GEN_TAC `i:num` THEN STRIP_TAC THEN
1014   X_GEN_TAC `j:num` THEN STRIP_TAC THEN
1015   REWRITE_TAC[det; GSYM SUM_LMUL] THEN
1016   MATCH_MP_TAC SUM_EQ THEN REWRITE_TAC[FORALL_IN_GSPEC] THEN
1017   X_GEN_TAC `p:num->num` THEN DISCH_TAC THEN
1018   ONCE_REWRITE_TAC[REAL_ARITH `a * b * c:real = b * a * c`] THEN
1019   AP_TERM_TAC THEN
1020   SUBGOAL_THEN
1021    `1..dimindex (:N) = i INSERT ((1..dimindex (:N)) DELETE i)`
1022   SUBST1_TAC THENL
1023    [REWRITE_TAC[EXTENSION; IN_INSERT; IN_NUMSEG; IN_DELETE] THEN ASM_ARITH_TAC;
1024     SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE]] THEN
1025   SUBGOAL_THEN
1026    `1 <= (p:num->num) i /\ p i <= dimindex(:N)`
1027   ASSUME_TAC THENL
1028    [FIRST_ASSUM(MP_TAC o MATCH_MP PERMUTES_IMAGE) THEN
1029     REWRITE_TAC[EXTENSION; IN_IMAGE; IN_NUMSEG] THEN ASM SET_TAC[];
1030     ASM_SIMP_TAC[LAMBDA_BETA]] THEN
1031   COND_CASES_TAC THEN ASM_REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO] THEN
1032   SUBGOAL_THEN
1033    `dimindex(:N) - 1 = CARD((1..dimindex(:N)) DELETE i)`
1034   SUBST1_TAC THENL
1035    [ASM_SIMP_TAC[CARD_DELETE; FINITE_NUMSEG; IN_NUMSEG; CARD_NUMSEG_1];
1036     ASM_SIMP_TAC[REAL_MUL_LID; GSYM PRODUCT_CONST; FINITE_NUMSEG;
1037                  FINITE_DELETE; GSYM PRODUCT_MUL]] THEN
1038   MATCH_MP_TAC PRODUCT_EQ THEN
1039   X_GEN_TAC `k:num` THEN REWRITE_TAC[IN_DELETE; IN_NUMSEG] THEN STRIP_TAC THEN
1040   SUBGOAL_THEN
1041    `1 <= (p:num->num) k /\ p k <= dimindex(:N)`
1042   ASSUME_TAC THENL
1043    [FIRST_ASSUM(MP_TAC o MATCH_MP PERMUTES_IMAGE) THEN
1044     REWRITE_TAC[EXTENSION; IN_IMAGE; IN_NUMSEG] THEN ASM SET_TAC[];
1045     ASM_SIMP_TAC[LAMBDA_BETA] THEN REAL_ARITH_TAC]);;
1046
1047 let COFACTOR_0 = prove
1048  (`cofactor(mat 0:real^N^N) = if dimindex(:N) = 1 then mat 1 else mat 0`,
1049   MP_TAC(ISPECL [`mat 1:real^N^N`; `&0`] COFACTOR_CMUL) THEN
1050   REWRITE_TAC[MATRIX_CMUL_LZERO; COFACTOR_I; REAL_POW_ZERO] THEN
1051   DISCH_THEN SUBST1_TAC THEN
1052   SIMP_TAC[DIMINDEX_GE_1; ARITH_RULE `1 <= n ==> (n - 1 = 0 <=> n = 1)`] THEN
1053   COND_CASES_TAC THEN REWRITE_TAC[MATRIX_CMUL_LZERO; MATRIX_CMUL_LID]);;
1054
1055 (* ------------------------------------------------------------------------- *)
1056 (* Explicit formulas for low dimensions.                                     *)
1057 (* ------------------------------------------------------------------------- *)
1058
1059 let PRODUCT_1 = prove
1060  (`product(1..1) f = f(1)`,
1061   REWRITE_TAC[PRODUCT_SING_NUMSEG]);;
1062
1063 let PRODUCT_2 = prove
1064  (`!t. product(1..2) t = t(1) * t(2)`,
1065   REWRITE_TAC[num_CONV `2`; PRODUCT_CLAUSES_NUMSEG] THEN
1066   REWRITE_TAC[PRODUCT_SING_NUMSEG; ARITH; REAL_MUL_ASSOC]);;
1067
1068 let PRODUCT_3 = prove
1069  (`!t. product(1..3) t = t(1) * t(2) * t(3)`,
1070   REWRITE_TAC[num_CONV `3`; num_CONV `2`; PRODUCT_CLAUSES_NUMSEG] THEN
1071   REWRITE_TAC[PRODUCT_SING_NUMSEG; ARITH; REAL_MUL_ASSOC]);;
1072
1073 let PRODUCT_4 = prove
1074  (`!t. product(1..4) t = t(1) * t(2) * t(3) * t(4)`,
1075   REWRITE_TAC[num_CONV `4`; num_CONV `3`; num_CONV `2`;
1076               PRODUCT_CLAUSES_NUMSEG] THEN
1077   REWRITE_TAC[PRODUCT_SING_NUMSEG; ARITH; REAL_MUL_ASSOC]);;
1078
1079 let DET_1 = prove
1080  (`!A:real^1^1. det A = A$1$1`,
1081   REWRITE_TAC[det; DIMINDEX_1; PERMUTES_SING; NUMSEG_SING] THEN
1082   REWRITE_TAC[SUM_SING; SET_RULE `{x | x = a} = {a}`; PRODUCT_SING] THEN
1083   REWRITE_TAC[SIGN_I; I_THM] THEN REAL_ARITH_TAC);;
1084
1085 let DET_2 = prove
1086  (`!A:real^2^2. det A = A$1$1 * A$2$2 - A$1$2 * A$2$1`,
1087   GEN_TAC THEN REWRITE_TAC[det; DIMINDEX_2] THEN
1088   CONV_TAC(LAND_CONV(RATOR_CONV(ONCE_DEPTH_CONV NUMSEG_CONV))) THEN
1089   SIMP_TAC[SUM_OVER_PERMUTATIONS_INSERT; FINITE_INSERT; FINITE_EMPTY;
1090            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1091   REWRITE_TAC[PERMUTES_EMPTY; SUM_SING; SET_RULE `{x | x = a} = {a}`] THEN
1092   REWRITE_TAC[SWAP_REFL; I_O_ID] THEN
1093   REWRITE_TAC[GSYM(NUMSEG_CONV `1..2`); SUM_2] THEN
1094   SIMP_TAC[SUM_CLAUSES; FINITE_INSERT; FINITE_EMPTY;
1095            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1096   SIMP_TAC[SIGN_COMPOSE; PERMUTATION_SWAP] THEN
1097   REWRITE_TAC[SIGN_SWAP; ARITH] THEN REWRITE_TAC[PRODUCT_2] THEN
1098   REWRITE_TAC[o_THM; swap; ARITH] THEN REAL_ARITH_TAC);;
1099
1100 let DET_3 = prove
1101  (`!A:real^3^3.
1102         det(A) = A$1$1 * A$2$2 * A$3$3 +
1103                  A$1$2 * A$2$3 * A$3$1 +
1104                  A$1$3 * A$2$1 * A$3$2 -
1105                  A$1$1 * A$2$3 * A$3$2 -
1106                  A$1$2 * A$2$1 * A$3$3 -
1107                  A$1$3 * A$2$2 * A$3$1`,
1108   GEN_TAC THEN REWRITE_TAC[det; DIMINDEX_3] THEN
1109   CONV_TAC(LAND_CONV(RATOR_CONV(ONCE_DEPTH_CONV NUMSEG_CONV))) THEN
1110   SIMP_TAC[SUM_OVER_PERMUTATIONS_INSERT; FINITE_INSERT; FINITE_EMPTY;
1111            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1112   REWRITE_TAC[PERMUTES_EMPTY; SUM_SING; SET_RULE `{x | x = a} = {a}`] THEN
1113   REWRITE_TAC[SWAP_REFL; I_O_ID] THEN
1114   REWRITE_TAC[GSYM(NUMSEG_CONV `1..3`); SUM_3] THEN
1115   SIMP_TAC[SUM_CLAUSES; FINITE_INSERT; FINITE_EMPTY;
1116            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1117   SIMP_TAC[SIGN_COMPOSE; PERMUTATION_SWAP] THEN
1118   REWRITE_TAC[SIGN_SWAP; ARITH] THEN REWRITE_TAC[PRODUCT_3] THEN
1119   REWRITE_TAC[o_THM; swap; ARITH] THEN REAL_ARITH_TAC);;
1120
1121 let DET_4 = prove
1122  (`!A:real^4^4.
1123         det(A) = A$1$1 * A$2$2 * A$3$3 * A$4$4 +
1124                  A$1$1 * A$2$3 * A$3$4 * A$4$2 +
1125                  A$1$1 * A$2$4 * A$3$2 * A$4$3 +
1126                  A$1$2 * A$2$1 * A$3$4 * A$4$3 +
1127                  A$1$2 * A$2$3 * A$3$1 * A$4$4 +
1128                  A$1$2 * A$2$4 * A$3$3 * A$4$1 +
1129                  A$1$3 * A$2$1 * A$3$2 * A$4$4 +
1130                  A$1$3 * A$2$2 * A$3$4 * A$4$1 +
1131                  A$1$3 * A$2$4 * A$3$1 * A$4$2 +
1132                  A$1$4 * A$2$1 * A$3$3 * A$4$2 +
1133                  A$1$4 * A$2$2 * A$3$1 * A$4$3 +
1134                  A$1$4 * A$2$3 * A$3$2 * A$4$1 -
1135                  A$1$1 * A$2$2 * A$3$4 * A$4$3 -
1136                  A$1$1 * A$2$3 * A$3$2 * A$4$4 -
1137                  A$1$1 * A$2$4 * A$3$3 * A$4$2 -
1138                  A$1$2 * A$2$1 * A$3$3 * A$4$4 -
1139                  A$1$2 * A$2$3 * A$3$4 * A$4$1 -
1140                  A$1$2 * A$2$4 * A$3$1 * A$4$3 -
1141                  A$1$3 * A$2$1 * A$3$4 * A$4$2 -
1142                  A$1$3 * A$2$2 * A$3$1 * A$4$4 -
1143                  A$1$3 * A$2$4 * A$3$2 * A$4$1 -
1144                  A$1$4 * A$2$1 * A$3$2 * A$4$3 -
1145                  A$1$4 * A$2$2 * A$3$3 * A$4$1 -
1146                  A$1$4 * A$2$3 * A$3$1 * A$4$2`,
1147   let lemma = prove
1148    (`(sum {3,4} f = f 3 + f 4) /\
1149      (sum {2,3,4} f = f 2 + f 3 + f 4)`,
1150     SIMP_TAC[SUM_CLAUSES; FINITE_INSERT; FINITE_EMPTY] THEN
1151     REWRITE_TAC[ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN REAL_ARITH_TAC) in
1152   GEN_TAC THEN REWRITE_TAC[det; DIMINDEX_4] THEN
1153   CONV_TAC(LAND_CONV(RATOR_CONV(ONCE_DEPTH_CONV NUMSEG_CONV))) THEN
1154   SIMP_TAC[SUM_OVER_PERMUTATIONS_INSERT; FINITE_INSERT; FINITE_EMPTY;
1155            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1156   REWRITE_TAC[PERMUTES_EMPTY; SUM_SING; SET_RULE `{x | x = a} = {a}`] THEN
1157   REWRITE_TAC[SWAP_REFL; I_O_ID] THEN
1158   REWRITE_TAC[GSYM(NUMSEG_CONV `1..4`); SUM_4; lemma] THEN
1159   SIMP_TAC[SIGN_COMPOSE; PERMUTATION_SWAP; PERMUTATION_COMPOSE] THEN
1160   REWRITE_TAC[SIGN_SWAP; ARITH] THEN REWRITE_TAC[PRODUCT_4] THEN
1161   REWRITE_TAC[o_THM; swap; ARITH] THEN REAL_ARITH_TAC);;
1162
1163 (* ------------------------------------------------------------------------- *)
1164 (* Existence of the characteristic polynomial.                               *)
1165 (* ------------------------------------------------------------------------- *)
1166
1167 let CHARACTERISTIC_POLYNOMIAL = prove
1168  (`!A:real^N^N.
1169         ?a. a(dimindex(:N)) = &1 /\
1170             !x. det(x %% mat 1 - A) =
1171                 sum (0..dimindex(:N)) (\i. a i * x pow i)`,
1172   GEN_TAC THEN REWRITE_TAC[det] THEN
1173   SUBGOAL_THEN
1174    `!p n. IMAGE p (1..dimindex(:N)) SUBSET 1..dimindex(:N) /\
1175           n <= dimindex(:N)
1176           ==> ?a. a n = (if !i. 1 <= i /\ i <= n ==> p i = i then &1 else &0) /\
1177                   !x. product (1..n) (\i. (x %% mat 1 - A:real^N^N)$i$p i) =
1178                       sum (0..n) (\i. a i * x pow i)`
1179   MP_TAC THENL
1180    [GEN_TAC THEN REWRITE_TAC[IMP_CONJ; RIGHT_FORALL_IMP_THM] THEN
1181     DISCH_TAC THEN INDUCT_TAC THEN
1182     REWRITE_TAC[PRODUCT_CLAUSES_NUMSEG] THEN
1183     REWRITE_TAC[LE_0; ARITH_EQ; ARITH_RULE `1 <= SUC n`] THENL
1184      [EXISTS_TAC `\i. if i = 0 then &1 else &0` THEN
1185       SIMP_TAC[real_pow; REAL_MUL_LID; ARITH_RULE `1 <= i ==> ~(i <= 0)`;
1186                SUM_CLAUSES_NUMSEG];
1187       DISCH_TAC THEN FIRST_X_ASSUM(MP_TAC o check (is_imp o concl)) THEN
1188       ASM_SIMP_TAC[ARITH_RULE `SUC n <= N ==> n <= N`] THEN
1189       DISCH_THEN(X_CHOOSE_THEN `a:num->real` STRIP_ASSUME_TAC) THEN
1190       ASM_REWRITE_TAC[] THEN
1191       REWRITE_TAC[MATRIX_SUB_COMPONENT; MATRIX_CMUL_COMPONENT] THEN
1192       ASSUME_TAC(ARITH_RULE `1 <= SUC n`) THEN
1193       FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [SUBSET]) THEN
1194       REWRITE_TAC[FORALL_IN_IMAGE; IN_NUMSEG] THEN
1195       DISCH_THEN(MP_TAC o SPEC `SUC n`) THEN ASM_REWRITE_TAC[] THEN
1196       STRIP_TAC THEN ASM_SIMP_TAC[MAT_COMPONENT] THEN
1197       ASM_CASES_TAC `p(SUC n) = SUC n` THEN ASM_REWRITE_TAC[] THENL
1198        [ALL_TAC;
1199         EXISTS_TAC `\i. if i <= n
1200                         then --((A:real^N^N)$(SUC n)$(p(SUC n))) * a i
1201                         else &0` THEN
1202         SIMP_TAC[SUM_CLAUSES_NUMSEG; LE_0; ARITH_RULE `~(SUC n <= n)`] THEN
1203         CONJ_TAC THENL
1204          [COND_CASES_TAC THEN REWRITE_TAC[] THEN
1205           FIRST_X_ASSUM(MP_TAC o SPEC `SUC n`) THEN
1206           ASM_REWRITE_TAC[] THEN ASM_ARITH_TAC;
1207           REWRITE_TAC[REAL_MUL_LZERO; REAL_ADD_RID; GSYM SUM_RMUL] THEN
1208           GEN_TAC THEN MATCH_MP_TAC SUM_EQ_NUMSEG THEN REWRITE_TAC[] THEN
1209           REAL_ARITH_TAC]] THEN
1210       REWRITE_TAC[REAL_SUB_LDISTRIB; REAL_MUL_RID] THEN
1211       REWRITE_TAC[GSYM SUM_RMUL] THEN EXISTS_TAC
1212       `\i. (if i = 0 then &0 else a(i - 1)) -
1213            (if i = SUC n then &0 else (A:real^N^N)$(SUC n)$(SUC n) * a i)` THEN
1214       ASM_REWRITE_TAC[NOT_SUC; LE; SUC_SUB1; REAL_SUB_RZERO] THEN
1215       CONJ_TAC THENL [ASM_MESON_TAC[LE_REFL]; ALL_TAC] THEN
1216       REWRITE_TAC[REAL_SUB_RDISTRIB; SUM_SUB_NUMSEG] THEN
1217       GEN_TAC THEN BINOP_TAC THENL
1218        [SIMP_TAC[SUM_CLAUSES_LEFT; ARITH_RULE `0 <= SUC n`] THEN
1219         REWRITE_TAC[ADD1; SUM_OFFSET; ARITH_RULE `~(i + 1 = 0)`; ADD_SUB] THEN
1220         REWRITE_TAC[REAL_MUL_LZERO; REAL_POW_ADD; REAL_POW_1; REAL_ADD_LID];
1221         SIMP_TAC[SUM_CLAUSES_NUMSEG; LE_0; REAL_MUL_LZERO; REAL_ADD_RID] THEN
1222         SIMP_TAC[ARITH_RULE `i <= n ==> ~(i = SUC n)`]] THEN
1223       MATCH_MP_TAC SUM_EQ_NUMSEG THEN REWRITE_TAC[REAL_ADD_LID; REAL_MUL_AC]];
1224     GEN_REWRITE_TAC LAND_CONV [SWAP_FORALL_THM] THEN
1225     DISCH_THEN(MP_TAC o SPEC `dimindex(:N)`) THEN REWRITE_TAC[LE_REFL] THEN
1226     GEN_REWRITE_TAC (LAND_CONV o BINDER_CONV) [RIGHT_IMP_EXISTS_THM] THEN
1227     REWRITE_TAC[SKOLEM_THM; LEFT_IMP_EXISTS_THM] THEN
1228     X_GEN_TAC `a:(num->num)->num->real` THEN DISCH_TAC] THEN
1229   EXISTS_TAC
1230    `\i:num. sum {p | p permutes 1..dimindex(:N)} (\p. sign p * a p i)` THEN
1231   REWRITE_TAC[] THEN CONJ_TAC THENL
1232    [MP_TAC(ISPECL
1233      [`\p:num->num. sign p * a p (dimindex(:N))`;
1234       `{p | p permutes 1..dimindex(:N)}`;
1235       `I:num->num`] SUM_DELETE) THEN
1236     SIMP_TAC[IN_ELIM_THM; PERMUTES_I; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
1237     MATCH_MP_TAC(REAL_ARITH `k = &1 /\ s' = &0 ==> s' = s - k ==> s = &1`) THEN
1238     CONJ_TAC THENL
1239      [FIRST_X_ASSUM(MP_TAC o SPEC `I:num->num`) THEN
1240       SIMP_TAC[IMAGE_I; SUBSET_REFL; SIGN_I; I_THM; REAL_MUL_LID];
1241       MATCH_MP_TAC SUM_EQ_0 THEN X_GEN_TAC `p:num->num` THEN
1242       REWRITE_TAC[IN_ELIM_THM; IN_DELETE] THEN STRIP_TAC THEN
1243       FIRST_X_ASSUM(MP_TAC o SPEC `p:num->num`) THEN ANTS_TAC THENL
1244        [ASM_MESON_TAC[PERMUTES_IMAGE; SUBSET_REFL]; ALL_TAC] THEN
1245       COND_CASES_TAC THEN SIMP_TAC[REAL_MUL_RZERO] THEN
1246       FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [permutes]) THEN
1247       FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV [FUN_EQ_THM]) THEN
1248       REWRITE_TAC[IN_NUMSEG; I_THM] THEN ASM_MESON_TAC[]];
1249     X_GEN_TAC `x:real` THEN REWRITE_TAC[GSYM SUM_RMUL] THEN
1250     W(MP_TAC o PART_MATCH (lhs o rand) SUM_SWAP o rand o snd) THEN
1251     SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
1252     DISCH_THEN SUBST1_TAC THEN MATCH_MP_TAC SUM_EQ THEN
1253     X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
1254     REWRITE_TAC[GSYM REAL_MUL_ASSOC; SUM_LMUL] THEN AP_TERM_TAC THEN
1255     FIRST_X_ASSUM(MP_TAC o SPEC `p:num->num`) THEN ANTS_TAC THENL
1256      [ASM_MESON_TAC[PERMUTES_IMAGE; SUBSET_REFL]; SIMP_TAC[]]]);;
1257
1258 (* ------------------------------------------------------------------------- *)
1259 (* Grassmann-Plucker relations for n = 2, n = 3 and n = 4.                   *)
1260 (* I have a proof of the general n case but the proof is a bit long and the  *)
1261 (* result doesn't seem generally useful enough to go in the main theories.   *)
1262 (* ------------------------------------------------------------------------- *)
1263
1264 let GRASSMANN_PLUCKER_2 = prove
1265  (`!x1 x2 y1 y2:real^2.
1266         det(vector[x1;x2]) * det(vector[y1;y2]) =
1267           det(vector[y1;x2]) * det(vector[x1;y2]) +
1268           det(vector[y2;x2]) * det(vector[y1;x1])`,
1269   REWRITE_TAC[DET_2; VECTOR_2] THEN REAL_ARITH_TAC);;
1270
1271 let GRASSMANN_PLUCKER_3 = prove
1272  (`!x1 x2 x3 y1 y2 y3:real^3.
1273         det(vector[x1;x2;x3]) * det(vector[y1;y2;y3]) =
1274           det(vector[y1;x2;x3]) * det(vector[x1;y2;y3]) +
1275           det(vector[y2;x2;x3]) * det(vector[y1;x1;y3]) +
1276           det(vector[y3;x2;x3]) * det(vector[y1;y2;x1])`,
1277   REWRITE_TAC[DET_3; VECTOR_3] THEN REAL_ARITH_TAC);;
1278
1279 let GRASSMANN_PLUCKER_4 = prove
1280  (`!x1 x2 x3 x4:real^4 y1 y2 y3 y4:real^4.
1281         det(vector[x1;x2;x3;x4]) * det(vector[y1;y2;y3;y4]) =
1282           det(vector[y1;x2;x3;x4]) * det(vector[x1;y2;y3;y4]) +
1283           det(vector[y2;x2;x3;x4]) * det(vector[y1;x1;y3;y4]) +
1284           det(vector[y3;x2;x3;x4]) * det(vector[y1;y2;x1;y4]) +
1285           det(vector[y4;x2;x3;x4]) * det(vector[y1;y2;y3;x1])`,
1286   REWRITE_TAC[DET_4; VECTOR_4] THEN REAL_ARITH_TAC);;
1287
1288 (* ------------------------------------------------------------------------- *)
1289 (* Determinants of integer matrices.                                         *)
1290 (* ------------------------------------------------------------------------- *)
1291
1292 let INTEGER_PRODUCT = prove
1293  (`!f s. (!x. x IN s ==> integer(f x)) ==> integer(product s f)`,
1294   REPEAT STRIP_TAC THEN MATCH_MP_TAC PRODUCT_CLOSED THEN
1295   ASM_REWRITE_TAC[INTEGER_CLOSED]);;
1296
1297 let INTEGER_SIGN = prove
1298  (`!p. integer(sign p)`,
1299   SIMP_TAC[sign; COND_RAND; INTEGER_CLOSED; COND_ID]);;
1300
1301 let INTEGER_DET = prove
1302  (`!M:real^N^N.
1303         (!i j. 1 <= i /\ i <= dimindex(:N) /\
1304                1 <= j /\ j <= dimindex(:N)
1305                ==> integer(M$i$j))
1306         ==> integer(det M)`,
1307   REPEAT STRIP_TAC THEN REWRITE_TAC[det] THEN
1308   MATCH_MP_TAC INTEGER_SUM THEN X_GEN_TAC `p:num->num` THEN
1309   REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
1310   MATCH_MP_TAC INTEGER_MUL THEN REWRITE_TAC[INTEGER_SIGN] THEN
1311   MATCH_MP_TAC INTEGER_PRODUCT THEN REWRITE_TAC[IN_NUMSEG] THEN
1312   REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN
1313   ASM_MESON_TAC[IN_NUMSEG; permutes]);;
1314
1315 (* ------------------------------------------------------------------------- *)
1316 (* Diagonal matrices (for arbitrary rectangular matrix, not just square).    *)
1317 (* ------------------------------------------------------------------------- *)
1318
1319 let diagonal_matrix = new_definition
1320  `diagonal_matrix(A:real^N^M) <=>
1321         !i j. 1 <= i /\ i <= dimindex(:M) /\
1322               1 <= j /\ j <= dimindex(:N) /\
1323               ~(i = j)
1324               ==> A$i$j = &0`;;
1325
1326 let TRANSP_DIAGONAL_MATRIX = prove
1327  (`!A:real^N^N. diagonal_matrix A ==> transp A = A`,
1328   GEN_TAC THEN REWRITE_TAC[diagonal_matrix; CART_EQ; TRANSP_COMPONENT] THEN
1329   STRIP_TAC THEN X_GEN_TAC `i:num` THEN STRIP_TAC THEN X_GEN_TAC `j:num` THEN
1330   STRIP_TAC THEN ASM_CASES_TAC `i:num = j` THEN ASM_SIMP_TAC[]);;
1331
1332 (* ------------------------------------------------------------------------- *)
1333 (* Orthogonality of a transformation and matrix.                             *)
1334 (* ------------------------------------------------------------------------- *)
1335
1336 let orthogonal_transformation = new_definition
1337  `orthogonal_transformation(f:real^N->real^N) <=>
1338         linear f /\ !v w. f(v) dot f(w) = v dot w`;;
1339
1340 let ORTHOGONAL_TRANSFORMATION = prove
1341  (`!f. orthogonal_transformation f <=> linear f /\ !v. norm(f v) = norm(v)`,
1342   GEN_TAC THEN REWRITE_TAC[orthogonal_transformation] THEN EQ_TAC THENL
1343    [MESON_TAC[vector_norm]; SIMP_TAC[DOT_NORM] THEN MESON_TAC[LINEAR_ADD]]);;
1344
1345 let ORTHOGONAL_TRANSFORMATION_COMPOSE = prove
1346  (`!f g. orthogonal_transformation f /\ orthogonal_transformation g
1347          ==> orthogonal_transformation(f o g)`,
1348   SIMP_TAC[orthogonal_transformation; LINEAR_COMPOSE; o_THM]);;
1349
1350 let orthogonal_matrix = new_definition
1351  `orthogonal_matrix(Q:real^N^N) <=>
1352       transp(Q) ** Q = mat 1 /\ Q ** transp(Q) = mat 1`;;
1353
1354 let ORTHOGONAL_MATRIX = prove
1355  (`orthogonal_matrix(Q:real^N^N) <=> transp(Q) ** Q = mat 1`,
1356   MESON_TAC[MATRIX_LEFT_RIGHT_INVERSE; orthogonal_matrix]);;
1357
1358 let ORTHOGONAL_MATRIX_ALT = prove
1359  (`!A:real^N^N. orthogonal_matrix A <=> A ** transp A = mat 1`,
1360   MESON_TAC[MATRIX_LEFT_RIGHT_INVERSE; orthogonal_matrix]);;
1361
1362 let ORTHOGONAL_MATRIX_ID = prove
1363  (`orthogonal_matrix(mat 1)`,
1364   REWRITE_TAC[orthogonal_matrix; TRANSP_MAT; MATRIX_MUL_LID]);;
1365
1366 let ORTHOGONAL_MATRIX_MUL = prove
1367  (`!A B. orthogonal_matrix A /\ orthogonal_matrix B
1368          ==> orthogonal_matrix(A ** B)`,
1369   REWRITE_TAC[orthogonal_matrix; MATRIX_TRANSP_MUL] THEN
1370   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[GSYM MATRIX_MUL_ASSOC] THEN
1371   GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [MATRIX_MUL_ASSOC] THEN
1372   ASM_REWRITE_TAC[MATRIX_MUL_LID; MATRIX_MUL_RID]);;
1373
1374 let ORTHOGONAL_TRANSFORMATION_MATRIX = prove
1375  (`!f:real^N->real^N.
1376      orthogonal_transformation f <=> linear f /\ orthogonal_matrix(matrix f)`,
1377   REPEAT STRIP_TAC THEN EQ_TAC THENL
1378    [REWRITE_TAC[orthogonal_transformation; ORTHOGONAL_MATRIX] THEN
1379     STRIP_TAC THEN ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
1380     X_GEN_TAC `i:num` THEN STRIP_TAC THEN
1381     X_GEN_TAC `j:num` THEN STRIP_TAC THEN
1382     FIRST_X_ASSUM(MP_TAC o SPECL [`basis i:real^N`; `basis j:real^N`]) THEN
1383     FIRST_ASSUM(fun th -> REWRITE_TAC[GSYM(MATCH_MP MATRIX_WORKS th)]) THEN
1384     REWRITE_TAC[DOT_MATRIX_VECTOR_MUL] THEN
1385     ABBREV_TAC `A = transp (matrix f) ** matrix(f:real^N->real^N)` THEN
1386     ASM_SIMP_TAC[matrix_mul; columnvector; rowvector; basis; LAMBDA_BETA;
1387              SUM_DELTA; DIMINDEX_1; LE_REFL; dot; IN_NUMSEG; mat;
1388              MESON[REAL_MUL_LID; REAL_MUL_LZERO; REAL_MUL_RID; REAL_MUL_RZERO]
1389               `(if b then &1 else &0) * x = (if b then x else &0) /\
1390                x * (if b then &1 else &0) = (if b then x else &0)`];
1391     REWRITE_TAC[orthogonal_matrix; ORTHOGONAL_TRANSFORMATION; NORM_EQ] THEN
1392     DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
1393     FIRST_ASSUM(fun th -> REWRITE_TAC[GSYM(MATCH_MP MATRIX_WORKS th)]) THEN
1394     ASM_REWRITE_TAC[DOT_MATRIX_VECTOR_MUL] THEN
1395     SIMP_TAC[DOT_MATRIX_PRODUCT; MATRIX_MUL_LID]]);;
1396
1397 let ORTHOGONAL_MATRIX_TRANSFORMATION = prove
1398  (`!A:real^N^N. orthogonal_matrix A <=> orthogonal_transformation(\x. A ** x)`,
1399   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX; MATRIX_VECTOR_MUL_LINEAR] THEN
1400   REWRITE_TAC[MATRIX_OF_MATRIX_VECTOR_MUL]);;
1401
1402 let ORTHOGONAL_MATRIX_MATRIX = prove
1403  (`!f:real^N->real^N.
1404     orthogonal_transformation f ==> orthogonal_matrix(matrix f)`,
1405   SIMP_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX]);;
1406
1407 let DET_ORTHOGONAL_MATRIX = prove
1408  (`!Q. orthogonal_matrix Q ==> det(Q) = &1 \/ det(Q) = -- &1`,
1409   GEN_TAC THEN REWRITE_TAC[REAL_RING `x = &1 \/ x = -- &1 <=> x * x = &1`] THEN
1410   GEN_REWRITE_TAC (RAND_CONV o LAND_CONV o RAND_CONV) [GSYM DET_TRANSP] THEN
1411   SIMP_TAC[GSYM DET_MUL; orthogonal_matrix; DET_I]);;
1412
1413 let ORTHOGONAL_MATRIX_TRANSP = prove
1414  (`!A:real^N^N. orthogonal_matrix(transp A) <=> orthogonal_matrix A`,
1415   REWRITE_TAC[orthogonal_matrix; TRANSP_TRANSP; CONJ_ACI]);;
1416
1417 let MATRIX_MUL_LTRANSP_DOT_COLUMN = prove
1418  (`!A:real^N^M. transp A ** A = (lambda i j. (column i A) dot (column j A))`,
1419   SIMP_TAC[matrix_mul; CART_EQ; LAMBDA_BETA; transp; dot; column]);;
1420
1421 let MATRIX_MUL_RTRANSP_DOT_ROW = prove
1422  (`!A:real^N^M. A ** transp A = (lambda i j. (row i A) dot (row j A))`,
1423   SIMP_TAC[matrix_mul; CART_EQ; LAMBDA_BETA; transp; dot; row]);;
1424
1425 let ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS = prove
1426  (`!A:real^N^N.
1427         orthogonal_matrix A <=>
1428         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(column i A) = &1) /\
1429         (!i j. 1 <= i /\ i <= dimindex(:N) /\
1430                1 <= j /\ j <= dimindex(:N) /\ ~(i = j)
1431                ==> orthogonal (column i A) (column j A))`,
1432   REWRITE_TAC[ORTHOGONAL_MATRIX] THEN
1433   SIMP_TAC[MATRIX_MUL_LTRANSP_DOT_COLUMN; CART_EQ; mat; LAMBDA_BETA] THEN
1434   REWRITE_TAC[orthogonal; NORM_EQ_1] THEN MESON_TAC[]);;
1435
1436 let ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS = prove
1437  (`!A:real^N^N.
1438         orthogonal_matrix A <=>
1439         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(row i A) = &1) /\
1440         (!i j. 1 <= i /\ i <= dimindex(:N) /\
1441                1 <= j /\ j <= dimindex(:N) /\ ~(i = j)
1442                ==> orthogonal (row i A) (row j A))`,
1443   ONCE_REWRITE_TAC[GSYM ORTHOGONAL_MATRIX_TRANSP] THEN
1444   SIMP_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS; COLUMN_TRANSP]);;
1445
1446 let ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_INDEXED = prove
1447  (`!A:real^N^N.
1448         orthogonal_matrix A <=>
1449         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(row i A) = &1) /\
1450         pairwise (\i j. orthogonal (row i A) (row j A)) (1..dimindex(:N))`,
1451   REPEAT GEN_TAC THEN REWRITE_TAC[ORTHOGONAL_MATRIX_ALT] THEN
1452   SIMP_TAC[CART_EQ; LAMBDA_BETA; pairwise; MAT_COMPONENT] THEN
1453   SIMP_TAC[MATRIX_MUL_RTRANSP_DOT_ROW; IN_NUMSEG; LAMBDA_BETA] THEN
1454   REWRITE_TAC[NORM_EQ_SQUARE; REAL_POS] THEN
1455   CONV_TAC REAL_RAT_REDUCE_CONV THEN REWRITE_TAC[orthogonal] THEN
1456   MESON_TAC[]);;
1457
1458 let ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_PAIRWISE = prove
1459  (`!A:real^N^N.
1460         orthogonal_matrix A <=>
1461         CARD(rows A) = dimindex(:N) /\
1462         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(row i A) = &1) /\
1463         pairwise orthogonal (rows A)`,
1464   REWRITE_TAC[rows; ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_INDEXED] THEN
1465   GEN_TAC THEN ONCE_REWRITE_TAC[SIMPLE_IMAGE_GEN] THEN
1466   REWRITE_TAC[PAIRWISE_IMAGE; GSYM numseg] THEN
1467   MATCH_MP_TAC(TAUT `(p ==> (q <=> r /\ s)) ==> (p /\ q <=> r /\ p /\ s)`) THEN
1468   DISCH_TAC THEN GEN_REWRITE_TAC (RAND_CONV o LAND_CONV o RAND_CONV)
1469     [GSYM CARD_NUMSEG_1] THEN
1470   SIMP_TAC[CARD_IMAGE_EQ_INJ; FINITE_NUMSEG] THEN
1471   REWRITE_TAC[pairwise; IN_NUMSEG] THEN
1472   ASM_MESON_TAC[ORTHOGONAL_REFL; NORM_ARITH `~(norm(vec 0:real^N) = &1)`]);;
1473
1474 let ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_SPAN = prove
1475  (`!A:real^N^N.
1476         orthogonal_matrix A <=>
1477         span(rows A) = (:real^N) /\
1478         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(row i A) = &1) /\
1479         pairwise orthogonal (rows A)`,
1480   GEN_TAC THEN REWRITE_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_PAIRWISE] THEN
1481   EQ_TAC THEN STRIP_TAC THEN ASM_REWRITE_TAC[] THENL
1482    [MATCH_MP_TAC(SET_RULE `UNIV SUBSET s ==> s = UNIV`) THEN
1483     MATCH_MP_TAC CARD_GE_DIM_INDEPENDENT THEN
1484     ASM_REWRITE_TAC[DIM_UNIV; SUBSET_UNIV; LE_REFL];
1485     CONV_TAC SYM_CONV THEN REWRITE_TAC[GSYM DIM_UNIV] THEN
1486     FIRST_X_ASSUM(SUBST1_TAC o SYM) THEN REWRITE_TAC[DIM_SPAN] THEN
1487     MATCH_MP_TAC DIM_EQ_CARD] THEN
1488   MATCH_MP_TAC PAIRWISE_ORTHOGONAL_INDEPENDENT THEN
1489   ASM_REWRITE_TAC[rows; IN_ELIM_THM] THEN
1490   ASM_MESON_TAC[NORM_ARITH `~(norm(vec 0:real^N) = &1)`]);;
1491
1492 let ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS_INDEXED = prove
1493  (`!A:real^N^N.
1494       orthogonal_matrix A <=>
1495       (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(column i A) = &1) /\
1496       pairwise (\i j. orthogonal (column i A) (column j A)) (1..dimindex(:N))`,
1497   ONCE_REWRITE_TAC[GSYM ORTHOGONAL_MATRIX_TRANSP] THEN
1498   REWRITE_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_INDEXED] THEN
1499   SIMP_TAC[ROW_TRANSP; ROWS_TRANSP; pairwise; IN_NUMSEG]);;
1500
1501 let ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS_PAIRWISE = prove
1502  (`!A:real^N^N.
1503         orthogonal_matrix A <=>
1504         CARD(columns A) = dimindex(:N) /\
1505         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(column i A) = &1) /\
1506         pairwise orthogonal (columns A)`,
1507   ONCE_REWRITE_TAC[GSYM ORTHOGONAL_MATRIX_TRANSP] THEN
1508   REWRITE_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_PAIRWISE] THEN
1509   SIMP_TAC[ROW_TRANSP; ROWS_TRANSP]);;
1510
1511 let ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS_SPAN = prove
1512  (`!A:real^N^N.
1513         orthogonal_matrix A <=>
1514         span(columns A) = (:real^N) /\
1515         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(column i A) = &1) /\
1516         pairwise orthogonal (columns A)`,
1517   ONCE_REWRITE_TAC[GSYM ORTHOGONAL_MATRIX_TRANSP] THEN
1518   REWRITE_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS_SPAN] THEN
1519   SIMP_TAC[ROW_TRANSP; ROWS_TRANSP]);;
1520
1521 let ORTHOGONAL_MATRIX_2 = prove
1522  (`!A:real^2^2. orthogonal_matrix A <=>
1523                 A$1$1 pow 2 + A$2$1 pow 2 = &1 /\
1524                 A$1$2 pow 2 + A$2$2 pow 2 = &1 /\
1525                 A$1$1 * A$1$2 + A$2$1 * A$2$2 = &0`,
1526   SIMP_TAC[orthogonal_matrix; CART_EQ; matrix_mul; LAMBDA_BETA;
1527            TRANSP_COMPONENT; MAT_COMPONENT] THEN
1528   REWRITE_TAC[DIMINDEX_2; FORALL_2; SUM_2] THEN
1529   CONV_TAC NUM_REDUCE_CONV THEN CONV_TAC REAL_RING);;
1530
1531 let ORTHOGONAL_MATRIX_2_ALT = prove
1532  (`!A:real^2^2. orthogonal_matrix A <=>
1533                 A$1$1 pow 2 + A$2$1 pow 2 = &1 /\
1534                 (A$1$1 = A$2$2 /\ A$1$2 = --(A$2$1) \/
1535                  A$1$1 = --(A$2$2) /\ A$1$2 = A$2$1)`,
1536   REWRITE_TAC[ORTHOGONAL_MATRIX_2] THEN CONV_TAC REAL_RING);;
1537
1538 let ORTHOGONAL_MATRIX_INV = prove
1539  (`!A:real^N^N. orthogonal_matrix A ==> matrix_inv A = transp A`,
1540   MESON_TAC[orthogonal_matrix; MATRIX_INV_UNIQUE]);;
1541
1542 let ORTHOGONAL_TRANSFORMATION_ORTHOGONAL_EIGENVECTORS = prove
1543  (`!f:real^N->real^N v w a b.
1544         orthogonal_transformation f /\ f v = a % v /\ f w = b % w /\ ~(a = b)
1545         ==> orthogonal v w`,
1546   REWRITE_TAC[orthogonal_transformation] THEN REPEAT STRIP_TAC THEN
1547   FIRST_X_ASSUM(fun th ->
1548     MP_TAC(SPECL [`v:real^N`; `v:real^N`] th) THEN
1549     MP_TAC(SPECL [`v:real^N`; `w:real^N`] th) THEN
1550     MP_TAC(SPECL [`w:real^N`; `w:real^N`] th)) THEN
1551   ASM_REWRITE_TAC[DOT_LMUL; DOT_RMUL; orthogonal] THEN
1552   REWRITE_TAC[REAL_MUL_ASSOC; REAL_RING `x * y = y <=> x = &1 \/ y = &0`] THEN
1553   REWRITE_TAC[DOT_EQ_0] THEN
1554   ASM_CASES_TAC `v:real^N = vec 0` THEN ASM_REWRITE_TAC[DOT_LZERO] THEN
1555   ASM_CASES_TAC `w:real^N = vec 0` THEN ASM_REWRITE_TAC[DOT_RZERO] THEN
1556   ASM_CASES_TAC `(v:real^N) dot w = &0` THEN ASM_REWRITE_TAC[] THEN
1557   UNDISCH_TAC `~(a:real = b)` THEN CONV_TAC REAL_RING);;
1558
1559 let ORTHOGONAL_MATRIX_ORTHOGONAL_EIGENVECTORS = prove
1560  (`!A:real^N^N v w a b.
1561         orthogonal_matrix A /\ A ** v = a % v /\ A ** w = b % w /\ ~(a = b)
1562         ==> orthogonal v w`,
1563   REWRITE_TAC[ORTHOGONAL_MATRIX_TRANSFORMATION;
1564               ORTHOGONAL_TRANSFORMATION_ORTHOGONAL_EIGENVECTORS]);;
1565
1566 (* ------------------------------------------------------------------------- *)
1567 (* Linearity of scaling, and hence isometry, that preserves origin.          *)
1568 (* ------------------------------------------------------------------------- *)
1569
1570 let SCALING_LINEAR = prove
1571  (`!f:real^M->real^N c.
1572         (f(vec 0) = vec 0) /\ (!x y. dist(f x,f y) = c * dist(x,y))
1573         ==> linear(f)`,
1574   REPEAT STRIP_TAC THEN
1575   SUBGOAL_THEN `!v w. ((f:real^M->real^N) v) dot (f w) = c pow 2 * (v dot w)`
1576   ASSUME_TAC THENL
1577    [FIRST_ASSUM(MP_TAC o GEN `v:real^M` o
1578       SPECL [`v:real^M`; `vec 0 :real^M`]) THEN
1579     REWRITE_TAC[dist] THEN ASM_REWRITE_TAC[VECTOR_SUB_RZERO] THEN
1580     DISCH_TAC THEN ASM_REWRITE_TAC[DOT_NORM_NEG; GSYM dist] THEN
1581     REAL_ARITH_TAC;
1582     ALL_TAC] THEN
1583   REWRITE_TAC[linear; VECTOR_EQ] THEN
1584   ASM_REWRITE_TAC[DOT_LADD; DOT_RADD; DOT_LMUL; DOT_RMUL] THEN
1585   REAL_ARITH_TAC);;
1586
1587 let ISOMETRY_LINEAR = prove
1588  (`!f:real^M->real^N.
1589         (f(vec 0) = vec 0) /\ (!x y. dist(f x,f y) = dist(x,y))
1590         ==> linear(f)`,
1591   MESON_TAC[SCALING_LINEAR; REAL_MUL_LID]);;
1592
1593 let ISOMETRY_IMP_AFFINITY = prove
1594  (`!f:real^M->real^N.
1595         (!x y. dist(f x,f y) = dist(x,y))
1596         ==> ?h. linear h /\ !x. f(x) = f(vec 0) + h(x)`,
1597   REPEAT STRIP_TAC THEN
1598   EXISTS_TAC `\x. (f:real^M->real^N) x - f(vec 0)` THEN
1599   REWRITE_TAC[VECTOR_ARITH `a + (x - a):real^N = x`] THEN
1600   MATCH_MP_TAC ISOMETRY_LINEAR THEN REWRITE_TAC[VECTOR_SUB_REFL] THEN
1601   ASM_REWRITE_TAC[NORM_ARITH `dist(x - a:real^N,y - a) = dist(x,y)`]);;
1602
1603 (* ------------------------------------------------------------------------- *)
1604 (* Hence another formulation of orthogonal transformation.                   *)
1605 (* ------------------------------------------------------------------------- *)
1606
1607 let ORTHOGONAL_TRANSFORMATION_ISOMETRY = prove
1608  (`!f:real^N->real^N.
1609         orthogonal_transformation f <=>
1610         (f(vec 0) = vec 0) /\ (!x y. dist(f x,f y) = dist(x,y))`,
1611   GEN_TAC THEN REWRITE_TAC[ORTHOGONAL_TRANSFORMATION] THEN EQ_TAC THENL
1612    [MESON_TAC[LINEAR_0; LINEAR_SUB; dist]; STRIP_TAC] THEN
1613   ASM_SIMP_TAC[ISOMETRY_LINEAR] THEN X_GEN_TAC `x:real^N` THEN
1614   FIRST_X_ASSUM(MP_TAC o SPECL [`x:real^N`; `vec 0:real^N`]) THEN
1615   ASM_REWRITE_TAC[dist; VECTOR_SUB_RZERO]);;
1616
1617 (* ------------------------------------------------------------------------- *)
1618 (* Can extend an isometry from unit sphere.                                  *)
1619 (* ------------------------------------------------------------------------- *)
1620
1621 let ISOMETRY_SPHERE_EXTEND = prove
1622  (`!f:real^N->real^N.
1623         (!x. norm(x) = &1 ==> norm(f x) = &1) /\
1624         (!x y. norm(x) = &1 /\ norm(y) = &1 ==> dist(f x,f y) = dist(x,y))
1625         ==> ?g. orthogonal_transformation g /\
1626                 (!x. norm(x) = &1 ==> g(x) = f(x))`,
1627   let lemma = prove
1628    (`!x:real^N y:real^N x':real^N y':real^N x0 y0 x0' y0'.
1629           x = norm(x) % x0 /\ y = norm(y) % y0 /\
1630           x' = norm(x) % x0' /\ y' = norm(y) % y0' /\
1631           norm(x0) = &1 /\ norm(x0') = &1 /\ norm(y0) = &1 /\ norm(y0') = &1 /\
1632           norm(x0' - y0') = norm(x0 - y0)
1633           ==> norm(x' - y') = norm(x - y)`,
1634     REPEAT GEN_TAC THEN
1635     MAP_EVERY ABBREV_TAC [`a = norm(x:real^N)`; `b = norm(y:real^N)`] THEN
1636     REPLICATE_TAC 4 (DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC)) THEN
1637     ASM_REWRITE_TAC[] THEN REWRITE_TAC[NORM_EQ; NORM_EQ_1] THEN
1638     REWRITE_TAC[DOT_LSUB; DOT_RSUB; DOT_LMUL; DOT_RMUL] THEN
1639     REWRITE_TAC[DOT_SYM] THEN CONV_TAC REAL_RING) in
1640   REPEAT STRIP_TAC THEN
1641   EXISTS_TAC `\x. if x = vec 0 then vec 0
1642                   else norm(x) % (f:real^N->real^N)(inv(norm x) % x)` THEN
1643   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_ISOMETRY] THEN
1644   SIMP_TAC[VECTOR_MUL_LID; REAL_INV_1] THEN CONJ_TAC THENL
1645    [ALL_TAC; MESON_TAC[NORM_0; REAL_ARITH `~(&1 = &0)`]] THEN
1646   REPEAT GEN_TAC THEN REPEAT(COND_CASES_TAC THEN ASM_REWRITE_TAC[]) THEN
1647   REWRITE_TAC[dist; VECTOR_SUB_LZERO; VECTOR_SUB_RZERO; NORM_NEG; NORM_MUL;
1648               REAL_ABS_NORM] THEN
1649   ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
1650   ASM_SIMP_TAC[GSYM REAL_EQ_RDIV_EQ; NORM_POS_LT] THEN
1651   ASM_SIMP_TAC[REAL_DIV_REFL; REAL_LT_IMP_NZ; NORM_EQ_0] THEN
1652   TRY(FIRST_X_ASSUM MATCH_MP_TAC) THEN
1653   REWRITE_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
1654   ASM_SIMP_TAC[REAL_MUL_LINV; NORM_EQ_0] THEN
1655   MATCH_MP_TAC lemma THEN MAP_EVERY EXISTS_TAC
1656    [`inv(norm x) % x:real^N`; `inv(norm y) % y:real^N`;
1657     `(f:real^N->real^N) (inv (norm x) % x)`;
1658     `(f:real^N->real^N) (inv (norm y) % y)`] THEN
1659   REWRITE_TAC[NORM_MUL; VECTOR_MUL_ASSOC; REAL_ABS_INV; REAL_ABS_NORM] THEN
1660   ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_RINV; NORM_EQ_0] THEN
1661   ASM_REWRITE_TAC[GSYM dist; VECTOR_MUL_LID] THEN
1662   REPEAT CONJ_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN
1663   REWRITE_TAC[NORM_MUL; VECTOR_MUL_ASSOC; REAL_ABS_INV; REAL_ABS_NORM] THEN
1664   ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_RINV; NORM_EQ_0]);;
1665
1666 let ORTHOGONAL_TRANSFORMATION_LINEAR = prove
1667  (`!f:real^N->real^N. orthogonal_transformation f ==> linear f`,
1668   SIMP_TAC[orthogonal_transformation]);;
1669
1670 let ORTHOGONAL_TRANSFORMATION_INJECTIVE = prove
1671  (`!f:real^N->real^N.
1672         orthogonal_transformation f ==> !x y. f x = f y ==> x = y`,
1673   SIMP_TAC[LINEAR_INJECTIVE_0; ORTHOGONAL_TRANSFORMATION; GSYM NORM_EQ_0]);;
1674
1675 let ORTHOGONAL_TRANSFORMATION_SURJECTIVE = prove
1676  (`!f:real^N->real^N.
1677         orthogonal_transformation f ==> !y. ?x. f x = y`,
1678   MESON_TAC[LINEAR_INJECTIVE_IMP_SURJECTIVE;
1679             ORTHOGONAL_TRANSFORMATION_INJECTIVE; orthogonal_transformation]);;
1680
1681 let ORTHOGONAL_TRANSFORMATION_INVERSE_o = prove
1682  (`!f:real^N->real^N.
1683         orthogonal_transformation f
1684         ==> ?g. orthogonal_transformation g /\ g o f = I /\ f o g = I`,
1685   REPEAT STRIP_TAC THEN
1686   FIRST_ASSUM(ASSUME_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
1687   FIRST_ASSUM(ASSUME_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_INJECTIVE) THEN
1688   MP_TAC(ISPEC `f:real^N->real^N` LINEAR_INJECTIVE_LEFT_INVERSE) THEN
1689   ASM_REWRITE_TAC[] THEN MATCH_MP_TAC MONO_EXISTS THEN
1690   X_GEN_TAC `g:real^N->real^N` THEN  STRIP_TAC THEN
1691   MP_TAC(ISPECL [`f:real^N->real^N`; `g:real^N->real^N`]
1692     LINEAR_INVERSE_LEFT) THEN
1693   ASM_REWRITE_TAC[] THEN DISCH_TAC THEN ASM_REWRITE_TAC[] THEN
1694   ASM_REWRITE_TAC[ORTHOGONAL_TRANSFORMATION] THEN X_GEN_TAC `v:real^N` THEN
1695   MATCH_MP_TAC EQ_TRANS THEN
1696   EXISTS_TAC `norm((f:real^N->real^N)((g:real^N->real^N) v))` THEN
1697   CONJ_TAC THENL [ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION]; ALL_TAC] THEN
1698   RULE_ASSUM_TAC(REWRITE_RULE[FUN_EQ_THM; o_THM; I_THM]) THEN
1699   ASM_REWRITE_TAC[]);;
1700
1701 let ORTHOGONAL_TRANSFORMATION_INVERSE = prove
1702  (`!f:real^N->real^N.
1703         orthogonal_transformation f
1704         ==> ?g. orthogonal_transformation g /\
1705                 (!x. g(f x) = x) /\ (!y. f(g y) = y)`,
1706   GEN_TAC THEN
1707   DISCH_THEN(MP_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_INVERSE_o) THEN
1708   REWRITE_TAC[FUN_EQ_THM; o_THM; I_THM]);;
1709
1710 let ORTHOGONAL_TRANSFORMATION_ID = prove
1711  (`orthogonal_transformation(\x. x)`,
1712   REWRITE_TAC[orthogonal_transformation; LINEAR_ID]);;
1713
1714 let ORTHOGONAL_TRANSFORMATION_I = prove
1715  (`orthogonal_transformation I`,
1716   REWRITE_TAC[I_DEF; ORTHOGONAL_TRANSFORMATION_ID]);;
1717
1718 (* ------------------------------------------------------------------------- *)
1719 (* We can find an orthogonal matrix taking any unit vector to any other.     *)
1720 (* ------------------------------------------------------------------------- *)
1721
1722 let FINITE_INDEX_NUMSEG_SPECIAL = prove
1723  (`!s a:A.
1724         FINITE s /\ a IN s
1725         ==> ?f. (!i j. i IN 1..CARD s /\ j IN 1..CARD s /\ f i = f j
1726                        ==> i = j) /\
1727                 s = IMAGE f (1..CARD s) /\
1728                 f 1 = a`,
1729   REPEAT STRIP_TAC THEN
1730   FIRST_ASSUM(MP_TAC o GEN_REWRITE_RULE I [FINITE_INDEX_NUMSEG]) THEN
1731   DISCH_THEN(X_CHOOSE_THEN `f:num->A` STRIP_ASSUME_TAC) THEN
1732   SUBGOAL_THEN `?k. k IN 1..CARD(s:A->bool) /\ (a:A) = f k`
1733   STRIP_ASSUME_TAC THENL[ASM SET_TAC[]; ALL_TAC] THEN
1734   EXISTS_TAC `(f:num->A) o swap(1,k)` THEN
1735   SUBGOAL_THEN `1 IN 1..CARD(s:A->bool)` ASSUME_TAC THENL
1736    [REWRITE_TAC[IN_NUMSEG; LE_REFL; ARITH_RULE `1 <= x <=> ~(x = 0)`] THEN
1737     ASM_SIMP_TAC[CARD_EQ_0; ARITH_EQ] THEN ASM SET_TAC[];
1738     ALL_TAC] THEN
1739   ASM_REWRITE_TAC[o_THM; swap] THEN
1740   CONJ_TAC THENL [ASM SET_TAC[]; ALL_TAC] THEN
1741   UNDISCH_THEN `s = IMAGE (f:num->A) (1..CARD(s:A->bool))`
1742    (fun th -> GEN_REWRITE_TAC LAND_CONV [th]) THEN
1743   REWRITE_TAC[EXTENSION; IN_IMAGE; o_THM] THEN
1744   X_GEN_TAC `b:A` THEN EQ_TAC THEN
1745   DISCH_THEN(X_CHOOSE_THEN `i:num` STRIP_ASSUME_TAC) THEN
1746   EXISTS_TAC `swap(1,k) i` THEN
1747   REWRITE_TAC[swap] THEN ASM_MESON_TAC[swap]);;
1748
1749 let ORTHOGONAL_MATRIX_EXISTS_BASIS = prove
1750  (`!a:real^N.
1751         norm(a) = &1
1752         ==> ?A. orthogonal_matrix A /\ A**(basis 1) = a`,
1753   REPEAT STRIP_TAC THEN
1754   FIRST_ASSUM(MP_TAC o MATCH_MP VECTOR_IN_ORTHONORMAL_BASIS) THEN
1755   REWRITE_TAC[HAS_SIZE] THEN
1756   DISCH_THEN(X_CHOOSE_THEN `s:real^N->bool` STRIP_ASSUME_TAC) THEN
1757   MP_TAC(ISPECL [`s:real^N->bool`; `a:real^N`]
1758    FINITE_INDEX_NUMSEG_SPECIAL) THEN ASM_REWRITE_TAC[IN_NUMSEG] THEN
1759   REWRITE_TAC[TAUT `a /\ b ==> c <=> c \/ ~a \/ ~b`] THEN
1760   DISCH_THEN(X_CHOOSE_THEN `f:num->real^N`
1761    (CONJUNCTS_THEN2 ASSUME_TAC (CONJUNCTS_THEN2 (ASSUME_TAC o SYM)
1762      ASSUME_TAC))) THEN
1763   EXISTS_TAC `(lambda i j. ((f j):real^N)$i):real^N^N` THEN
1764   SIMP_TAC[CART_EQ; LAMBDA_BETA; matrix_vector_mul; BASIS_COMPONENT;
1765            IN_NUMSEG] THEN
1766   ONCE_REWRITE_TAC[COND_RAND] THEN SIMP_TAC[REAL_MUL_RZERO; SUM_DELTA] THEN
1767   ASM_REWRITE_TAC[IN_NUMSEG; REAL_MUL_RID; LE_REFL; DIMINDEX_GE_1] THEN
1768   REWRITE_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS] THEN
1769   SIMP_TAC[column; LAMBDA_BETA] THEN CONJ_TAC THENL
1770    [X_GEN_TAC `i:num` THEN REPEAT STRIP_TAC THEN
1771     MATCH_MP_TAC EQ_TRANS THEN
1772     EXISTS_TAC `norm((f:num->real^N) i)` THEN CONJ_TAC THENL
1773      [AP_TERM_TAC THEN ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA];
1774       ASM_MESON_TAC[IN_IMAGE; IN_NUMSEG]];
1775     MAP_EVERY X_GEN_TAC [`i:num`; `j:num`] THEN STRIP_TAC THEN
1776     SUBGOAL_THEN `orthogonal ((f:num->real^N) i) (f j)` MP_TAC THENL
1777      [ASM_MESON_TAC[pairwise; IN_IMAGE; IN_NUMSEG]; ALL_TAC] THEN
1778     MATCH_MP_TAC EQ_IMP THEN BINOP_TAC THEN
1779     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA]]);;
1780
1781 let ORTHOGONAL_TRANSFORMATION_EXISTS_1 = prove
1782  (`!a b:real^N.
1783         norm(a) = &1 /\ norm(b) = &1
1784         ==> ?f. orthogonal_transformation f /\ f a = b`,
1785   REPEAT STRIP_TAC THEN
1786   MP_TAC(ISPEC `b:real^N` ORTHOGONAL_MATRIX_EXISTS_BASIS) THEN
1787   MP_TAC(ISPEC `a:real^N` ORTHOGONAL_MATRIX_EXISTS_BASIS) THEN
1788   ASM_REWRITE_TAC[] THEN
1789   DISCH_THEN(X_CHOOSE_THEN `A:real^N^N` (STRIP_ASSUME_TAC o GSYM)) THEN
1790   DISCH_THEN(X_CHOOSE_THEN `B:real^N^N` (STRIP_ASSUME_TAC o GSYM)) THEN
1791   EXISTS_TAC `\x:real^N. ((B:real^N^N) ** transp(A:real^N^N)) ** x` THEN
1792   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX; MATRIX_VECTOR_MUL_LINEAR;
1793               MATRIX_OF_MATRIX_VECTOR_MUL] THEN
1794   ASM_SIMP_TAC[ORTHOGONAL_MATRIX_MUL; ORTHOGONAL_MATRIX_TRANSP] THEN
1795   REWRITE_TAC[GSYM MATRIX_VECTOR_MUL_ASSOC] THEN AP_TERM_TAC THEN
1796   RULE_ASSUM_TAC(REWRITE_RULE[ORTHOGONAL_MATRIX]) THEN
1797   ASM_REWRITE_TAC[MATRIX_VECTOR_MUL_ASSOC; MATRIX_VECTOR_MUL_LID]);;
1798
1799 let ORTHOGONAL_TRANSFORMATION_EXISTS = prove
1800  (`!a b:real^N.
1801         norm(a) = norm(b) ==> ?f. orthogonal_transformation f /\ f a = b`,
1802   REPEAT GEN_TAC THEN ASM_CASES_TAC `b:real^N = vec 0` THEN
1803   ASM_SIMP_TAC[NORM_0; NORM_EQ_0] THENL
1804    [MESON_TAC[ORTHOGONAL_TRANSFORMATION_ID]; ALL_TAC] THEN
1805   ASM_CASES_TAC `a:real^N = vec 0` THENL
1806    [ASM_MESON_TAC[NORM_0; NORM_EQ_0]; ALL_TAC] THEN
1807   DISCH_TAC THEN
1808   MP_TAC(ISPECL [`inv(norm a) % a:real^N`; `inv(norm b) % b:real^N`]
1809                 ORTHOGONAL_TRANSFORMATION_EXISTS_1) THEN
1810   REWRITE_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
1811   ASM_SIMP_TAC[NORM_EQ_0; REAL_MUL_LINV] THEN
1812   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
1813   DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
1814   ASM_REWRITE_TAC[] THEN
1815   FIRST_ASSUM(ASSUME_TAC o MATCH_MP LINEAR_CMUL o
1816               MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
1817   ASM_REWRITE_TAC[VECTOR_ARITH
1818     `a % x:real^N = a % y <=> a % (x - y) = vec 0`] THEN
1819   ASM_REWRITE_TAC[VECTOR_MUL_EQ_0; REAL_INV_EQ_0; NORM_EQ_0; VECTOR_SUB_EQ]);;
1820
1821 (* ------------------------------------------------------------------------- *)
1822 (* Or indeed, taking any subspace to another of suitable dimension.          *)
1823 (* ------------------------------------------------------------------------- *)
1824
1825 let ORTHOGONAL_TRANSFORMATION_INTO_SUBSPACE = prove
1826  (`!s t:real^N->bool.
1827         subspace s /\ subspace t /\ dim s <= dim t
1828         ==> ?f. orthogonal_transformation f /\ IMAGE f s SUBSET t`,
1829   REPEAT STRIP_TAC THEN
1830   MP_TAC(ISPEC `t:real^N->bool` ORTHONORMAL_BASIS_SUBSPACE) THEN
1831   MP_TAC(ISPEC `s:real^N->bool` ORTHONORMAL_BASIS_SUBSPACE) THEN
1832   ASM_REWRITE_TAC[HAS_SIZE] THEN
1833   DISCH_THEN(X_CHOOSE_THEN `b:real^N->bool` STRIP_ASSUME_TAC) THEN
1834   DISCH_THEN(X_CHOOSE_THEN `c:real^N->bool` STRIP_ASSUME_TAC) THEN
1835   MP_TAC(ISPECL [`c:real^N->bool`; `(:real^N)`] ORTHONORMAL_EXTENSION) THEN
1836   MP_TAC(ISPECL [`b:real^N->bool`; `(:real^N)`] ORTHONORMAL_EXTENSION) THEN
1837   ASM_REWRITE_TAC[UNION_UNIV; SPAN_UNIV; LEFT_IMP_EXISTS_THM] THEN
1838   X_GEN_TAC `b':real^N->bool` THEN STRIP_TAC THEN
1839   X_GEN_TAC `c':real^N->bool` THEN STRIP_TAC THEN
1840   SUBGOAL_THEN
1841    `independent(b UNION b':real^N->bool) /\
1842     independent(c UNION c':real^N->bool)`
1843   STRIP_ASSUME_TAC THENL
1844    [CONJ_TAC THEN MATCH_MP_TAC PAIRWISE_ORTHOGONAL_INDEPENDENT THEN
1845     ASM_REWRITE_TAC[IN_UNION] THEN
1846     ASM_MESON_TAC[NORM_ARITH `~(norm(vec 0:real^N) = &1)`];
1847     ALL_TAC] THEN
1848   SUBGOAL_THEN `FINITE(b UNION b':real^N->bool) /\
1849                 FINITE(c UNION c':real^N->bool)`
1850   MP_TAC THENL
1851    [ASM_SIMP_TAC[PAIRWISE_ORTHOGONAL_IMP_FINITE];
1852     REWRITE_TAC[FINITE_UNION] THEN STRIP_TAC] THEN
1853   SUBGOAL_THEN
1854    `?f:real^N->real^N.
1855         (!x y. x IN b UNION b' /\ y IN b UNION b' ==> (f x = f y <=> x = y)) /\
1856         IMAGE f b SUBSET c /\
1857         IMAGE f (b UNION b') SUBSET c UNION c'`
1858    (X_CHOOSE_THEN `fb:real^N->real^N` STRIP_ASSUME_TAC)
1859   THENL
1860    [MP_TAC(ISPECL [`b:real^N->bool`; `c:real^N->bool`]
1861         CARD_LE_INJ) THEN
1862     ASM_REWRITE_TAC[LEFT_IMP_EXISTS_THM; INJECTIVE_ON_ALT] THEN
1863     X_GEN_TAC `f:real^N->real^N` THEN STRIP_TAC THEN
1864     MP_TAC(ISPECL [`b':real^N->bool`;
1865                    `(c UNION c') DIFF IMAGE (f:real^N->real^N) b`]
1866         CARD_LE_INJ) THEN
1867     ANTS_TAC THENL
1868      [ASM_SIMP_TAC[FINITE_UNION; FINITE_DIFF] THEN
1869       W(MP_TAC o PART_MATCH (lhs o rand) CARD_DIFF o rand o snd) THEN
1870       ASM_REWRITE_TAC[FINITE_UNION] THEN
1871       ANTS_TAC THENL [ASM SET_TAC[]; DISCH_THEN SUBST1_TAC] THEN
1872       MATCH_MP_TAC(ARITH_RULE `a + b:num = c ==> a <= c - b`) THEN
1873       W(MP_TAC o PART_MATCH (lhs o rand) CARD_IMAGE_INJ o
1874         rand o lhs o snd) THEN
1875       ANTS_TAC THENL [ASM_MESON_TAC[]; DISCH_THEN SUBST1_TAC] THEN
1876       W(MP_TAC o PART_MATCH (rhs o rand) CARD_UNION o lhs o snd) THEN
1877       ANTS_TAC THENL [ASM SET_TAC[]; DISCH_THEN(SUBST1_TAC o SYM)] THEN
1878       GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [UNION_COMM] THEN
1879       MATCH_MP_TAC(MESON[LE_ANTISYM]
1880        `(FINITE s /\ CARD s <= CARD t) /\
1881         (FINITE t /\ CARD t <= CARD s) ==> CARD s = CARD t`) THEN
1882       CONJ_TAC THEN MATCH_MP_TAC INDEPENDENT_SPAN_BOUND THEN
1883       ASM_REWRITE_TAC[FINITE_UNION; SUBSET_UNIV];
1884       DISCH_THEN(X_CHOOSE_THEN `g:real^N->real^N` STRIP_ASSUME_TAC) THEN
1885       EXISTS_TAC `\x. if x IN b then (f:real^N->real^N) x else g x` THEN
1886       REWRITE_TAC[SUBSET; FORALL_IN_IMAGE] THEN ASM SET_TAC[]];
1887     ALL_TAC] THEN
1888   MP_TAC(ISPECL [`fb:real^N->real^N`; `b UNION b':real^N->bool`]
1889     LINEAR_INDEPENDENT_EXTEND) THEN ASM_REWRITE_TAC[] THEN
1890   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
1891   STRIP_TAC THEN ASM_REWRITE_TAC[] THEN CONJ_TAC THENL
1892    [ASM_REWRITE_TAC[ORTHOGONAL_TRANSFORMATION];
1893     REWRITE_TAC[SYM(ASSUME `span b:real^N->bool = s`)] THEN
1894     ASM_SIMP_TAC[GSYM SPAN_LINEAR_IMAGE] THEN
1895     REWRITE_TAC[SYM(ASSUME `span c:real^N->bool = t`)] THEN
1896     MATCH_MP_TAC SPAN_MONO THEN ASM SET_TAC[]] THEN
1897   SUBGOAL_THEN
1898    `!v. v IN UNIV ==> norm((f:real^N->real^N) v) = norm v`
1899    (fun th -> ASM_MESON_TAC[th; IN_UNIV]) THEN
1900   UNDISCH_THEN `span (b UNION b') = (:real^N)` (SUBST1_TAC o SYM) THEN
1901   ASM_SIMP_TAC[SPAN_FINITE; FINITE_UNION; IN_ELIM_THM; LEFT_IMP_EXISTS_THM] THEN
1902   MAP_EVERY X_GEN_TAC [`z:real^N`; `u:real^N->real`] THEN
1903   DISCH_THEN(SUBST1_TAC o SYM) THEN ASM_SIMP_TAC[LINEAR_VSUM; FINITE_UNION] THEN
1904   REWRITE_TAC[o_DEF; NORM_EQ_SQUARE; NORM_POS_LE; GSYM NORM_POW_2] THEN
1905   ASM_SIMP_TAC[LINEAR_CMUL] THEN
1906   W(MP_TAC o PART_MATCH (lhand o rand)
1907     NORM_VSUM_PYTHAGOREAN o rand o snd) THEN
1908   W(MP_TAC o PART_MATCH (lhand o rand)
1909     NORM_VSUM_PYTHAGOREAN o lhand o rand o snd) THEN
1910   RULE_ASSUM_TAC(REWRITE_RULE[pairwise]) THEN
1911   ASM_SIMP_TAC[pairwise; ORTHOGONAL_CLAUSES; FINITE_UNION] THEN ANTS_TAC THENL
1912    [REPEAT STRIP_TAC THEN REWRITE_TAC[ORTHOGONAL_MUL] THEN
1913     REPEAT DISJ2_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM SET_TAC[];
1914     REPEAT(DISCH_THEN SUBST1_TAC) THEN ASM_SIMP_TAC[NORM_MUL] THEN
1915     MATCH_MP_TAC SUM_EQ THEN ASM SET_TAC[]]);;
1916
1917 let ORTHOGONAL_TRANSFORMATION_ONTO_SUBSPACE = prove
1918  (`!s t:real^N->bool.
1919         subspace s /\ subspace t /\ dim s = dim t
1920         ==> ?f. orthogonal_transformation f /\ IMAGE f s = t`,
1921   REPEAT STRIP_TAC THEN
1922   MP_TAC(ISPECL [`s:real^N->bool`; `t:real^N->bool`]
1923         ORTHOGONAL_TRANSFORMATION_INTO_SUBSPACE) THEN
1924   ASM_REWRITE_TAC[LE_REFL] THEN MATCH_MP_TAC MONO_EXISTS THEN
1925   X_GEN_TAC `f:real^N->real^N` THEN STRIP_TAC THEN
1926   ASM_REWRITE_TAC[] THEN
1927   SUBGOAL_THEN `span(IMAGE (f:real^N->real^N) s) = span t` MP_TAC THENL
1928    [MATCH_MP_TAC DIM_EQ_SPAN THEN ASM_REWRITE_TAC[] THEN
1929     W(MP_TAC o PART_MATCH (lhs o rand) DIM_INJECTIVE_LINEAR_IMAGE o
1930       rand o snd) THEN
1931     ASM_MESON_TAC[LE_REFL; orthogonal_transformation;
1932                   ORTHOGONAL_TRANSFORMATION_INJECTIVE];
1933     ASM_SIMP_TAC[SPAN_LINEAR_IMAGE; ORTHOGONAL_TRANSFORMATION_LINEAR] THEN
1934     ASM_SIMP_TAC[SPAN_OF_SUBSPACE]]);;
1935
1936 (* ------------------------------------------------------------------------- *)
1937 (* Rotation, reflection, rotoinversion.                                      *)
1938 (* ------------------------------------------------------------------------- *)
1939
1940 let rotation_matrix = new_definition
1941  `rotation_matrix Q <=> orthogonal_matrix Q /\ det(Q) = &1`;;
1942
1943 let rotoinversion_matrix = new_definition
1944  `rotoinversion_matrix Q <=> orthogonal_matrix Q /\ det(Q) = -- &1`;;
1945
1946 let ORTHOGONAL_ROTATION_OR_ROTOINVERSION = prove
1947  (`!Q. orthogonal_matrix Q <=> rotation_matrix Q \/ rotoinversion_matrix Q`,
1948   MESON_TAC[rotation_matrix; rotoinversion_matrix; DET_ORTHOGONAL_MATRIX]);;
1949
1950 let ROTATION_MATRIX_2 = prove
1951  (`!A:real^2^2. rotation_matrix A <=>
1952                 A$1$1 pow 2 + A$2$1 pow 2 = &1 /\
1953                 A$1$1 = A$2$2 /\ A$1$2 = --(A$2$1)`,
1954   REWRITE_TAC[rotation_matrix; ORTHOGONAL_MATRIX_2; DET_2] THEN
1955   CONV_TAC REAL_RING);;
1956
1957 (* ------------------------------------------------------------------------- *)
1958 (* Slightly stronger results giving rotation, but only in >= 2 dimensions.   *)
1959 (* ------------------------------------------------------------------------- *)
1960
1961 let ROTATION_MATRIX_EXISTS_BASIS = prove
1962  (`!a:real^N.
1963         2 <= dimindex(:N) /\ norm(a) = &1
1964         ==> ?A. rotation_matrix A /\ A**(basis 1) = a`,
1965   REPEAT STRIP_TAC THEN
1966   FIRST_X_ASSUM(X_CHOOSE_THEN `A:real^N^N` STRIP_ASSUME_TAC o
1967    MATCH_MP ORTHOGONAL_MATRIX_EXISTS_BASIS) THEN
1968   FIRST_ASSUM(DISJ_CASES_TAC o GEN_REWRITE_RULE I
1969    [ORTHOGONAL_ROTATION_OR_ROTOINVERSION])
1970   THENL [ASM_MESON_TAC[]; ALL_TAC] THEN
1971   EXISTS_TAC `transp(lambda i. if i = dimindex(:N) then -- &1 % transp A$i
1972                                else (transp A:real^N^N)$i):real^N^N` THEN
1973   REWRITE_TAC[rotation_matrix; DET_TRANSP] THEN REPEAT CONJ_TAC THENL
1974    [REWRITE_TAC[ORTHOGONAL_MATRIX_TRANSP];
1975     SIMP_TAC[DET_ROW_MUL; DIMINDEX_GE_1; LE_REFL] THEN
1976     MATCH_MP_TAC(REAL_ARITH `x = -- &1 ==> -- &1 * x = &1`) THEN
1977     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [rotoinversion_matrix]) THEN
1978     DISCH_THEN(SUBST1_TAC o SYM o CONJUNCT2) THEN
1979     GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN
1980     AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN MESON_TAC[];
1981     FIRST_X_ASSUM(SUBST1_TAC o SYM) THEN
1982     SIMP_TAC[matrix_vector_mul; LAMBDA_BETA; CART_EQ; transp;
1983              BASIS_COMPONENT] THEN
1984     ONCE_REWRITE_TAC[REAL_ARITH
1985       `x * (if p then &1 else &0) = if p then x else &0`] THEN
1986     ASM_SIMP_TAC[ARITH_RULE `2 <= n ==> ~(1 = n)`; LAMBDA_BETA]] THEN
1987   FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I
1988    [GSYM ORTHOGONAL_MATRIX_TRANSP]) THEN
1989   SPEC_TAC(`transp(A:real^N^N)`,`B:real^N^N`) THEN GEN_TAC THEN
1990   SUBGOAL_THEN `!i. 1 <= i /\ i <= dimindex(:N)
1991                     ==> row i ((lambda i. if i = dimindex(:N) then -- &1 % B$i
1992                                 else (B:real^N^N)$i):real^N^N) =
1993                         if i = dimindex(:N) then --(row i B) else row i B`
1994   ASSUME_TAC THENL
1995    [SIMP_TAC[row; LAMBDA_BETA; LAMBDA_ETA; VECTOR_MUL_LID; VECTOR_MUL_LNEG];
1996     ASM_SIMP_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS] THEN
1997     ASM_MESON_TAC[ORTHOGONAL_LNEG; ORTHOGONAL_RNEG; NORM_NEG]]);;
1998
1999 let ROTATION_EXISTS_1 = prove
2000  (`!a b:real^N.
2001         2 <= dimindex(:N) /\ norm(a) = &1 /\ norm(b) = &1
2002         ==> ?f. orthogonal_transformation f /\ det(matrix f) = &1 /\ f a = b`,
2003   REPEAT STRIP_TAC THEN
2004   MP_TAC(ISPEC `b:real^N` ROTATION_MATRIX_EXISTS_BASIS) THEN
2005   MP_TAC(ISPEC `a:real^N` ROTATION_MATRIX_EXISTS_BASIS) THEN
2006   ASM_REWRITE_TAC[rotation_matrix] THEN
2007   DISCH_THEN(X_CHOOSE_THEN `A:real^N^N`
2008    (CONJUNCTS_THEN2 STRIP_ASSUME_TAC (ASSUME_TAC o SYM))) THEN
2009   DISCH_THEN(X_CHOOSE_THEN `B:real^N^N`
2010    (CONJUNCTS_THEN2 STRIP_ASSUME_TAC (ASSUME_TAC o SYM))) THEN
2011   EXISTS_TAC `\x:real^N. ((B:real^N^N) ** transp(A:real^N^N)) ** x` THEN
2012   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX; MATRIX_VECTOR_MUL_LINEAR;
2013               MATRIX_OF_MATRIX_VECTOR_MUL; DET_MUL; DET_TRANSP] THEN
2014   ASM_SIMP_TAC[ORTHOGONAL_MATRIX_MUL; ORTHOGONAL_MATRIX_TRANSP] THEN
2015   REWRITE_TAC[GSYM MATRIX_VECTOR_MUL_ASSOC; REAL_MUL_LID] THEN AP_TERM_TAC THEN
2016   RULE_ASSUM_TAC(REWRITE_RULE[ORTHOGONAL_MATRIX]) THEN
2017   ASM_REWRITE_TAC[MATRIX_VECTOR_MUL_ASSOC; MATRIX_VECTOR_MUL_LID]);;
2018
2019 let ROTATION_EXISTS = prove
2020  (`!a b:real^N.
2021         2 <= dimindex(:N) /\ norm(a) = norm(b)
2022         ==> ?f. orthogonal_transformation f /\ det(matrix f) = &1 /\ f a = b`,
2023   REPEAT GEN_TAC THEN ASM_CASES_TAC `b:real^N = vec 0` THEN
2024   ASM_SIMP_TAC[NORM_0; NORM_EQ_0] THENL
2025    [MESON_TAC[ORTHOGONAL_TRANSFORMATION_ID; MATRIX_ID; DET_I]; ALL_TAC] THEN
2026   ASM_CASES_TAC `a:real^N = vec 0` THENL
2027    [ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_ID; MATRIX_ID; DET_I; NORM_0;
2028                   NORM_EQ_0]; ALL_TAC] THEN
2029   DISCH_TAC THEN
2030   MP_TAC(ISPECL [`inv(norm a) % a:real^N`; `inv(norm b) % b:real^N`]
2031                 ROTATION_EXISTS_1) THEN
2032   REWRITE_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
2033   ASM_SIMP_TAC[NORM_EQ_0; REAL_MUL_LINV] THEN
2034   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
2035   DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
2036   ASM_REWRITE_TAC[] THEN
2037   FIRST_ASSUM(ASSUME_TAC o MATCH_MP LINEAR_CMUL o
2038               MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
2039   ASM_REWRITE_TAC[VECTOR_ARITH
2040    `a % x:real^N = a % y <=> a % (x - y) = vec 0`] THEN
2041   ASM_REWRITE_TAC[VECTOR_MUL_EQ_0; REAL_INV_EQ_0; NORM_EQ_0; VECTOR_SUB_EQ]);;
2042
2043 let ROTATION_RIGHTWARD_LINE = prove
2044  (`!a:real^N k.
2045         1 <= k /\ k <= dimindex(:N)
2046         ==> ?b f. orthogonal_transformation f /\
2047                   (2 <= dimindex(:N) ==> det(matrix f) = &1) /\
2048                   f(b % basis k) = a /\
2049                   &0 <= b`,
2050   REPEAT STRIP_TAC THEN EXISTS_TAC `norm(a:real^N)` THEN
2051   ASM_SIMP_TAC[VECTOR_MUL_COMPONENT; BASIS_COMPONENT; LE_REFL; DIMINDEX_GE_1;
2052                REAL_MUL_RID; NORM_POS_LE; LT_IMP_LE; LTE_ANTISYM] THEN
2053   REWRITE_TAC[ARITH_RULE `2 <= n <=> 1 <= n /\ ~(n = 1)`; DIMINDEX_GE_1] THEN
2054   ASM_CASES_TAC `dimindex(:N) = 1` THEN ASM_REWRITE_TAC[] THENL
2055    [MATCH_MP_TAC ORTHOGONAL_TRANSFORMATION_EXISTS;
2056     MATCH_MP_TAC ROTATION_EXISTS] THEN
2057    ASM_SIMP_TAC[NORM_MUL; NORM_BASIS; LE_REFL; DIMINDEX_GE_1] THEN
2058    REWRITE_TAC[REAL_ABS_NORM; REAL_MUL_RID] THEN
2059    MATCH_MP_TAC(ARITH_RULE `~(n = 1) /\ 1 <= n ==> 2 <= n`) THEN
2060    ASM_REWRITE_TAC[DIMINDEX_GE_1]);;
2061
2062 (* ------------------------------------------------------------------------- *)
2063 (* In 3 dimensions, a rotation is indeed about an "axis".                    *)
2064 (* ------------------------------------------------------------------------- *)
2065
2066 let EULER_ROTATION_THEOREM = prove
2067  (`!A:real^3^3. rotation_matrix A ==> ?v:real^3. ~(v = vec 0) /\ A ** v = v`,
2068   REPEAT STRIP_TAC THEN
2069   MP_TAC(ISPEC `A - mat 1:real^3^3` HOMOGENEOUS_LINEAR_EQUATIONS_DET) THEN
2070   REWRITE_TAC[MATRIX_VECTOR_MUL_SUB_RDISTRIB;
2071               VECTOR_SUB_EQ; MATRIX_VECTOR_MUL_LID] THEN
2072   DISCH_THEN SUBST1_TAC THEN POP_ASSUM MP_TAC THEN
2073   REWRITE_TAC[rotation_matrix; orthogonal_matrix; DET_3] THEN
2074   SIMP_TAC[CART_EQ; FORALL_3; MAT_COMPONENT; DIMINDEX_3; LAMBDA_BETA; ARITH;
2075            MATRIX_SUB_COMPONENT; MAT_COMPONENT; SUM_3;
2076            matrix_mul; transp; matrix_vector_mul] THEN
2077   CONV_TAC REAL_RING);;
2078
2079 let EULER_ROTOINVERSION_THEOREM = prove
2080  (`!A:real^3^3.
2081      rotoinversion_matrix A ==> ?v:real^3. ~(v = vec 0) /\ A ** v = --v`,
2082   REPEAT STRIP_TAC THEN
2083   REWRITE_TAC[VECTOR_ARITH `a:real^N = --v <=> a + v = vec 0`] THEN
2084   MP_TAC(ISPEC `A + mat 1:real^3^3` HOMOGENEOUS_LINEAR_EQUATIONS_DET) THEN
2085   REWRITE_TAC[MATRIX_VECTOR_MUL_ADD_RDISTRIB; MATRIX_VECTOR_MUL_LID] THEN
2086   DISCH_THEN SUBST1_TAC THEN POP_ASSUM MP_TAC THEN
2087   REWRITE_TAC[rotoinversion_matrix; orthogonal_matrix; DET_3] THEN
2088   SIMP_TAC[CART_EQ; FORALL_3; MAT_COMPONENT; DIMINDEX_3; LAMBDA_BETA; ARITH;
2089            MATRIX_ADD_COMPONENT; MAT_COMPONENT; SUM_3;
2090            matrix_mul; transp; matrix_vector_mul] THEN
2091   CONV_TAC REAL_RING);;
2092
2093 (* ------------------------------------------------------------------------- *)
2094 (* We can always rotate so that a hyperplane is "horizontal".                *)
2095 (* ------------------------------------------------------------------------- *)
2096
2097 let ROTATION_LOWDIM_HORIZONTAL = prove
2098  (`!s:real^N->bool.
2099         dim s < dimindex(:N)
2100         ==> ?f. orthogonal_transformation f /\ det(matrix f) = &1 /\
2101                (IMAGE f s) SUBSET {z | z$(dimindex(:N)) = &0}`,
2102   GEN_TAC THEN ASM_CASES_TAC `dim(s:real^N->bool) = 0` THENL
2103    [RULE_ASSUM_TAC(REWRITE_RULE[DIM_EQ_0]) THEN DISCH_TAC THEN
2104     EXISTS_TAC `\x:real^N. x` THEN
2105     REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_ID; MATRIX_ID; DET_I] THEN
2106     FIRST_X_ASSUM(MATCH_MP_TAC o MATCH_MP (SET_RULE
2107      `s SUBSET {a} ==> a IN t ==> IMAGE (\x. x) s SUBSET t`)) THEN
2108     SIMP_TAC[IN_ELIM_THM; VEC_COMPONENT; LE_REFL; DIMINDEX_GE_1];
2109     DISCH_TAC] THEN
2110   SUBGOAL_THEN `2 <= dimindex(:N)` ASSUME_TAC THENL
2111    [ASM_ARITH_TAC; ALL_TAC] THEN
2112   FIRST_X_ASSUM(X_CHOOSE_THEN `a:real^N` STRIP_ASSUME_TAC o MATCH_MP
2113     LOWDIM_SUBSET_HYPERPLANE) THEN
2114   MP_TAC(ISPECL [`a:real^N`; `norm(a:real^N) % basis(dimindex(:N)):real^N`]
2115         ROTATION_EXISTS) THEN
2116   ASM_SIMP_TAC[NORM_MUL; NORM_BASIS; LE_REFL; DIMINDEX_GE_1] THEN
2117   REWRITE_TAC[REAL_ABS_NORM; REAL_MUL_RID] THEN
2118   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
2119   STRIP_TAC THEN ASM_REWRITE_TAC[FORALL_IN_IMAGE; SUBSET; IN_ELIM_THM] THEN
2120   X_GEN_TAC `x:real^N` THEN DISCH_TAC THEN
2121   SUBGOAL_THEN `(f:real^N->real^N) x dot (f a) = &0` MP_TAC THENL
2122    [RULE_ASSUM_TAC(REWRITE_RULE[orthogonal_transformation]) THEN
2123     ASM_REWRITE_TAC[] THEN ONCE_REWRITE_TAC[DOT_SYM] THEN
2124     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [SUBSET]) THEN
2125     ASM_SIMP_TAC[SPAN_SUPERSET; IN_ELIM_THM];
2126     ASM_SIMP_TAC[DOT_BASIS; LE_REFL; DIMINDEX_GE_1; DOT_RMUL] THEN
2127     ASM_REWRITE_TAC[REAL_ENTIRE; NORM_EQ_0]]);;
2128
2129 let ORTHOGONAL_TRANSFORMATION_LOWDIM_HORIZONTAL = prove
2130  (`!s:real^N->bool.
2131         dim s < dimindex(:N)
2132         ==> ?f. orthogonal_transformation f /\
2133                (IMAGE f s) SUBSET {z | z$(dimindex(:N)) = &0}`,
2134   GEN_TAC THEN DISCH_THEN(MP_TAC o MATCH_MP ROTATION_LOWDIM_HORIZONTAL) THEN
2135   MATCH_MP_TAC MONO_EXISTS THEN SIMP_TAC[]);;
2136
2137 let ORTHOGONAL_TRANSFORMATION_BETWEEN_ORTHOGONAL_SETS = prove
2138  (`!v:num->real^N w k.
2139         pairwise (\i j. orthogonal (v i) (v j)) k /\
2140         pairwise (\i j. orthogonal (w i) (w j)) k /\
2141         (!i. i IN k ==> norm(v i) = norm(w i))
2142         ==> ?f. orthogonal_transformation f /\
2143                 (!i. i IN k ==> f(v i) = w i)`,
2144   let lemma1 = prove
2145    (`!v:num->real^N n.
2146           pairwise (\i j. orthogonal (v i) (v j)) (1..n) /\
2147           (!i. 1 <= i /\ i <= n ==> norm(v i) = &1)
2148           ==> ?f. orthogonal_transformation f /\
2149                   (!i. 1 <= i /\ i <= n ==> f(basis i) = v i)`,
2150     REWRITE_TAC[pairwise; IN_NUMSEG; GSYM CONJ_ASSOC] THEN
2151     REPEAT STRIP_TAC THEN
2152     SUBGOAL_THEN `pairwise orthogonal (IMAGE (v:num->real^N) (1..n))`
2153     ASSUME_TAC THENL
2154      [REWRITE_TAC[PAIRWISE_IMAGE] THEN ASM_SIMP_TAC[pairwise; IN_NUMSEG];
2155       ALL_TAC] THEN
2156     FIRST_ASSUM(MP_TAC o MATCH_MP (REWRITE_RULE[IMP_CONJ]
2157           PAIRWISE_ORTHOGONAL_INDEPENDENT)) THEN
2158     REWRITE_TAC[SET_RULE
2159      `~(a IN IMAGE f s) <=> !x. x IN s ==> ~(f x = a)`] THEN
2160     ANTS_TAC THENL
2161      [REWRITE_TAC[IN_NUMSEG] THEN
2162       ASM_MESON_TAC[NORM_0; REAL_ARITH `~(&1 = &0)`];
2163       DISCH_THEN(MP_TAC o CONJUNCT2 o MATCH_MP INDEPENDENT_BOUND)] THEN
2164     SUBGOAL_THEN
2165      `!i j. 1 <= i /\ i <= n /\ 1 <= j /\ j <= n /\ ~(i = j)
2166             ==> ~(v i:real^N = v j)`
2167     ASSUME_TAC THENL
2168      [ASM_MESON_TAC[ORTHOGONAL_REFL; NORM_0; REAL_ARITH `~(&1 = &0)`];
2169       ALL_TAC] THEN
2170     SUBGOAL_THEN `CARD(IMAGE (v:num->real^N) (1..n)) = n` ASSUME_TAC THENL
2171      [W(MP_TAC o PART_MATCH (lhs o rand) CARD_IMAGE_INJ o lhs o snd) THEN
2172       ASM_REWRITE_TAC[CARD_NUMSEG_1; IN_NUMSEG; FINITE_NUMSEG] THEN
2173       ASM_MESON_TAC[];
2174       ASM_REWRITE_TAC[] THEN DISCH_TAC] THEN
2175     SUBGOAL_THEN
2176      `?w:num->real^N.
2177           pairwise (\i j. orthogonal (w i) (w j)) (1..dimindex(:N)) /\
2178           (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(w i) = &1) /\
2179           (!i. 1 <= i /\ i <= n ==> w i = v i)`
2180     STRIP_ASSUME_TAC THENL
2181      [ALL_TAC;
2182       EXISTS_TAC
2183        `(\x. vsum(1..dimindex(:N)) (\i. x$i % w i)):real^N->real^N` THEN
2184       SIMP_TAC[BASIS_COMPONENT; IN_NUMSEG; COND_RATOR; COND_RAND] THEN
2185       REWRITE_TAC[VECTOR_MUL_LID; VECTOR_MUL_LZERO; VSUM_DELTA] THEN
2186       ASM_SIMP_TAC[IN_NUMSEG] THEN CONJ_TAC THENL
2187        [ALL_TAC; ASM_MESON_TAC[LE_TRANS]] THEN
2188       REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX] THEN
2189       CONJ_TAC THENL
2190        [MATCH_MP_TAC LINEAR_COMPOSE_VSUM THEN
2191         REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG] THEN
2192         REWRITE_TAC[linear; VECTOR_ADD_COMPONENT; VECTOR_MUL_COMPONENT] THEN
2193         REPEAT STRIP_TAC THEN VECTOR_ARITH_TAC;
2194         ALL_TAC] THEN
2195       REWRITE_TAC[matrix; column; ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS] THEN
2196       SIMP_TAC[LAMBDA_BETA; LAMBDA_ETA; BASIS_COMPONENT; IN_NUMSEG] THEN
2197       SIMP_TAC[COND_RATOR; COND_RAND; VECTOR_MUL_LZERO; VSUM_DELTA] THEN
2198       SIMP_TAC[IN_NUMSEG; orthogonal; dot; LAMBDA_BETA; NORM_EQ_SQUARE] THEN
2199       REWRITE_TAC[VECTOR_MUL_LID; GSYM dot; GSYM NORM_EQ_SQUARE] THEN
2200       RULE_ASSUM_TAC(REWRITE_RULE[pairwise; IN_NUMSEG; orthogonal]) THEN
2201       ASM_SIMP_TAC[]] THEN
2202     FIRST_ASSUM(MP_TAC o SPEC `(:real^N)` o MATCH_MP
2203      (REWRITE_RULE[IMP_CONJ] ORTHONORMAL_EXTENSION)) THEN
2204     ASM_REWRITE_TAC[FORALL_IN_IMAGE; IN_NUMSEG; UNION_UNIV; SPAN_UNIV] THEN
2205     DISCH_THEN(X_CHOOSE_THEN `t:real^N->bool` STRIP_ASSUME_TAC) THEN
2206     MP_TAC(ISPECL [`n+1..dimindex(:N)`; `t:real^N->bool`]
2207           CARD_EQ_BIJECTION) THEN
2208     ANTS_TAC THENL
2209      [REWRITE_TAC[FINITE_NUMSEG] THEN
2210       MP_TAC(ISPECL [`(:real^N)`; `IMAGE v (1..n) UNION t:real^N->bool`]
2211           BASIS_CARD_EQ_DIM) THEN
2212       ASM_REWRITE_TAC[SUBSET_UNIV] THEN ANTS_TAC THENL
2213        [MATCH_MP_TAC PAIRWISE_ORTHOGONAL_INDEPENDENT THEN
2214         ASM_REWRITE_TAC[IN_UNION; DE_MORGAN_THM; IN_NUMSEG] THEN
2215         ASM_REWRITE_TAC[FORALL_IN_IMAGE; IN_NUMSEG; SET_RULE
2216          `~(x IN s) <=> !y. y IN s ==> ~(y = x)`] THEN
2217         ASM_MESON_TAC[NORM_0; REAL_ARITH `~(&1 = &0)`];
2218         ALL_TAC] THEN
2219       ASM_SIMP_TAC[FINITE_UNION; IMP_CONJ; FINITE_IMAGE; CARD_UNION;
2220                    SET_RULE `t INTER s = {} <=> DISJOINT s t`] THEN
2221       DISCH_TAC THEN DISCH_TAC THEN REWRITE_TAC[CARD_NUMSEG; DIM_UNIV] THEN
2222       ARITH_TAC;
2223       ALL_TAC] THEN
2224     REWRITE_TAC[CONJ_ASSOC; SET_RULE
2225      `(!x. x IN s ==> f x IN t) /\ (!y. y IN t ==> ?x. x IN s /\ f x = y) <=>
2226       t = IMAGE f s`] THEN
2227     REWRITE_TAC[GSYM CONJ_ASSOC; LEFT_IMP_EXISTS_THM; IN_NUMSEG] THEN
2228     X_GEN_TAC `w:num->real^N` THEN
2229     DISCH_THEN(CONJUNCTS_THEN2 SUBST_ALL_TAC MP_TAC) THEN
2230     REWRITE_TAC[ARITH_RULE `n + 1 <= x <=> n < x`; CONJ_ASSOC] THEN
2231     ONCE_REWRITE_TAC[TAUT `p /\ q ==> r <=> p /\ ~r ==> ~q`] THEN
2232     REWRITE_TAC[GSYM CONJ_ASSOC] THEN STRIP_TAC THEN
2233     REWRITE_TAC[TAUT `p /\ ~r ==> ~q <=> p /\ q ==> r`] THEN
2234     EXISTS_TAC `\i. if i <= n then (v:num->real^N) i else w i` THEN
2235     SIMP_TAC[] THEN
2236     RULE_ASSUM_TAC(REWRITE_RULE[FORALL_IN_IMAGE; IN_NUMSEG]) THEN
2237     CONJ_TAC THENL
2238      [ALL_TAC; ASM_MESON_TAC[ARITH_RULE `~(i <= n) ==> n + 1 <= i`]] THEN
2239     REWRITE_TAC[pairwise] THEN MATCH_MP_TAC WLOG_LT THEN REWRITE_TAC[] THEN
2240     CONJ_TAC THENL [MESON_TAC[ORTHOGONAL_SYM]; ALL_TAC] THEN
2241     MAP_EVERY X_GEN_TAC [`i:num`; `j:num`] THEN DISCH_TAC THEN
2242     ASM_CASES_TAC `j:num <= n` THEN ASM_REWRITE_TAC[IN_NUMSEG] THENL
2243      [COND_CASES_TAC THEN ASM_SIMP_TAC[] THEN ASM_ARITH_TAC; ALL_TAC] THEN
2244     ASM_CASES_TAC `i:num <= n` THEN ASM_REWRITE_TAC[] THEN STRIP_TAC THEN
2245     UNDISCH_TAC
2246      `pairwise orthogonal
2247         (IMAGE (v:num->real^N) (1..n) UNION IMAGE w (n+1..dimindex (:N)))` THEN
2248     REWRITE_TAC[pairwise] THEN ONCE_REWRITE_TAC[SWAP_FORALL_THM] THEN
2249     DISCH_THEN(MP_TAC o SPEC `(w:num->real^N) j`) THENL
2250      [DISCH_THEN(MP_TAC o SPEC `(v:num->real^N) i`);
2251       DISCH_THEN(MP_TAC o SPEC `(w:num->real^N) i`)] THEN
2252     ASM_REWRITE_TAC[IN_UNION; IN_IMAGE; IN_NUMSEG] THEN
2253     DISCH_THEN MATCH_MP_TAC THENL
2254      [CONJ_TAC THENL [ASM_MESON_TAC[]; ALL_TAC] THEN
2255       CONJ_TAC THENL
2256        [ASM_MESON_TAC[ARITH_RULE `~(x <= n) ==> n + 1 <= x`]; ALL_TAC];
2257       ASM_MESON_TAC[ARITH_RULE `~(x <= n) ==> n + 1 <= x /\ n < x`]] THEN
2258     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [DISJOINT]) THEN
2259     REWRITE_TAC[SET_RULE `IMAGE w t INTER IMAGE v s = {} <=>
2260       !i j. i IN s /\ j IN t ==> ~(v i = w j)`] THEN
2261     DISCH_THEN MATCH_MP_TAC THEN ASM_REWRITE_TAC[IN_NUMSEG] THEN
2262     ASM_ARITH_TAC) in
2263   let lemma2 = prove
2264    (`!v:num->real^N w k.
2265           pairwise (\i j. orthogonal (v i) (v j)) k /\
2266           pairwise (\i j. orthogonal (w i) (w j)) k /\
2267           (!i. i IN k ==> norm(v i) = norm(w i)) /\
2268           (!i. i IN k ==> ~(v i = vec 0) /\ ~(w i = vec 0))
2269           ==> ?f. orthogonal_transformation f /\
2270                   (!i. i IN k ==> f(v i) = w i)`,
2271     REPEAT STRIP_TAC THEN
2272     SUBGOAL_THEN `FINITE(k:num->bool)` MP_TAC THENL
2273      [SUBGOAL_THEN `pairwise orthogonal (IMAGE (v:num->real^N) k)`
2274       ASSUME_TAC THENL
2275        [REWRITE_TAC[PAIRWISE_IMAGE] THEN
2276         RULE_ASSUM_TAC(REWRITE_RULE[pairwise]) THEN ASM_SIMP_TAC[pairwise];
2277         ALL_TAC] THEN
2278       FIRST_ASSUM(MP_TAC o MATCH_MP (REWRITE_RULE[IMP_CONJ]
2279         PAIRWISE_ORTHOGONAL_INDEPENDENT)) THEN
2280       ANTS_TAC THENL [ASM SET_TAC[]; ALL_TAC] THEN
2281       DISCH_THEN(MP_TAC o MATCH_MP INDEPENDENT_IMP_FINITE) THEN
2282       MATCH_MP_TAC EQ_IMP THEN MATCH_MP_TAC FINITE_IMAGE_INJ_EQ THEN
2283       RULE_ASSUM_TAC(REWRITE_RULE[pairwise]) THEN
2284       ASM_MESON_TAC[ORTHOGONAL_REFL];
2285       ALL_TAC] THEN
2286     DISCH_THEN(MP_TAC o GEN_REWRITE_RULE I [FINITE_INDEX_NUMSEG]) THEN
2287     ONCE_REWRITE_TAC[TAUT `p /\ q /\ r ==> s <=> p /\ q /\ ~s ==> ~r`] THEN
2288     DISCH_THEN(X_CHOOSE_THEN `n:num->num` MP_TAC) THEN
2289     REWRITE_TAC[IN_NUMSEG] THEN GEN_REWRITE_TAC I [IMP_CONJ] THEN
2290     DISCH_THEN(fun th -> DISCH_THEN SUBST_ALL_TAC THEN ASSUME_TAC th) THEN
2291     RULE_ASSUM_TAC(REWRITE_RULE
2292      [PAIRWISE_IMAGE; FORALL_IN_IMAGE; IN_NUMSEG]) THEN
2293     MP_TAC(ISPECL
2294      [`\i. inv(norm(w(n i))) % (w:num->real^N) ((n:num->num) i)`;
2295       `CARD(k:num->bool)`] lemma1) THEN
2296     MP_TAC(ISPECL
2297      [`\i. inv(norm(v(n i))) % (v:num->real^N) ((n:num->num) i)`;
2298       `CARD(k:num->bool)`] lemma1) THEN
2299     ASM_SIMP_TAC[NORM_MUL; REAL_MUL_LINV; NORM_EQ_0; REAL_ABS_INV;
2300                  REAL_ABS_NORM; pairwise; orthogonal; IN_NUMSEG] THEN
2301     RULE_ASSUM_TAC(REWRITE_RULE[pairwise; orthogonal; IN_NUMSEG]) THEN
2302     ASM_SIMP_TAC[DOT_LMUL; DOT_RMUL; REAL_ENTIRE; FORALL_IN_IMAGE] THEN
2303     DISCH_THEN(X_CHOOSE_THEN `f:real^N->real^N` STRIP_ASSUME_TAC) THEN
2304     DISCH_THEN(X_CHOOSE_THEN `g:real^N->real^N` STRIP_ASSUME_TAC) THEN
2305     MP_TAC(ISPEC `f:real^N->real^N` ORTHOGONAL_TRANSFORMATION_INVERSE) THEN
2306     ASM_REWRITE_TAC[] THEN
2307     DISCH_THEN(X_CHOOSE_THEN `f':real^N->real^N` STRIP_ASSUME_TAC) THEN
2308     EXISTS_TAC `(g:real^N->real^N) o (f':real^N->real^N)` THEN
2309     ASM_SIMP_TAC[ORTHOGONAL_TRANSFORMATION_COMPOSE; IN_NUMSEG] THEN
2310     X_GEN_TAC `i:num` THEN DISCH_TAC THEN REWRITE_TAC[o_THM] THEN
2311     MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC
2312      `(g:real^N->real^N) (norm((w:num->real^N)(n(i:num))) % basis i)` THEN
2313     CONJ_TAC THENL
2314      [AP_TERM_TAC THEN FIRST_X_ASSUM(MATCH_MP_TAC o MATCH_MP (MESON[]
2315        `(!x. f'(f x) = x) ==> f x = y ==> f' y = x`));
2316       ALL_TAC] THEN
2317     RULE_ASSUM_TAC(REWRITE_RULE[orthogonal_transformation]) THEN
2318     ASM_SIMP_TAC[LINEAR_CMUL; VECTOR_MUL_ASSOC] THEN
2319     ASM_SIMP_TAC[REAL_MUL_RINV; NORM_EQ_0; VECTOR_MUL_LID]) in
2320   REPEAT STRIP_TAC THEN MP_TAC(ISPECL
2321    [`v:num->real^N`; `w:num->real^N`;
2322     `{i | i IN k /\ ~((v:num->real^N) i = vec 0)}`] lemma2) THEN
2323   ASM_SIMP_TAC[IN_ELIM_THM; CONJ_ASSOC] THEN ANTS_TAC THENL
2324    [CONJ_TAC THENL [ALL_TAC; ASM_MESON_TAC[NORM_EQ_0]] THEN
2325     CONJ_TAC THEN MATCH_MP_TAC PAIRWISE_MONO THEN EXISTS_TAC `k:num->bool` THEN
2326     ASM_REWRITE_TAC[] THEN SET_TAC[];
2327     MATCH_MP_TAC MONO_EXISTS THEN SIMP_TAC[orthogonal_transformation] THEN
2328     GEN_TAC THEN STRIP_TAC THEN X_GEN_TAC `i:num` THEN DISCH_TAC THEN
2329     ASM_CASES_TAC `(v:num->real^N) i = vec 0` THEN ASM_SIMP_TAC[] THEN
2330     ASM_MESON_TAC[LINEAR_0; NORM_EQ_0]]);;
2331
2332 (* ------------------------------------------------------------------------- *)
2333 (* Reflection of a vector about 0 along a line.                              *)
2334 (* ------------------------------------------------------------------------- *)
2335
2336 let reflect_along = new_definition
2337  `reflect_along v (x:real^N) = x - (&2 * (x dot v) / (v dot v)) % v`;;
2338
2339 let REFLECT_ALONG_ADD = prove
2340  (`!v x y:real^N.
2341       reflect_along v (x + y) = reflect_along v x + reflect_along v y`,
2342   REPEAT GEN_TAC THEN
2343   REWRITE_TAC[reflect_along; VECTOR_ARITH
2344    `x - a % v + y - b % v:real^N = (x + y) - (a + b) % v`] THEN
2345   AP_TERM_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN
2346   REWRITE_TAC[DOT_LADD] THEN REAL_ARITH_TAC);;
2347
2348 let REFLECT_ALONG_MUL = prove
2349  (`!v a x:real^N. reflect_along v (a % x) = a % reflect_along v x`,
2350   REWRITE_TAC[reflect_along; DOT_LMUL; REAL_ARITH
2351    `&2 * (a * x) / y = a * &2 * x / y`] THEN
2352   REWRITE_TAC[VECTOR_SUB_LDISTRIB; VECTOR_MUL_ASSOC]);;
2353
2354 let LINEAR_REFLECT_ALONG = prove
2355  (`!v:real^N. linear(reflect_along v)`,
2356   REWRITE_TAC[linear; REFLECT_ALONG_ADD; REFLECT_ALONG_MUL]);;
2357
2358 let REFLECT_ALONG_0 = prove
2359  (`!v:real^N. reflect_along v (vec 0) = vec 0`,
2360   REWRITE_TAC[MATCH_MP LINEAR_0 (SPEC_ALL LINEAR_REFLECT_ALONG)]);;
2361
2362 let REFLECT_ALONG_REFL = prove
2363  (`!v:real^N. reflect_along v v = --v`,
2364   GEN_TAC THEN ASM_CASES_TAC `v:real^N = vec 0` THEN
2365   ASM_REWRITE_TAC[VECTOR_NEG_0; REFLECT_ALONG_0] THEN
2366   REWRITE_TAC[reflect_along] THEN
2367   ASM_SIMP_TAC[REAL_DIV_REFL; DOT_EQ_0] THEN VECTOR_ARITH_TAC);;
2368
2369 let REFLECT_ALONG_INVOLUTION = prove
2370  (`!v x:real^N. reflect_along v (reflect_along v x) = x`,
2371   REWRITE_TAC[reflect_along; DOT_LSUB; VECTOR_MUL_EQ_0; VECTOR_ARITH
2372    `x - a % v - b % v:real^N = x <=> (a + b) % v = vec 0`] THEN
2373   REWRITE_TAC[DOT_LMUL; GSYM DOT_EQ_0] THEN CONV_TAC REAL_FIELD);;
2374
2375 let REFLECT_ALONG_EQ_0 = prove
2376  (`!v x:real^N. reflect_along v x = vec 0 <=> x = vec 0`,
2377   MESON_TAC[REFLECT_ALONG_0; REFLECT_ALONG_INVOLUTION]);;
2378
2379 let ORTHGOONAL_TRANSFORMATION_REFLECT_ALONG = prove
2380  (`!v:real^N. orthogonal_transformation(reflect_along v)`,
2381   GEN_TAC THEN ASM_CASES_TAC `v:real^N = vec 0` THENL
2382    [GEN_REWRITE_TAC RAND_CONV [GSYM ETA_AX] THEN
2383     ASM_REWRITE_TAC[reflect_along; VECTOR_MUL_RZERO; VECTOR_SUB_RZERO;
2384                     ORTHOGONAL_TRANSFORMATION_ID];
2385     REWRITE_TAC[ORTHOGONAL_TRANSFORMATION] THEN
2386     REWRITE_TAC[LINEAR_REFLECT_ALONG; NORM_EQ] THEN
2387     REWRITE_TAC[reflect_along; VECTOR_ARITH
2388       `(a - b:real^N) dot (a - b) = (a dot a + b dot b) - &2 * a dot b`] THEN
2389     REWRITE_TAC[DOT_LMUL; DOT_RMUL] THEN X_GEN_TAC `w:real^N` THEN
2390     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV [GSYM DOT_EQ_0]) THEN
2391     CONV_TAC REAL_FIELD]);;
2392
2393 let REFLECT_ALONG_EQ_SELF = prove
2394  (`!v x:real^N. reflect_along v x = x <=> orthogonal v x`,
2395   REPEAT GEN_TAC THEN REWRITE_TAC[reflect_along; orthogonal] THEN
2396   REWRITE_TAC[VECTOR_ARITH `x - a:real^N = x <=> a = vec 0`] THEN
2397   REWRITE_TAC[VECTOR_MUL_EQ_0] THEN
2398   ASM_CASES_TAC `v:real^N = vec 0` THEN ASM_SIMP_TAC[DOT_LZERO; DOT_SYM] THEN
2399   FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV [GSYM DOT_EQ_0]) THEN
2400   CONV_TAC REAL_FIELD);;
2401
2402 let REFLECT_ALONG_ZERO = prove
2403  (`!x:real^N. reflect_along (vec 0) = I`,
2404   REWRITE_TAC[FUN_EQ_THM; I_THM; REFLECT_ALONG_EQ_SELF; ORTHOGONAL_0]);;
2405
2406 let REFLECT_ALONG_LINEAR_IMAGE = prove
2407  (`!f:real^M->real^N v x.
2408         linear f /\ (!x. norm(f x) = norm x)
2409         ==> reflect_along (f v) (f x) = f(reflect_along v x)`,
2410   REWRITE_TAC[reflect_along] THEN
2411   SIMP_TAC[PRESERVES_NORM_PRESERVES_DOT; LINEAR_SUB; LINEAR_CMUL]);;
2412
2413 add_linear_invariants [REFLECT_ALONG_LINEAR_IMAGE];;
2414
2415 let REFLECT_ALONG_SCALE = prove
2416  (`!c v x:real^N. ~(c = &0) ==> reflect_along (c % v) x = reflect_along v x`,
2417   REPEAT STRIP_TAC THEN
2418   ASM_CASES_TAC `v:real^N = vec 0` THEN
2419   ASM_REWRITE_TAC[VECTOR_MUL_RZERO; REFLECT_ALONG_ZERO] THEN
2420   REWRITE_TAC[reflect_along; VECTOR_MUL_ASSOC] THEN
2421   AP_TERM_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN
2422   REWRITE_TAC[DOT_RMUL] THEN REWRITE_TAC[DOT_LMUL] THEN
2423   FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE RAND_CONV [GSYM DOT_EQ_0]) THEN
2424   POP_ASSUM MP_TAC THEN CONV_TAC REAL_FIELD);;
2425
2426 let REFLECT_ALONG_1D = prove
2427  (`!v x:real^N.
2428         dimindex(:N) = 1 ==> reflect_along v x = if v = vec 0 then x else --x`,
2429   REPEAT STRIP_TAC THEN
2430   ASM_SIMP_TAC[reflect_along; dot; SUM_1; CART_EQ; FORALL_1] THEN
2431   REWRITE_TAC[VEC_COMPONENT; COND_RATOR; COND_RAND] THEN
2432   SIMP_TAC[VECTOR_NEG_COMPONENT; VECTOR_MUL_COMPONENT;
2433            VECTOR_SUB_COMPONENT; REAL_MUL_RZERO] THEN
2434   CONV_TAC REAL_FIELD);;
2435
2436 let REFLECT_ALONG_BASIS = prove
2437  (`!x:real^N k.
2438         1 <= k /\ k <= dimindex(:N)
2439         ==> reflect_along (basis k) x = x - (&2 * x$k) % basis k`,
2440   SIMP_TAC[reflect_along; DOT_BASIS; BASIS_COMPONENT; REAL_DIV_1]);;
2441
2442 let MATRIX_REFLECT_ALONG_BASIS = prove
2443  (`!k. 1 <= k /\ k <= dimindex(:N)
2444        ==> matrix(reflect_along (basis k)):real^N^N =
2445            lambda i j. if i = k /\ j = k then --(&1)
2446                        else if i = j then &1
2447                        else &0`,
2448   SIMP_TAC[CART_EQ; LAMBDA_BETA; matrix; REFLECT_ALONG_BASIS;
2449            VECTOR_SUB_COMPONENT; BASIS_COMPONENT; VECTOR_MUL_COMPONENT] THEN
2450   X_GEN_TAC `k:num` THEN STRIP_TAC THEN
2451   X_GEN_TAC `i:num` THEN STRIP_TAC THEN
2452   X_GEN_TAC `j:num` THEN STRIP_TAC THEN
2453   ASM_CASES_TAC `i:num = j` THEN ASM_REWRITE_TAC[] THEN
2454   REPEAT(COND_CASES_TAC THEN ASM_REWRITE_TAC[]) THEN ASM_REAL_ARITH_TAC);;
2455
2456 let ROTOINVERSION_MATRIX_REFLECT_ALONG = prove
2457  (`!v:real^N. ~(v = vec 0) ==> rotoinversion_matrix(matrix(reflect_along v))`,
2458   REPEAT STRIP_TAC THEN REWRITE_TAC[rotoinversion_matrix] THEN
2459   CONJ_TAC THENL
2460    [ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX;
2461                 ORTHGOONAL_TRANSFORMATION_REFLECT_ALONG];
2462     ALL_TAC] THEN
2463   ABBREV_TAC `w:real^N = inv(norm v) % v` THEN
2464   SUBGOAL_THEN `reflect_along (v:real^N) = reflect_along w` SUBST1_TAC THENL
2465    [EXPAND_TAC "w" THEN REWRITE_TAC[FUN_EQ_THM] THEN
2466     ASM_SIMP_TAC[REFLECT_ALONG_SCALE; REAL_INV_EQ_0; NORM_EQ_0];
2467     SUBGOAL_THEN `norm(w:real^N) = &1` MP_TAC THENL
2468      [EXPAND_TAC "w" THEN SIMP_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
2469       MATCH_MP_TAC REAL_MUL_LINV THEN ASM_REWRITE_TAC[NORM_EQ_0];
2470       POP_ASSUM_LIST(K ALL_TAC) THEN SPEC_TAC(`w:real^N`,`v:real^N`)]] THEN
2471   X_GEN_TAC `v:real^N` THEN ASM_CASES_TAC `v:real^N = vec 0` THEN
2472   ASM_REWRITE_TAC[NORM_0; REAL_OF_NUM_EQ; ARITH_EQ] THEN DISCH_TAC THEN
2473   MP_TAC(ISPECL [`v:real^N`; `basis 1:real^N`]
2474         ORTHOGONAL_TRANSFORMATION_EXISTS) THEN
2475   ASM_SIMP_TAC[NORM_BASIS; DIMINDEX_GE_1; LE_REFL] THEN
2476   DISCH_THEN(X_CHOOSE_THEN `f:real^N->real^N` STRIP_ASSUME_TAC) THEN
2477   SUBGOAL_THEN
2478    `matrix(reflect_along v) =
2479     transp(matrix(f:real^N->real^N)) ** matrix(reflect_along (f v)) ** matrix f`
2480   SUBST1_TAC THENL
2481    [UNDISCH_THEN `(f:real^N->real^N) v = basis 1` (K ALL_TAC) THEN
2482     REWRITE_TAC[MATRIX_EQ; GSYM MATRIX_VECTOR_MUL_ASSOC] THEN
2483     ASM_SIMP_TAC[MATRIX_WORKS; LINEAR_REFLECT_ALONG;
2484                  ORTHOGONAL_TRANSFORMATION_LINEAR] THEN
2485     X_GEN_TAC `x:real^N` THEN MATCH_MP_TAC EQ_TRANS THEN
2486     EXISTS_TAC `(transp(matrix(f:real^N->real^N)) ** matrix f) **
2487                 (reflect_along v x:real^N)` THEN
2488     CONJ_TAC THENL
2489      [ASM_MESON_TAC[ORTHOGONAL_MATRIX; MATRIX_VECTOR_MUL_LID;
2490                 ORTHOGONAL_TRANSFORMATION_MATRIX];
2491       REWRITE_TAC[GSYM MATRIX_VECTOR_MUL_ASSOC] THEN
2492       ASM_SIMP_TAC[MATRIX_WORKS; ORTHOGONAL_TRANSFORMATION_LINEAR] THEN
2493       AP_TERM_TAC THEN CONV_TAC SYM_CONV THEN
2494       MATCH_MP_TAC REFLECT_ALONG_LINEAR_IMAGE THEN
2495       ASM_REWRITE_TAC[GSYM ORTHOGONAL_TRANSFORMATION]];
2496     ASM_REWRITE_TAC[DET_MUL; DET_TRANSP] THEN
2497     MATCH_MP_TAC(REAL_RING
2498      `(x = &1 \/ x = -- &1) /\ y = a ==> x * y * x = a`) THEN
2499     CONJ_TAC THENL
2500      [ASM_MESON_TAC[DET_ORTHOGONAL_MATRIX; ORTHOGONAL_TRANSFORMATION_MATRIX];
2501       ALL_TAC] THEN
2502     W(MP_TAC o PART_MATCH (lhs o rand) DET_UPPERTRIANGULAR o lhand o snd) THEN
2503     SIMP_TAC[MATRIX_REFLECT_ALONG_BASIS; DIMINDEX_GE_1; LE_REFL] THEN
2504     SIMP_TAC[LAMBDA_BETA; ARITH_RULE
2505      `j < i ==> ~(i = j) /\ ~(i = 1 /\ j = 1)`] THEN
2506     DISCH_THEN(K ALL_TAC) THEN
2507     SIMP_TAC[PRODUCT_CLAUSES_LEFT; DIMINDEX_GE_1] THEN
2508     MATCH_MP_TAC(REAL_RING `x = &1 ==> a * x = a`) THEN
2509     MATCH_MP_TAC PRODUCT_EQ_1 THEN
2510     REWRITE_TAC[IN_NUMSEG] THEN REPEAT STRIP_TAC THEN
2511     COND_CASES_TAC THEN ASM_REWRITE_TAC[] THEN ASM_ARITH_TAC]);;
2512
2513 let DET_MATRIX_REFLECT_ALONG = prove
2514  (`!v:real^N. det(matrix(reflect_along v)) =
2515                 if v = vec 0 then &1 else --(&1)`,
2516   GEN_TAC THEN COND_CASES_TAC THEN ASM_REWRITE_TAC[REFLECT_ALONG_ZERO] THEN
2517   REWRITE_TAC[MATRIX_I; DET_I] THEN
2518   FIRST_ASSUM(MP_TAC o MATCH_MP ROTOINVERSION_MATRIX_REFLECT_ALONG) THEN
2519   SIMP_TAC[rotoinversion_matrix]);;
2520
2521 (* ------------------------------------------------------------------------- *)
2522 (* All orthogonal transformations are a composition of reflections.          *)
2523 (* ------------------------------------------------------------------------- *)
2524
2525 let ORTHOGONAL_TRANSFORMATION_GENERATED_BY_REFLECTIONS = prove
2526  (`!f:real^N->real^N n.
2527         orthogonal_transformation f /\
2528         dimindex(:N) <= dim {x | f x = x} + n
2529         ==> ?l. LENGTH l <= n /\ ALL (\v. ~(v = vec 0)) l /\
2530                 f = ITLIST (\v h. reflect_along v o h) l I`,
2531   ONCE_REWRITE_TAC[GSYM SWAP_FORALL_THM] THEN INDUCT_TAC THENL
2532    [REWRITE_TAC[CONJUNCT1 LE; LENGTH_EQ_NIL; ADD_CLAUSES; UNWIND_THM2] THEN
2533     SIMP_TAC[DIM_SUBSET_UNIV; ARITH_RULE `a:num <= b ==> (b <= a <=> a = b)`;
2534              ITLIST; DIM_EQ_FULL; orthogonal_transformation] THEN
2535     SIMP_TAC[SPAN_OF_SUBSPACE; SUBSPACE_LINEAR_FIXED_POINTS; IMP_CONJ] THEN
2536     REWRITE_TAC[EXTENSION; IN_UNIV; IN_ELIM_THM] THEN
2537     SIMP_TAC[FUN_EQ_THM; I_THM; ALL];
2538     REPEAT STRIP_TAC THEN ASM_CASES_TAC `!x:real^N. f x = x` THENL
2539      [EXISTS_TAC `[]:(real^N) list` THEN
2540       ASM_REWRITE_TAC[ITLIST; FUN_EQ_THM; I_THM; ALL; LENGTH; LE_0];
2541       FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [NOT_FORALL_THM])] THEN
2542     DISCH_THEN(X_CHOOSE_TAC `a:real^N`) THEN
2543     ABBREV_TAC `v:real^N = inv(&2) % (f a - a)` THEN FIRST_X_ASSUM
2544       (MP_TAC o SPEC `reflect_along v o (f:real^N->real^N)`) THEN
2545     ASM_SIMP_TAC[ORTHGOONAL_TRANSFORMATION_REFLECT_ALONG;
2546                  ORTHOGONAL_TRANSFORMATION_COMPOSE] THEN
2547     ANTS_TAC THENL
2548      [FIRST_X_ASSUM(MATCH_MP_TAC o MATCH_MP (ARITH_RULE
2549        `a <= d + SUC n ==> d < d' ==> a <= d' + n`)) THEN
2550       MATCH_MP_TAC DIM_PSUBSET THEN REWRITE_TAC[PSUBSET_ALT] THEN
2551       SUBGOAL_THEN
2552        `!y:real^N. dist(y,f a) = dist(y,a) ==> reflect_along v y = y`
2553       ASSUME_TAC THENL
2554        [REWRITE_TAC[dist; NORM_EQ_SQUARE; NORM_POS_LE; NORM_POW_2] THEN
2555         REWRITE_TAC[VECTOR_ARITH
2556          `(y - b:real^N) dot (y - b) =
2557           (y dot y + b dot b) - &2 * y dot b`] THEN
2558         REWRITE_TAC[REAL_ARITH `(y + aa) - &2 * a = (y + bb) - &2 * b <=>
2559                                 a - b = inv(&2) * (aa - bb)`] THEN
2560         RULE_ASSUM_TAC(REWRITE_RULE[orthogonal_transformation]) THEN
2561         ASM_REWRITE_TAC[REAL_SUB_REFL; REAL_MUL_RZERO] THEN
2562         EXPAND_TAC "v" THEN REWRITE_TAC[GSYM DOT_RSUB; reflect_along] THEN
2563         SIMP_TAC[DOT_RMUL; real_div; REAL_MUL_LZERO; REAL_MUL_RZERO] THEN
2564         REWRITE_TAC[VECTOR_MUL_LZERO; VECTOR_SUB_RZERO];
2565         ALL_TAC] THEN
2566       CONJ_TAC THENL
2567        [MATCH_MP_TAC SPAN_MONO THEN SIMP_TAC[SUBSET; IN_ELIM_THM; o_THM] THEN
2568         ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_ISOMETRY];
2569         ALL_TAC] THEN
2570       EXISTS_TAC `a:real^N` THEN
2571       ASM_SIMP_TAC[SUBSPACE_LINEAR_FIXED_POINTS; SPAN_OF_SUBSPACE;
2572                    ORTHOGONAL_TRANSFORMATION_LINEAR; IN_ELIM_THM] THEN
2573       MATCH_MP_TAC SPAN_SUPERSET THEN REWRITE_TAC[IN_ELIM_THM; o_THM] THEN
2574       MATCH_MP_TAC EQ_TRANS THEN
2575       EXISTS_TAC `reflect_along (v:real^N) (midpoint(f a,a) + v)` THEN
2576       CONJ_TAC THENL
2577        [AP_TERM_TAC;
2578         REWRITE_TAC[REFLECT_ALONG_ADD] THEN
2579         ASM_SIMP_TAC[DIST_MIDPOINT; REFLECT_ALONG_REFL]] THEN
2580       EXPAND_TAC "v" THEN REWRITE_TAC[midpoint] THEN VECTOR_ARITH_TAC;
2581       DISCH_THEN(X_CHOOSE_THEN `l:(real^N)list` STRIP_ASSUME_TAC) THEN
2582       EXISTS_TAC `CONS (v:real^N) l` THEN
2583       ASM_REWRITE_TAC[ALL; LENGTH; LE_SUC; VECTOR_SUB_EQ; ITLIST] THEN
2584       EXPAND_TAC "v" THEN ASM_REWRITE_TAC[VECTOR_MUL_EQ_0] THEN
2585       CONV_TAC REAL_RAT_REDUCE_CONV THEN ASM_REWRITE_TAC[VECTOR_SUB_EQ] THEN
2586       FIRST_X_ASSUM(MP_TAC o AP_TERM
2587        `(o)(reflect_along (v:real^N)):(real^N->real^N)->(real^N->real^N)`) THEN
2588       REWRITE_TAC[FUN_EQ_THM; o_THM; REFLECT_ALONG_INVOLUTION]]]);;
2589
2590 (* ------------------------------------------------------------------------- *)
2591 (* Extract scaling, translation and linear invariance theorems.              *)
2592 (* For the linear case, chain through some basic consequences automatically, *)
2593 (* e.g. norm-preserving and linear implies injective.                        *)
2594 (* ------------------------------------------------------------------------- *)
2595
2596 let SCALING_THEOREMS v =
2597   let th1 = UNDISCH(snd(EQ_IMP_RULE(ISPEC v NORM_POS_LT))) in
2598   let t = rand(concl th1) in
2599   end_itlist CONJ (map (C MP th1 o SPEC t) (!scaling_theorems));;
2600
2601 let TRANSLATION_INVARIANTS x =
2602   end_itlist CONJ (mapfilter (ISPEC x) (!invariant_under_translation));;
2603
2604 let USABLE_CONCLUSION f ths th =
2605   let ith = PURE_REWRITE_RULE[RIGHT_FORALL_IMP_THM] (ISPEC f th) in
2606   let bod = concl ith in
2607   let cjs = conjuncts(fst(dest_imp bod)) in
2608   let ths = map (fun t -> find(fun th -> aconv (concl th) t) ths) cjs in
2609   GEN_ALL(MP ith (end_itlist CONJ ths));;
2610
2611 let LINEAR_INVARIANTS =
2612   let sths = (CONJUNCTS o prove)
2613    (`(!f:real^M->real^N.
2614          linear f /\ (!x. norm(f x) = norm x)
2615          ==> (!x y. f x = f y ==> x = y)) /\
2616      (!f:real^N->real^N.
2617          linear f /\ (!x. norm(f x) = norm x) ==> (!y. ?x. f x = y)) /\
2618      (!f:real^N->real^N. linear f /\ (!x y. f x = f y ==> x = y)
2619                          ==> (!y. ?x. f x = y)) /\
2620      (!f:real^N->real^N. linear f /\ (!y. ?x. f x = y)
2621                          ==> (!x y. f x = f y ==> x = y))`,
2622     CONJ_TAC THENL
2623      [ONCE_REWRITE_TAC[GSYM VECTOR_SUB_EQ] THEN
2624       SIMP_TAC[GSYM LINEAR_SUB; GSYM NORM_EQ_0];
2625       MESON_TAC[ORTHOGONAL_TRANSFORMATION_SURJECTIVE;
2626                 ORTHOGONAL_TRANSFORMATION_INJECTIVE; ORTHOGONAL_TRANSFORMATION;
2627                 LINEAR_SURJECTIVE_IFF_INJECTIVE]]) in
2628   fun f ths ->
2629     let ths' = ths @ mapfilter (USABLE_CONCLUSION f ths) sths in
2630     end_itlist CONJ
2631      (mapfilter (USABLE_CONCLUSION f ths') (!invariant_under_linear));;
2632
2633 (* ------------------------------------------------------------------------- *)
2634 (* Tactic to pick WLOG a particular point as the origin. The conversion form *)
2635 (* assumes it's the outermost universal variable; the tactic is more general *)
2636 (* and allows any free or outer universally quantified variable. The list    *)
2637 (* "avoid" is the points not to translate. There is also a tactic to help in *)
2638 (* proving new translation theorems, which uses similar machinery.           *)
2639 (* ------------------------------------------------------------------------- *)
2640
2641 let GEOM_ORIGIN_CONV,GEOM_TRANSLATE_CONV =
2642   let pth = prove
2643    (`!a:real^N. a = a + vec 0 /\
2644                 {} = IMAGE (\x. a + x) {} /\
2645                 {} = IMAGE (IMAGE (\x. a + x)) {} /\
2646                 (:real^N) = IMAGE (\x. a + x) (:real^N) /\
2647                 (:real^N->bool) = IMAGE (IMAGE (\x. a + x)) (:real^N->bool) /\
2648                 [] = MAP (\x. a + x) []`,
2649     REWRITE_TAC[IMAGE_CLAUSES; VECTOR_ADD_RID; MAP] THEN
2650     REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
2651     REWRITE_TAC[SURJECTIVE_IMAGE] THEN
2652     REWRITE_TAC[VECTOR_ARITH `a + y:real^N = x <=> y = x - a`; EXISTS_REFL])
2653   and qth = prove
2654    (`!a:real^N.
2655         ((!P. (!x. P x) <=> (!x. P (a + x))) /\
2656          (!P. (?x. P x) <=> (?x. P (a + x))) /\
2657          (!Q. (!s. Q s) <=> (!s. Q(IMAGE (\x. a + x) s))) /\
2658          (!Q. (?s. Q s) <=> (?s. Q(IMAGE (\x. a + x) s))) /\
2659          (!Q. (!s. Q s) <=> (!s. Q(IMAGE (IMAGE (\x. a + x)) s))) /\
2660          (!Q. (?s. Q s) <=> (?s. Q(IMAGE (IMAGE (\x. a + x)) s))) /\
2661          (!P. (!g:real^1->real^N. P g) <=> (!g. P ((\x. a + x) o g))) /\
2662          (!P. (?g:real^1->real^N. P g) <=> (?g. P ((\x. a + x) o g))) /\
2663          (!P. (!g:num->real^N. P g) <=> (!g. P ((\x. a + x) o g))) /\
2664          (!P. (?g:num->real^N. P g) <=> (?g. P ((\x. a + x) o g))) /\
2665          (!Q. (!l. Q l) <=> (!l. Q(MAP (\x. a + x) l))) /\
2666          (!Q. (?l. Q l) <=> (?l. Q(MAP (\x. a + x) l)))) /\
2667         ((!P. {x | P x} = IMAGE (\x. a + x) {x | P(a + x)}) /\
2668          (!Q. {s | Q s} =
2669               IMAGE (IMAGE (\x. a + x)) {s | Q(IMAGE (\x. a + x) s)}) /\
2670          (!R. {l | R l} = IMAGE (MAP (\x. a + x)) {l | R(MAP (\x. a + x) l)}))`,
2671     GEN_TAC THEN MATCH_MP_TAC QUANTIFY_SURJECTION_HIGHER_THM THEN
2672     X_GEN_TAC `y:real^N` THEN EXISTS_TAC `y - a:real^N` THEN
2673     VECTOR_ARITH_TAC) in
2674   let GEOM_ORIGIN_CONV avoid tm =
2675     let x,tm0 = dest_forall tm in
2676     let th0 = ISPEC x pth in
2677     let x' = genvar(type_of x) in
2678     let ith = ISPEC x' qth in
2679     let th1 = PARTIAL_EXPAND_QUANTS_CONV avoid (ASSUME(concl ith)) tm0 in
2680     let th2 = CONV_RULE(RAND_CONV(SUBS_CONV(CONJUNCTS th0))) th1 in
2681     let th3 = INST[x,x'] (PROVE_HYP ith th2) in
2682     let ths = TRANSLATION_INVARIANTS x in
2683     let thr = REFL x in
2684     let th4 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV)
2685       [BETA_THM;ADD_ASSUM(concl thr) ths] th3 in
2686     let th5 = MK_FORALL x (PROVE_HYP thr th4) in
2687     GEN_REWRITE_RULE (RAND_CONV o TRY_CONV) [FORALL_SIMP] th5
2688   and GEOM_TRANSLATE_CONV avoid a tm =
2689     let cth = CONJUNCT2(ISPEC a pth)
2690     and vth = ISPEC a qth in
2691     let th1 = PARTIAL_EXPAND_QUANTS_CONV avoid (ASSUME(concl vth)) tm in
2692     let th2 = CONV_RULE(RAND_CONV(SUBS_CONV(CONJUNCTS cth))) th1 in
2693     let th3 = PROVE_HYP vth th2 in
2694     let ths = TRANSLATION_INVARIANTS a in
2695     let thr = REFL a in
2696     let th4 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV)
2697         [BETA_THM;ADD_ASSUM(concl thr) ths] th3 in
2698     PROVE_HYP thr th4 in
2699   GEOM_ORIGIN_CONV,GEOM_TRANSLATE_CONV;;
2700
2701 let GEN_GEOM_ORIGIN_TAC x avoid (asl,w as gl) =
2702   let avs,bod = strip_forall w
2703   and avs' = subtract (frees w) (freesl(map (concl o snd) asl)) in
2704   (MAP_EVERY X_GEN_TAC avs THEN
2705    MAP_EVERY (fun t -> SPEC_TAC(t,t)) (rev(subtract (avs@avs') [x])) THEN
2706    SPEC_TAC(x,x) THEN CONV_TAC(GEOM_ORIGIN_CONV avoid)) gl;;
2707
2708 let GEOM_ORIGIN_TAC x = GEN_GEOM_ORIGIN_TAC x [];;
2709
2710 let GEOM_TRANSLATE_TAC avoid (asl,w) =
2711   let a,bod = dest_forall w in
2712   let n = length(fst(strip_forall bod)) in
2713   (X_GEN_TAC a THEN
2714    CONV_TAC(funpow n BINDER_CONV (LAND_CONV(GEOM_TRANSLATE_CONV avoid a))) THEN
2715    REWRITE_TAC[]) (asl,w);;
2716
2717 (* ------------------------------------------------------------------------- *)
2718 (* Rename existential variables in conclusion to fresh genvars.              *)
2719 (* ------------------------------------------------------------------------- *)
2720
2721 let EXISTS_GENVAR_RULE =
2722   let rec rule vs th =
2723     match vs with
2724       [] -> th
2725     | v::ovs -> let x,bod = dest_exists(concl th) in
2726                 let th1 = rule ovs (ASSUME bod) in
2727                 let th2 = SIMPLE_CHOOSE x (SIMPLE_EXISTS x th1) in
2728                 PROVE_HYP th (CONV_RULE (GEN_ALPHA_CONV v) th2) in
2729   fun th -> rule (map (genvar o type_of) (fst(strip_exists(concl th)))) th;;
2730
2731 (* ------------------------------------------------------------------------- *)
2732 (* Rotate so that WLOG some point is a +ve multiple of basis vector k.       *)
2733 (* For general N, it's better to use k = 1 so the side-condition can be      *)
2734 (* discharged. For dimensions 1, 2 and 3 anything will work automatically.   *)
2735 (* Could generalize by asking the user to prove theorem 1 <= k <= N.         *)
2736 (* ------------------------------------------------------------------------- *)
2737
2738 let GEOM_BASIS_MULTIPLE_RULE =
2739   let pth = prove
2740    (`!f. orthogonal_transformation (f:real^N->real^N)
2741          ==> (vec 0 = f(vec 0) /\
2742               {} = IMAGE f {} /\
2743               {} = IMAGE (IMAGE f) {} /\
2744               (:real^N) = IMAGE f (:real^N) /\
2745               (:real^N->bool) = IMAGE (IMAGE f) (:real^N->bool) /\
2746               [] = MAP f []) /\
2747              ((!P. (!x. P x) <=> (!x. P (f x))) /\
2748               (!P. (?x. P x) <=> (?x. P (f x))) /\
2749               (!Q. (!s. Q s) <=> (!s. Q (IMAGE f s))) /\
2750               (!Q. (?s. Q s) <=> (?s. Q (IMAGE f s))) /\
2751               (!Q. (!s. Q s) <=> (!s. Q (IMAGE (IMAGE f) s))) /\
2752               (!Q. (?s. Q s) <=> (?s. Q (IMAGE (IMAGE f) s))) /\
2753               (!P. (!g:real^1->real^N. P g) <=> (!g. P (f o g))) /\
2754               (!P. (?g:real^1->real^N. P g) <=> (?g. P (f o g))) /\
2755               (!P. (!g:num->real^N. P g) <=> (!g. P (f o g))) /\
2756               (!P. (?g:num->real^N. P g) <=> (?g. P (f o g))) /\
2757               (!Q. (!l. Q l) <=> (!l. Q(MAP f l))) /\
2758               (!Q. (?l. Q l) <=> (?l. Q(MAP f l)))) /\
2759              ((!P. {x | P x} = IMAGE f {x | P(f x)}) /\
2760               (!Q. {s | Q s} = IMAGE (IMAGE f) {s | Q(IMAGE f s)}) /\
2761               (!R. {l | R l} = IMAGE (MAP f) {l | R(MAP f l)}))`,
2762     REPEAT GEN_TAC THEN DISCH_TAC THEN
2763     FIRST_ASSUM(ASSUME_TAC o
2764           MATCH_MP ORTHOGONAL_TRANSFORMATION_SURJECTIVE) THEN
2765     CONJ_TAC THENL
2766      [REWRITE_TAC[IMAGE_CLAUSES; MAP] THEN
2767       FIRST_ASSUM(ASSUME_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
2768       CONJ_TAC THENL [ASM_MESON_TAC[LINEAR_0]; ALL_TAC] THEN
2769       REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
2770       ASM_REWRITE_TAC[SURJECTIVE_IMAGE];
2771       MATCH_MP_TAC QUANTIFY_SURJECTION_HIGHER_THM THEN ASM_REWRITE_TAC[]])
2772   and oth = prove
2773    (`!f:real^N->real^N.
2774         orthogonal_transformation f /\
2775         (2 <= dimindex(:N) ==> det(matrix f) = &1)
2776         ==> linear f /\
2777             (!x y. f x = f y ==> x = y) /\
2778             (!y. ?x. f x = y) /\
2779             (!x. norm(f x) = norm x) /\
2780             (2 <= dimindex(:N) ==> det(matrix f) = &1)`,
2781     GEN_TAC THEN STRIP_TAC THEN ASM_REWRITE_TAC[] THEN REPEAT CONJ_TAC THENL
2782      [ASM_SIMP_TAC[ORTHOGONAL_TRANSFORMATION_LINEAR];
2783       ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_INJECTIVE];
2784       ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_SURJECTIVE];
2785       ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION]])
2786   and arithconv = REWRITE_CONV[DIMINDEX_1; DIMINDEX_2; DIMINDEX_3;
2787                                ARITH_RULE `1 <= 1`; DIMINDEX_GE_1] THENC
2788                   NUM_REDUCE_CONV in
2789   fun k tm ->
2790     let x,bod = dest_forall tm in
2791     let th0 = ISPECL [x; mk_small_numeral k] ROTATION_RIGHTWARD_LINE in
2792     let th1 = EXISTS_GENVAR_RULE
2793      (MP th0 (EQT_ELIM(arithconv(lhand(concl th0))))) in
2794     let [a;f],tm1 = strip_exists(concl th1) in
2795     let th_orth,th2 = CONJ_PAIR(ASSUME tm1) in
2796     let th_det,th2a = CONJ_PAIR th2 in
2797     let th_works,th_zero = CONJ_PAIR th2a in
2798     let thc,thq = CONJ_PAIR(PROVE_HYP th2 (UNDISCH(ISPEC f pth))) in
2799     let th3 = CONV_RULE(RAND_CONV(SUBS_CONV(GSYM th_works::CONJUNCTS thc)))
2800                (EXPAND_QUANTS_CONV(ASSUME(concl thq)) bod) in
2801     let th4 = PROVE_HYP thq th3 in
2802     let thps = CONJUNCTS(MATCH_MP oth (CONJ th_orth th_det)) in
2803     let th5 = LINEAR_INVARIANTS f thps in
2804     let th6 = PROVE_HYP th_orth
2805      (GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV) [BETA_THM; th5] th4) in
2806     let ntm = mk_forall(a,mk_imp(concl th_zero,rand(concl th6))) in
2807     let th7 = MP(SPEC a (ASSUME ntm)) th_zero in
2808     let th8 = DISCH ntm (EQ_MP (SYM th6) th7) in
2809     if intersect (frees(concl th8)) [a;f] = [] then
2810       let th9 = PROVE_HYP th1 (itlist SIMPLE_CHOOSE [a;f] th8) in
2811       let th10 = DISCH ntm (GEN x (UNDISCH th9)) in
2812       let a' = variant (frees(concl th10))
2813                 (mk_var(fst(dest_var x),snd(dest_var a))) in
2814       CONV_RULE(LAND_CONV (GEN_ALPHA_CONV a')) th10
2815     else
2816       let mtm = list_mk_forall([a;f],mk_imp(hd(hyp th8),rand(concl th6))) in
2817       let th9 = EQ_MP (SYM th6) (UNDISCH(SPECL [a;f] (ASSUME mtm))) in
2818       let th10 = itlist SIMPLE_CHOOSE [a;f] (DISCH mtm th9) in
2819       let th11 = GEN x (PROVE_HYP th1 th10) in
2820       MATCH_MP MONO_FORALL th11;;
2821
2822 let GEOM_BASIS_MULTIPLE_TAC k l (asl,w as gl) =
2823   let avs,bod = strip_forall w
2824   and avs' = subtract (frees w) (freesl(map (concl o snd) asl)) in
2825   (MAP_EVERY X_GEN_TAC avs THEN
2826    MAP_EVERY (fun t -> SPEC_TAC(t,t)) (rev(subtract (avs@avs') [l])) THEN
2827    SPEC_TAC(l,l) THEN
2828    W(MATCH_MP_TAC o GEOM_BASIS_MULTIPLE_RULE k o snd)) gl;;
2829
2830 (* ------------------------------------------------------------------------- *)
2831 (* Create invariance theorems automatically, in simple cases. If there are   *)
2832 (* any nested quantifiers, this will need surjectivity. It's often possible  *)
2833 (* to prove a stronger theorem by more delicate manual reasoning, so this    *)
2834 (* isn't used nearly as often as GEOM_TRANSLATE_CONV / GEOM_TRANSLATE_TAC.   *)
2835 (* As a small step, some ad-hoc rewrites analogous to FORALL_IN_IMAGE are    *)
2836 (* tried if the first step doesn't finish the goal, but it's very ad hoc.    *)
2837 (* ------------------------------------------------------------------------- *)
2838
2839 let GEOM_TRANSFORM_TAC =
2840   let cth0 = prove
2841    (`!f:real^M->real^N.
2842           linear f
2843           ==> vec 0 = f(vec 0) /\
2844               {} = IMAGE f {} /\
2845               {} = IMAGE (IMAGE f) {}`,
2846     REWRITE_TAC[IMAGE_CLAUSES] THEN MESON_TAC[LINEAR_0])
2847   and cth1 = prove
2848    (`!f:real^M->real^N.
2849           (!y. ?x. f x = y)
2850           ==> (:real^N) = IMAGE f (:real^M) /\
2851               (:real^N->bool) = IMAGE (IMAGE f) (:real^M->bool)`,
2852     REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
2853     REWRITE_TAC[SURJECTIVE_IMAGE])
2854   and sths = (CONJUNCTS o prove)
2855    (`(!f:real^M->real^N.
2856          linear f /\ (!x. norm(f x) = norm x)
2857          ==> (!x y. f x = f y ==> x = y)) /\
2858      (!f:real^N->real^N.
2859          linear f /\ (!x. norm(f x) = norm x) ==> (!y. ?x. f x = y)) /\
2860      (!f:real^N->real^N. linear f /\ (!x y. f x = f y ==> x = y)
2861                          ==> (!y. ?x. f x = y)) /\
2862      (!f:real^N->real^N. linear f /\ (!y. ?x. f x = y)
2863                          ==> (!x y. f x = f y ==> x = y))`,
2864     CONJ_TAC THENL
2865      [ONCE_REWRITE_TAC[GSYM VECTOR_SUB_EQ] THEN
2866       SIMP_TAC[GSYM LINEAR_SUB; GSYM NORM_EQ_0];
2867       MESON_TAC[ORTHOGONAL_TRANSFORMATION_SURJECTIVE;
2868                 ORTHOGONAL_TRANSFORMATION_INJECTIVE; ORTHOGONAL_TRANSFORMATION;
2869                 LINEAR_SURJECTIVE_IFF_INJECTIVE]])
2870   and aths = (CONJUNCTS o prove)
2871    (`(!f s P. (!y. y IN IMAGE f s ==> P y) <=> (!x. x IN s ==> P(f x))) /\
2872      (!f s P. (!u. u IN IMAGE (IMAGE f) s ==> P u) <=>
2873               (!t. t IN s ==> P(IMAGE f t))) /\
2874      (!f s P. (?y. y IN IMAGE f s /\ P y) <=> (?x. x IN s /\ P(f x))) /\
2875      (!f s P. (?u. u IN IMAGE (IMAGE f) s /\ P u) <=>
2876               (?t. t IN s /\ P(IMAGE f t)))`,
2877     SET_TAC[]) in
2878   fun avoid (asl,w as gl) ->
2879     let f,wff = dest_forall w in
2880     let vs,bod = strip_forall wff in
2881     let ant,cons = dest_imp bod in
2882     let hths = CONJUNCTS(ASSUME ant) in
2883     let fths = hths @ mapfilter (USABLE_CONCLUSION f hths) sths in
2884     let cths = mapfilter (USABLE_CONCLUSION f fths) [cth0; cth1]
2885     and vconv =
2886       try let vth = USABLE_CONCLUSION f fths QUANTIFY_SURJECTION_HIGHER_THM in
2887           PROVE_HYP vth o PARTIAL_EXPAND_QUANTS_CONV avoid (ASSUME(concl vth))
2888       with Failure _ -> ALL_CONV
2889     and bths = LINEAR_INVARIANTS f fths in
2890     (MAP_EVERY X_GEN_TAC (f::vs) THEN DISCH_TAC THEN
2891      GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) cths THEN
2892      CONV_TAC(LAND_CONV vconv) THEN
2893      GEN_REWRITE_TAC (LAND_CONV o REDEPTH_CONV) [bths] THEN
2894      REWRITE_TAC[] THEN
2895      REWRITE_TAC(mapfilter (ADD_ASSUM ant o ISPEC f) aths) THEN
2896      GEN_REWRITE_TAC (LAND_CONV o REDEPTH_CONV) [bths] THEN
2897      REWRITE_TAC[]) gl;;
2898
2899 (* ------------------------------------------------------------------------- *)
2900 (* Scale so that a chosen vector has size 1. Generates a conjunction of      *)
2901 (* two formulas, one for the zero case (which one hopes is trivial) and      *)
2902 (* one just like the original goal but with a norm(...) = 1 assumption.      *)
2903 (* ------------------------------------------------------------------------- *)
2904
2905 let GEOM_NORMALIZE_RULE =
2906   let pth = prove
2907    (`!a:real^N. ~(a = vec 0)
2908                 ==> vec 0 = norm(a) % vec 0 /\
2909                     a = norm(a) % inv(norm a) % a /\
2910                     {} = IMAGE (\x. norm(a) % x) {} /\
2911                     {} = IMAGE (IMAGE (\x. norm(a) % x)) {} /\
2912                     (:real^N) = IMAGE (\x. norm(a) % x) (:real^N) /\
2913                     (:real^N->bool) =
2914                     IMAGE (IMAGE (\x. norm(a) % x)) (:real^N->bool) /\
2915                     [] = MAP (\x. norm(a) % x) []`,
2916     REWRITE_TAC[IMAGE_CLAUSES; VECTOR_MUL_ASSOC; VECTOR_MUL_RZERO; MAP] THEN
2917     SIMP_TAC[NORM_EQ_0; REAL_MUL_RINV; VECTOR_MUL_LID] THEN
2918     GEN_TAC THEN DISCH_TAC THEN
2919     REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
2920     ASM_REWRITE_TAC[SURJECTIVE_IMAGE] THEN
2921     X_GEN_TAC `y:real^N` THEN EXISTS_TAC `inv(norm(a:real^N)) % y:real^N` THEN
2922     ASM_SIMP_TAC[VECTOR_MUL_ASSOC; NORM_EQ_0; REAL_MUL_RINV; VECTOR_MUL_LID])
2923   and qth = prove
2924    (`!a:real^N.
2925         ~(a = vec 0)
2926         ==> ((!P. (!r:real. P r) <=> (!r. P(norm a * r))) /\
2927              (!P. (?r:real. P r) <=> (?r. P(norm a * r))) /\
2928              (!P. (!x:real^N. P x) <=> (!x. P (norm(a) % x))) /\
2929              (!P. (?x:real^N. P x) <=> (?x. P (norm(a) % x))) /\
2930              (!Q. (!s:real^N->bool. Q s) <=>
2931                   (!s. Q(IMAGE (\x. norm(a) % x) s))) /\
2932              (!Q. (?s:real^N->bool. Q s) <=>
2933                   (?s. Q(IMAGE (\x. norm(a) % x) s))) /\
2934              (!Q. (!s:(real^N->bool)->bool. Q s) <=>
2935                   (!s. Q(IMAGE (IMAGE (\x. norm(a) % x)) s))) /\
2936              (!Q. (?s:(real^N->bool)->bool. Q s) <=>
2937                   (?s. Q(IMAGE (IMAGE (\x. norm(a) % x)) s))) /\
2938              (!P. (!g:real^1->real^N. P g) <=>
2939                   (!g. P ((\x. norm(a) % x) o g))) /\
2940              (!P. (?g:real^1->real^N. P g) <=>
2941                   (?g. P ((\x. norm(a) % x) o g))) /\
2942              (!P. (!g:num->real^N. P g) <=>
2943                   (!g. P ((\x. norm(a) % x) o g))) /\
2944              (!P. (?g:num->real^N. P g) <=>
2945                   (?g. P ((\x. norm(a) % x) o g))) /\
2946              (!Q. (!l. Q l) <=> (!l. Q(MAP (\x:real^N. norm(a) % x) l))) /\
2947              (!Q. (?l. Q l) <=> (?l. Q(MAP (\x:real^N. norm(a) % x) l)))) /\
2948             ((!P. {x:real^N | P x} =
2949                   IMAGE (\x. norm(a) % x) {x | P(norm(a) % x)}) /\
2950              (!Q. {s:real^N->bool | Q s} =
2951                   IMAGE (IMAGE (\x. norm(a) % x))
2952                        {s | Q(IMAGE (\x. norm(a) % x) s)}) /\
2953              (!R. {l:(real^N)list | R l} =
2954                   IMAGE (MAP (\x:real^N. norm(a) % x))
2955                         {l | R(MAP (\x:real^N. norm(a) % x) l)}))`,
2956     GEN_TAC THEN DISCH_TAC THEN MATCH_MP_TAC(TAUT
2957      `(a /\ b) /\ c /\ d ==> (a /\ b /\ c) /\ d`) THEN
2958     CONJ_TAC THENL
2959      [ASM_MESON_TAC[NORM_EQ_0; REAL_FIELD `~(x = &0) ==> x * inv x * a = a`];
2960       MP_TAC(ISPEC `\x:real^N. norm(a:real^N) % x`
2961         (INST_TYPE [`:real^1`,`:C`] QUANTIFY_SURJECTION_HIGHER_THM)) THEN
2962       ASM_REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN
2963       ASM_SIMP_TAC[SURJECTIVE_SCALING; NORM_EQ_0]])
2964   and lth = prove
2965    (`(!b:real^N. ~(b = vec 0) ==> (P(b) <=> Q(inv(norm b) % b)))
2966      ==> P(vec 0) /\ (!b. norm(b) = &1 ==> Q b) ==> (!b. P b)`,
2967     REPEAT STRIP_TAC THEN
2968     ASM_CASES_TAC `b:real^N = vec 0` THEN ASM_SIMP_TAC[] THEN
2969     FIRST_X_ASSUM MATCH_MP_TAC THEN
2970     ASM_SIMP_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM;
2971                  REAL_MUL_LINV; NORM_EQ_0]) in
2972   fun avoid tm ->
2973     let x,tm0 = dest_forall tm in
2974     let cth = UNDISCH(ISPEC x pth)
2975     and vth = UNDISCH(ISPEC x qth) in
2976     let th1 = ONCE_REWRITE_CONV[cth] tm0 in
2977     let th2 = CONV_RULE (RAND_CONV
2978      (PARTIAL_EXPAND_QUANTS_CONV avoid vth)) th1 in
2979     let th3 = SCALING_THEOREMS x in
2980     let th3' = (end_itlist CONJ (map
2981        (fun th -> let avs,_ = strip_forall(concl th) in
2982                   let gvs = map (genvar o type_of) avs in
2983                   GENL gvs (SPECL gvs th))
2984        (CONJUNCTS th3))) in
2985     let th4 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV)
2986                [BETA_THM; th3'] th2 in
2987     MATCH_MP lth (GEN x (DISCH_ALL th4));;
2988
2989 let GEN_GEOM_NORMALIZE_TAC x avoid (asl,w as gl) =
2990   let avs,bod = strip_forall w
2991   and avs' = subtract (frees w) (freesl(map (concl o snd) asl)) in
2992   (MAP_EVERY X_GEN_TAC avs THEN
2993    MAP_EVERY (fun t -> SPEC_TAC(t,t)) (rev(subtract (avs@avs') [x])) THEN
2994    SPEC_TAC(x,x) THEN
2995    W(MATCH_MP_TAC o GEOM_NORMALIZE_RULE avoid o snd)) gl;;
2996
2997 let GEOM_NORMALIZE_TAC x = GEN_GEOM_NORMALIZE_TAC x [];;
2998
2999 (* ------------------------------------------------------------------------- *)
3000 (* Add invariance theorems for collinearity.                                 *)
3001 (* ------------------------------------------------------------------------- *)
3002
3003 let COLLINEAR_TRANSLATION_EQ = prove
3004  (`!a s. collinear (IMAGE (\x. a + x) s) <=> collinear s`,
3005   REWRITE_TAC[collinear] THEN GEOM_TRANSLATE_TAC["u"]);;
3006
3007 add_translation_invariants [COLLINEAR_TRANSLATION_EQ];;
3008
3009 let COLLINEAR_TRANSLATION = prove
3010  (`!s a. collinear s ==> collinear (IMAGE (\x. a + x) s)`,
3011   REWRITE_TAC[COLLINEAR_TRANSLATION_EQ]);;
3012
3013 let COLLINEAR_LINEAR_IMAGE = prove
3014  (`!f s. collinear s /\ linear f ==> collinear(IMAGE f s)`,
3015   REPEAT GEN_TAC THEN DISCH_THEN(CONJUNCTS_THEN2 MP_TAC ASSUME_TAC) THEN
3016   REWRITE_TAC[collinear; IMP_CONJ; RIGHT_FORALL_IMP_THM; FORALL_IN_IMAGE] THEN
3017   ASM_MESON_TAC[LINEAR_SUB; LINEAR_CMUL]);;
3018
3019 let COLLINEAR_LINEAR_IMAGE_EQ = prove
3020  (`!f s. linear f /\ (!x y. f x = f y ==> x = y)
3021          ==> (collinear (IMAGE f s) <=> collinear s)`,
3022   MATCH_ACCEPT_TAC(LINEAR_INVARIANT_RULE COLLINEAR_LINEAR_IMAGE));;
3023
3024 add_linear_invariants [COLLINEAR_LINEAR_IMAGE_EQ];;
3025
3026 (* ------------------------------------------------------------------------- *)
3027 (* Take a theorem "th" with outer universal quantifiers involving real^N     *)
3028 (* and a theorem "dth" asserting |- dimindex(:M) <= dimindex(:N) and         *)
3029 (* return a theorem replacing type :N by :M in th. Neither N or M need be a  *)
3030 (* type variable.                                                            *)
3031 (* ------------------------------------------------------------------------- *)
3032
3033 let GEOM_DROP_DIMENSION_RULE =
3034   let oth = prove
3035    (`!f:real^M->real^N.
3036           linear f /\ (!x. norm(f x) = norm x)
3037           ==> linear f /\
3038               (!x y. f x = f y ==> x = y) /\
3039               (!x. norm(f x) = norm x)`,
3040     MESON_TAC[PRESERVES_NORM_INJECTIVE])
3041   and cth = prove
3042    (`linear(f:real^M->real^N)
3043      ==> vec 0 = f(vec 0) /\
3044          {} = IMAGE f {} /\
3045          {} = IMAGE (IMAGE f) {} /\
3046          [] = MAP f []`,
3047     REWRITE_TAC[IMAGE_CLAUSES; MAP; GSYM LINEAR_0]) in
3048   fun dth th ->
3049     let ath = GEN_ALL th
3050     and eth = MATCH_MP ISOMETRY_UNIV_UNIV dth
3051     and avoid = variables(concl th) in
3052     let f,bod = dest_exists(concl eth) in
3053     let fimage = list_mk_icomb "IMAGE" [f]
3054     and fmap = list_mk_icomb "MAP" [f]
3055     and fcompose = list_mk_icomb "o" [f] in
3056     let fimage2 = list_mk_icomb "IMAGE" [fimage] in
3057     let lin,iso = CONJ_PAIR(ASSUME bod) in
3058     let olduniv = rand(rand(concl dth))
3059     and newuniv = rand(lhand(concl dth)) in
3060     let oldty = fst(dest_fun_ty(type_of olduniv))
3061     and newty = fst(dest_fun_ty(type_of newuniv)) in
3062     let newvar v =
3063        let n,t = dest_var v in
3064        variant avoid (mk_var(n,tysubst[newty,oldty] t)) in
3065     let newterm v =
3066       try let v' = newvar v in
3067           tryfind (fun f -> mk_comb(f,v')) [f;fimage;fmap;fcompose;fimage2]
3068       with Failure _ -> v in
3069     let specrule th =
3070       let v = fst(dest_forall(concl th)) in SPEC (newterm v) th in
3071     let sth = SUBS(CONJUNCTS(MATCH_MP cth lin)) ath in
3072     let fth = SUBS[SYM(MATCH_MP LINEAR_0 lin)] (repeat specrule sth) in
3073     let thps = CONJUNCTS(MATCH_MP oth (ASSUME bod)) in
3074     let th5 = LINEAR_INVARIANTS f thps in
3075     let th6 = GEN_REWRITE_RULE REDEPTH_CONV [th5] fth in
3076     let th7 = PROVE_HYP eth (SIMPLE_CHOOSE f th6) in
3077     GENL (map newvar (fst(strip_forall(concl ath)))) th7;;
3078
3079 (* ------------------------------------------------------------------------- *)
3080 (* Transfer theorems automatically between same-dimension spaces.            *)
3081 (* Given dth = A |- dimindex(:M) = dimindex(:N)                              *)
3082 (* and a theorem th involving variables of type real^N                       *)
3083 (* returns a corresponding theorem mapped to type real^M with assumptions A. *)
3084 (* ------------------------------------------------------------------------- *)
3085
3086 let GEOM_EQUAL_DIMENSION_RULE =
3087   let bth = prove
3088    (`dimindex(:M) = dimindex(:N)
3089      ==> ?f:real^M->real^N.
3090              (linear f /\ (!y. ?x. f x = y)) /\
3091              (!x. norm(f x) = norm x)`,
3092     REWRITE_TAC[SET_RULE `(!y. ?x. f x = y) <=> IMAGE f UNIV = UNIV`] THEN
3093     DISCH_TAC THEN REWRITE_TAC[GSYM CONJ_ASSOC] THEN
3094     MATCH_MP_TAC ISOMETRY_UNIV_SUBSPACE THEN
3095     REWRITE_TAC[SUBSPACE_UNIV; DIM_UNIV] THEN FIRST_ASSUM ACCEPT_TAC)
3096   and pth = prove
3097    (`!f:real^M->real^N.
3098         linear f /\ (!y. ?x. f x = y)
3099          ==> (vec 0 = f(vec 0) /\
3100               {} = IMAGE f {} /\
3101               {} = IMAGE (IMAGE f) {} /\
3102               (:real^N) = IMAGE f (:real^M) /\
3103               (:real^N->bool) = IMAGE (IMAGE f) (:real^M->bool) /\
3104               [] = MAP f []) /\
3105              ((!P. (!x. P x) <=> (!x. P (f x))) /\
3106               (!P. (?x. P x) <=> (?x. P (f x))) /\
3107               (!Q. (!s. Q s) <=> (!s. Q (IMAGE f s))) /\
3108               (!Q. (?s. Q s) <=> (?s. Q (IMAGE f s))) /\
3109               (!Q. (!s. Q s) <=> (!s. Q (IMAGE (IMAGE f) s))) /\
3110               (!Q. (?s. Q s) <=> (?s. Q (IMAGE (IMAGE f) s))) /\
3111               (!P. (!g:real^1->real^N. P g) <=> (!g. P (f o g))) /\
3112               (!P. (?g:real^1->real^N. P g) <=> (?g. P (f o g))) /\
3113               (!P. (!g:num->real^N. P g) <=> (!g. P (f o g))) /\
3114               (!P. (?g:num->real^N. P g) <=> (?g. P (f o g))) /\
3115               (!Q. (!l. Q l) <=> (!l. Q(MAP f l))) /\
3116               (!Q. (?l. Q l) <=> (?l. Q(MAP f l)))) /\
3117              ((!P. {x | P x} = IMAGE f {x | P(f x)}) /\
3118               (!Q. {s | Q s} = IMAGE (IMAGE f) {s | Q(IMAGE f s)}) /\
3119               (!R. {l | R l} = IMAGE (MAP f) {l | R(MAP f l)}))`,
3120     GEN_TAC THEN
3121     SIMP_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> (!y. ?x. f x = y)`;
3122              SURJECTIVE_IMAGE] THEN
3123     MATCH_MP_TAC MONO_AND THEN
3124     REWRITE_TAC[QUANTIFY_SURJECTION_HIGHER_THM] THEN
3125     REWRITE_TAC[IMAGE_CLAUSES; MAP] THEN MESON_TAC[LINEAR_0]) in
3126   fun dth th ->
3127     let eth = EXISTS_GENVAR_RULE (MATCH_MP bth dth) in
3128     let f,bod = dest_exists(concl eth) in
3129     let lsth,neth = CONJ_PAIR(ASSUME bod) in
3130     let cth,qth = CONJ_PAIR(MATCH_MP pth lsth) in
3131     let th1 = CONV_RULE
3132      (EXPAND_QUANTS_CONV qth THENC SUBS_CONV(CONJUNCTS cth)) th in
3133     let ith = LINEAR_INVARIANTS f (neth::CONJUNCTS lsth) in
3134     let th2 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV) [BETA_THM;ith] th1 in
3135     let th3 = GEN f (DISCH bod th2) in
3136     MP (CONV_RULE (REWR_CONV LEFT_FORALL_IMP_THM) th3) eth;;