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