development/thales/chaff/tmp/determinants_patch.ml

   1 (* ========================================================================= *)
   2 (* Determinant and trace of a square matrix.                                 *)
   3 (*                                                                           *)
   4 (*              (c) Copyright, John Harrison 1998-2008                       *)
   5 (* ========================================================================= *)
   6
   7 (* Feb 2, 2012, patch of a file not yet committed to SVN.
   8     Sent by email by Harrison.  Delete in a few weeks after Harrison commits changes *)
   9
  10
  11 overload_interface ("--",`(matrix_neg):real^N^M->real^N^M`);;
  12 overload_interface ("+",`(matrix_add):real^N^M->real^N^M->real^N^M`);;
  13 overload_interface ("-",`(matrix_sub):real^N^M->real^N^M->real^N^M`);;
  14
  15 make_overloadable "**" `:A->B->C`;;
  16
  17 overload_interface ("**",`(matrix_mul):real^N^M->real^P^N->real^P^M`);;
  18 overload_interface ("**",`(matrix_vector_mul):real^N^M->real^N->real^M`);;
  19 overload_interface ("**",`(vector_matrix_mul):real^M->real^N^M->real^N`);;
  20
  21 parse_as_infix("%%",(21,"right"));;
  22
  23
  24 needs "Multivariate/vectors.ml";;
  25 needs "Library/permutations.ml";;
  26 needs "Library/floor.ml";;
  27 needs "Library/products.ml";;
  28
  29
  30
  31 prioritize_real();;
  32
  33 (* ------------------------------------------------------------------------- *)
  34 (* Trace of a matrix (this is relatively easy).                              *)
  35 (* ------------------------------------------------------------------------- *)
  36
  37 let trace = new_definition
  38   `(trace:real^N^N->real) A = sum(1..dimindex(:N)) (\i. A$i$i)`;;
  39
  40 let TRACE_0 = prove
  41  (`trace(mat 0) = &0`,
  42   SIMP_TAC[trace; mat; LAMBDA_BETA; SUM_0]);;
  43
  44 let TRACE_I = prove
  45  (`trace(mat 1 :real^N^N) = &(dimindex(:N))`,
  46   SIMP_TAC[trace; mat; LAMBDA_BETA; SUM_CONST_NUMSEG; REAL_MUL_RID] THEN
  47   AP_TERM_TAC THEN ARITH_TAC);;
  48
  49 let TRACE_ADD = prove
  50  (`!A B:real^N^N. trace(A + B) = trace(A) + trace(B)`,
  51   SIMP_TAC[trace; matrix_add; SUM_ADD_NUMSEG; LAMBDA_BETA]);;
  52
  53 let TRACE_SUB = prove
  54  (`!A B:real^N^N. trace(A - B) = trace(A) - trace(B)`,
  55   SIMP_TAC[trace; matrix_sub; SUM_SUB_NUMSEG; LAMBDA_BETA]);;
  56
  57 let TRACE_MUL_SYM = prove
  58  (`!A B:real^N^N. trace(A ** B) = trace(B ** A)`,
  59   REPEAT GEN_TAC THEN SIMP_TAC[trace; matrix_mul; LAMBDA_BETA] THEN
  60   GEN_REWRITE_TAC RAND_CONV [SUM_SWAP_NUMSEG] THEN REWRITE_TAC[REAL_MUL_SYM]);;
  61
  62 (* ------------------------------------------------------------------------- *)
  63 (* Definition of determinant.                                                *)
  64 (* ------------------------------------------------------------------------- *)
  65
  66 let det = new_definition
  67  `det(A:real^N^N) =
  68         sum { p | p permutes 1..dimindex(:N) }
  69             (\p. sign(p) * product (1..dimindex(:N)) (\i. A$i$(p i)))`;;
  70
  71 (* ------------------------------------------------------------------------- *)
  72 (* A few general lemmas we need below.                                       *)
  73 (* ------------------------------------------------------------------------- *)
  74
  75 let IN_DIMINDEX_SWAP = prove
  76  (`!m n j. 1 <= m /\ m <= dimindex(:N) /\
  77              1 <= n /\ n <= dimindex(:N) /\
  78              1 <= j /\ j <= dimindex(:N)
  79            ==> 1 <= swap(m,n) j /\ swap(m,n) j <= dimindex(:N)`,
  80   REWRITE_TAC[swap] THEN ARITH_TAC);;
  81
  82 let LAMBDA_BETA_PERM = prove
  83  (`!p i. p permutes 1..dimindex(:N) /\ 1 <= i /\ i <= dimindex(:N)
  84          ==> ((lambda) g :A^N) $ p(i) = g(p i)`,
  85   ASM_MESON_TAC[LAMBDA_BETA; PERMUTES_IN_IMAGE; IN_NUMSEG]);;
  86
  87 let PRODUCT_PERMUTE = prove
  88  (`!f p s. p permutes s ==> product s f = product s (f o p)`,
  89   REWRITE_TAC[product] THEN MATCH_MP_TAC ITERATE_PERMUTE THEN
  90   REWRITE_TAC[MONOIDAL_REAL_MUL]);;
  91
  92 let PRODUCT_PERMUTE_NUMSEG = prove
  93  (`!f p m n. p permutes m..n ==> product(m..n) f = product(m..n) (f o p)`,
  94   MESON_TAC[PRODUCT_PERMUTE; FINITE_NUMSEG]);;
  95
  96 let REAL_MUL_SUM = prove
  97  (`!s t f g.
  98         FINITE s /\ FINITE t
  99         ==> sum s f * sum t g = sum s (\i. sum t (\j. f(i) * g(j)))`,
 100   SIMP_TAC[SUM_LMUL] THEN
 101   ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN SIMP_TAC[SUM_LMUL]);;
 102
 103 let REAL_MUL_SUM_NUMSEG = prove
 104  (`!m n p q. sum(m..n) f * sum(p..q) g =
 105              sum(m..n) (\i. sum(p..q) (\j. f(i) * g(j)))`,
 106   SIMP_TAC[REAL_MUL_SUM; FINITE_NUMSEG]);;
 107
 108 (* ------------------------------------------------------------------------- *)
 109 (* Basic determinant properties.                                             *)
 110 (* ------------------------------------------------------------------------- *)
 111
 112 let DET_TRANSP = prove
 113  (`!A:real^N^N. det(transp A) = det A`,
 114   GEN_TAC THEN REWRITE_TAC[det] THEN
 115   GEN_REWRITE_TAC LAND_CONV [SUM_PERMUTATIONS_INVERSE] THEN
 116   MATCH_MP_TAC SUM_EQ THEN
 117   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN X_GEN_TAC `p:num->num` THEN
 118   REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN BINOP_TAC THENL
 119    [ASM_MESON_TAC[SIGN_INVERSE; PERMUTATION_PERMUTES; FINITE_NUMSEG];
 120     ALL_TAC] THEN
 121   FIRST_ASSUM(fun th -> GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV)
 122     [GSYM(MATCH_MP PERMUTES_IMAGE th)]) THEN
 123   MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC
 124    `product(1..dimindex(:N))
 125        ((\i. (transp A:real^N^N)$i$inverse p(i)) o p)` THEN
 126   CONJ_TAC THENL
 127    [MATCH_MP_TAC PRODUCT_IMAGE THEN
 128     ASM_MESON_TAC[FINITE_NUMSEG; PERMUTES_INJECTIVE; PERMUTES_INVERSE];
 129     MATCH_MP_TAC PRODUCT_EQ THEN REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG] THEN
 130     SIMP_TAC[transp; LAMBDA_BETA; o_THM] THEN
 131     FIRST_ASSUM(MP_TAC o MATCH_MP PERMUTES_INVERSES_o) THEN
 132     SIMP_TAC[FUN_EQ_THM; I_THM; o_THM] THEN STRIP_TAC THEN
 133     ASM_SIMP_TAC[PERMUTES_IN_NUMSEG; LAMBDA_BETA_PERM; LAMBDA_BETA]]);;
 134
 135 let DET_LOWERTRIANGULAR = prove
 136  (`!A:real^N^N.
 137         (!i j. 1 <= i /\ i <= dimindex(:N) /\
 138                1 <= j /\ j <= dimindex(:N) /\ i < j ==> A$i$j = &0)
 139         ==> det(A) = product(1..dimindex(:N)) (\i. A$i$i)`,
 140   REPEAT STRIP_TAC THEN REWRITE_TAC[det] THEN MATCH_MP_TAC EQ_TRANS THEN
 141   EXISTS_TAC `sum {I}
 142      (\p. sign p * product(1..dimindex(:N)) (\i. (A:real^N^N)$i$p(i)))` THEN
 143   CONJ_TAC THENL
 144    [ALL_TAC; REWRITE_TAC[SUM_SING; SIGN_I; REAL_MUL_LID; I_THM]] THEN
 145   MATCH_MP_TAC SUM_SUPERSET THEN
 146   SIMP_TAC[IN_SING; FINITE_RULES; SUBSET; IN_ELIM_THM; PERMUTES_I] THEN
 147   X_GEN_TAC `p:num->num` THEN STRIP_TAC THEN
 148   ASM_REWRITE_TAC[PRODUCT_EQ_0_NUMSEG; REAL_ENTIRE; SIGN_NZ] THEN
 149   MP_TAC(SPECL [`p:num->num`; `1..dimindex(:N)`] PERMUTES_NUMSET_LE) THEN
 150   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; NOT_LT]);;
 151
 152 let DET_UPPERTRIANGULAR = prove
 153  (`!A:real^N^N.
 154         (!i j. 1 <= i /\ i <= dimindex(:N) /\
 155                1 <= j /\ j <= dimindex(:N) /\ j < i ==> A$i$j = &0)
 156         ==> det(A) = product(1..dimindex(:N)) (\i. A$i$i)`,
 157   REPEAT STRIP_TAC THEN REWRITE_TAC[det] THEN MATCH_MP_TAC EQ_TRANS THEN
 158   EXISTS_TAC `sum {I}
 159      (\p. sign p * product(1..dimindex(:N)) (\i. (A:real^N^N)$i$p(i)))` THEN
 160   CONJ_TAC THENL
 161    [ALL_TAC; REWRITE_TAC[SUM_SING; SIGN_I; REAL_MUL_LID; I_THM]] THEN
 162   MATCH_MP_TAC SUM_SUPERSET THEN
 163   SIMP_TAC[IN_SING; FINITE_RULES; SUBSET; IN_ELIM_THM; PERMUTES_I] THEN
 164   X_GEN_TAC `p:num->num` THEN STRIP_TAC THEN
 165   ASM_REWRITE_TAC[PRODUCT_EQ_0_NUMSEG; REAL_ENTIRE; SIGN_NZ] THEN
 166   MP_TAC(SPECL [`p:num->num`; `1..dimindex(:N)`] PERMUTES_NUMSET_GE) THEN
 167   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; NOT_LT]);;
 168
 169 let DET_DIAGONAL = prove
 170  (`!A:real^N^N.
 171         (!i j. 1 <= i /\ i <= dimindex(:N) /\
 172                1 <= j /\ j <= dimindex(:N) /\ ~(i = j) ==> A$i$j = &0)
 173         ==> det(A) = product(1..dimindex(:N)) (\i. A$i$i)`,
 174   REPEAT STRIP_TAC THEN MATCH_MP_TAC DET_LOWERTRIANGULAR THEN
 175   REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN ASM_REWRITE_TAC[] THEN
 176   ASM_MESON_TAC[LT_REFL]);;
 177
 178 let DET_I = prove
 179  (`det(mat 1 :real^N^N) = &1`,
 180   MATCH_MP_TAC EQ_TRANS THEN
 181   EXISTS_TAC `product(1..dimindex(:N)) (\i. (mat 1:real^N^N)$i$i)` THEN
 182   CONJ_TAC THENL
 183    [MATCH_MP_TAC DET_LOWERTRIANGULAR;
 184     MATCH_MP_TAC PRODUCT_EQ_1_NUMSEG] THEN
 185   SIMP_TAC[mat; LAMBDA_BETA] THEN MESON_TAC[LT_REFL]);;
 186
 187 let DET_0 = prove
 188  (`det(mat 0 :real^N^N) = &0`,
 189   MATCH_MP_TAC EQ_TRANS THEN
 190   EXISTS_TAC `product(1..dimindex(:N)) (\i. (mat 0:real^N^N)$i$i)` THEN
 191   CONJ_TAC THENL
 192    [MATCH_MP_TAC DET_LOWERTRIANGULAR;
 193     REWRITE_TAC[PRODUCT_EQ_0_NUMSEG] THEN EXISTS_TAC `1`] THEN
 194   SIMP_TAC[mat; LAMBDA_BETA; COND_ID; DIMINDEX_GE_1; LE_REFL]);;
 195
 196 let DET_PERMUTE_ROWS = prove
 197  (`!A:real^N^N p.
 198         p permutes 1..dimindex(:N)
 199         ==> det(lambda i. A$p(i)) = sign(p) * det(A)`,
 200   REWRITE_TAC[det] THEN SIMP_TAC[LAMBDA_BETA] THEN REPEAT STRIP_TAC THEN
 201   SIMP_TAC[GSYM SUM_LMUL; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 202   FIRST_ASSUM(fun th -> GEN_REWRITE_TAC LAND_CONV
 203     [MATCH_MP SUM_PERMUTATIONS_COMPOSE_R th]) THEN
 204   MATCH_MP_TAC SUM_EQ THEN
 205   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN X_GEN_TAC `q:num->num` THEN
 206   REWRITE_TAC[IN_ELIM_THM; REAL_MUL_ASSOC] THEN DISCH_TAC THEN BINOP_TAC THENL
 207    [ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
 208     ASM_MESON_TAC[SIGN_COMPOSE; PERMUTATION_PERMUTES; FINITE_NUMSEG];
 209     ALL_TAC] THEN
 210   MP_TAC(MATCH_MP PERMUTES_INVERSE (ASSUME `p permutes 1..dimindex(:N)`)) THEN
 211   DISCH_THEN(fun th -> GEN_REWRITE_TAC LAND_CONV
 212     [MATCH_MP PRODUCT_PERMUTE_NUMSEG th]) THEN
 213   MATCH_MP_TAC PRODUCT_EQ THEN REWRITE_TAC[IN_NUMSEG; FINITE_NUMSEG] THEN
 214   REPEAT STRIP_TAC THEN REWRITE_TAC[o_THM] THEN
 215   ASM_MESON_TAC[PERMUTES_INVERSES]);;
 216
 217 let DET_PERMUTE_COLUMNS = prove
 218  (`!A:real^N^N p.
 219         p permutes 1..dimindex(:N)
 220         ==> det((lambda i j. A$i$p(j)):real^N^N) = sign(p) * det(A)`,
 221   REPEAT STRIP_TAC THEN
 222   GEN_REWRITE_TAC (funpow 2 RAND_CONV) [GSYM DET_TRANSP] THEN
 223   FIRST_ASSUM(fun th -> ONCE_REWRITE_TAC
 224    [GSYM(MATCH_MP DET_PERMUTE_ROWS th)]) THEN
 225   GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN AP_TERM_TAC THEN
 226   ASM_SIMP_TAC[CART_EQ; transp; LAMBDA_BETA; LAMBDA_BETA_PERM]);;
 227
 228 let DET_IDENTICAL_ROWS = prove
 229  (`!A:real^N^N i j. 1 <= i /\ i <= dimindex(:N) /\
 230                     1 <= j /\ j <= dimindex(:N) /\ ~(i = j) /\
 231                     row i A = row j A
 232                     ==> det A = &0`,
 233   REPEAT STRIP_TAC THEN
 234   MP_TAC(SPECL [`A:real^N^N`; `swap(i:num,j:num)`] DET_PERMUTE_ROWS) THEN
 235   ASM_SIMP_TAC[PERMUTES_SWAP; IN_NUMSEG; SIGN_SWAP] THEN
 236   MATCH_MP_TAC(REAL_ARITH `a = b ==> b = -- &1 * a ==> a = &0`) THEN
 237   AP_TERM_TAC THEN FIRST_X_ASSUM(MP_TAC o SYM) THEN
 238   SIMP_TAC[row; CART_EQ; LAMBDA_BETA] THEN
 239   REWRITE_TAC[swap] THEN ASM_MESON_TAC[]);;
 240
 241 let DET_IDENTICAL_COLUMNS = prove
 242  (`!A:real^N^N i j. 1 <= i /\ i <= dimindex(:N) /\
 243                     1 <= j /\ j <= dimindex(:N) /\ ~(i = j) /\
 244                     column i A = column j A
 245                     ==> det A = &0`,
 246   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[GSYM DET_TRANSP] THEN
 247   MATCH_MP_TAC DET_IDENTICAL_ROWS THEN ASM_MESON_TAC[ROW_TRANSP]);;
 248
 249 let DET_ZERO_ROW = prove
 250  (`!A:real^N^N i.
 251        1 <= i /\ i <= dimindex(:N) /\ row i A = vec 0  ==> det A = &0`,
 252   SIMP_TAC[det; row; CART_EQ; LAMBDA_BETA; VEC_COMPONENT] THEN
 253   REPEAT STRIP_TAC THEN MATCH_MP_TAC SUM_EQ_0 THEN
 254   REWRITE_TAC[IN_ELIM_THM; REAL_ENTIRE; SIGN_NZ] THEN REPEAT STRIP_TAC THEN
 255   SIMP_TAC[PRODUCT_EQ_0; FINITE_NUMSEG; IN_NUMSEG] THEN
 256   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]);;
 257
 258 let DET_ZERO_COLUMN = prove
 259  (`!A:real^N^N i.
 260        1 <= i /\ i <= dimindex(:N) /\ column i A = vec 0  ==> det A = &0`,
 261   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[GSYM DET_TRANSP] THEN
 262   MATCH_MP_TAC DET_ZERO_ROW THEN ASM_MESON_TAC[ROW_TRANSP]);;
 263
 264 let DET_ROW_ADD = prove
 265  (`!a b c k.
 266          1 <= k /\ k <= dimindex(:N)
 267          ==> det((lambda i. if i = k then a + b else c i):real^N^N) =
 268              det((lambda i. if i = k then a else c i):real^N^N) +
 269              det((lambda i. if i = k then b else c i):real^N^N)`,
 270   SIMP_TAC[det; LAMBDA_BETA; GSYM SUM_ADD;
 271            FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 272   REPEAT STRIP_TAC THEN MATCH_MP_TAC SUM_EQ THEN
 273   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 274   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN
 275   DISCH_TAC THEN REWRITE_TAC[GSYM REAL_ADD_LDISTRIB] THEN AP_TERM_TAC THEN
 276   SUBGOAL_THEN `1..dimindex(:N) = k INSERT ((1..dimindex(:N)) DELETE k)`
 277   SUBST1_TAC THENL [ASM_MESON_TAC[INSERT_DELETE; IN_NUMSEG]; ALL_TAC] THEN
 278   SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE] THEN
 279   MATCH_MP_TAC(REAL_RING
 280    `c = a + b /\ y = x:real /\ z = x ==> c * x = a * y + b * z`) THEN
 281   REWRITE_TAC[VECTOR_ADD_COMPONENT] THEN
 282   CONJ_TAC THEN MATCH_MP_TAC PRODUCT_EQ THEN
 283   SIMP_TAC[IN_DELETE; FINITE_DELETE; FINITE_NUMSEG]);;
 284
 285 let DET_ROW_MUL = prove
 286  (`!a b c k.
 287         1 <= k /\ k <= dimindex(:N)
 288         ==> det((lambda i. if i = k then c % a else b i):real^N^N) =
 289             c * det((lambda i. if i = k then a else b i):real^N^N)`,
 290   SIMP_TAC[det; LAMBDA_BETA; GSYM SUM_LMUL;
 291            FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 292   REPEAT STRIP_TAC THEN MATCH_MP_TAC SUM_EQ THEN
 293   SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 294   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 295   SUBGOAL_THEN `1..dimindex(:N) = k INSERT ((1..dimindex(:N)) DELETE k)`
 296   SUBST1_TAC THENL [ASM_MESON_TAC[INSERT_DELETE; IN_NUMSEG]; ALL_TAC] THEN
 297   SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE] THEN
 298   MATCH_MP_TAC(REAL_RING
 299    `cp = c * p /\ p1 = p2:real ==> s * cp * p1 = c * s * p * p2`) THEN
 300   REWRITE_TAC[VECTOR_MUL_COMPONENT] THEN MATCH_MP_TAC PRODUCT_EQ THEN
 301   SIMP_TAC[IN_DELETE; FINITE_DELETE; FINITE_NUMSEG]);;
 302
 303 let DET_ROW_OPERATION = prove
 304  (`!A:real^N^N i.
 305         1 <= i /\ i <= dimindex(:N) /\
 306         1 <= j /\ j <= dimindex(:N) /\ ~(i = j)
 307         ==> det(lambda k. if k = i then row i A + c % row j A else row k A) =
 308             det A`,
 309   REPEAT STRIP_TAC THEN ASM_SIMP_TAC[DET_ROW_ADD; DET_ROW_MUL] THEN
 310   MATCH_MP_TAC(REAL_RING `a = b /\ d = &0 ==> a + c * d = b`) THEN
 311   CONJ_TAC THENL
 312    [AP_TERM_TAC THEN ASM_SIMP_TAC[LAMBDA_BETA; CART_EQ] THEN
 313     REPEAT STRIP_TAC THEN COND_CASES_TAC THEN
 314     ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ];
 315     MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 316     MAP_EVERY EXISTS_TAC [`i:num`; `j:num`] THEN
 317     ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ]]);;
 318
 319 let DET_ROW_SPAN = prove
 320  (`!A:real^N^N i x.
 321         1 <= i /\ i <= dimindex(:N) /\
 322         x IN span {row j A | 1 <= j /\ j <= dimindex(:N) /\ ~(j = i)}
 323         ==> det(lambda k. if k = i then row i A + x else row k A) =
 324             det A`,
 325   GEN_TAC THEN GEN_TAC THEN
 326   REWRITE_TAC[IMP_CONJ; RIGHT_FORALL_IMP_THM] THEN REPEAT DISCH_TAC THEN
 327   MATCH_MP_TAC SPAN_INDUCT_ALT THEN CONJ_TAC THENL
 328    [AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA; VECTOR_ADD_RID] THEN
 329     REPEAT STRIP_TAC THEN COND_CASES_TAC THEN ASM_SIMP_TAC[row; LAMBDA_BETA];
 330     ALL_TAC] THEN
 331   REPEAT GEN_TAC THEN REWRITE_TAC[IN_ELIM_THM] THEN
 332   DISCH_THEN(CONJUNCTS_THEN2 (X_CHOOSE_TAC `j:num`) (SUBST_ALL_TAC o SYM)) THEN
 333   ONCE_REWRITE_TAC[VECTOR_ARITH
 334      `a + c % x + y:real^N = (a + y) + c % x`] THEN
 335   ABBREV_TAC `z = row i (A:real^N^N) + y` THEN
 336   ASM_SIMP_TAC[DET_ROW_MUL; DET_ROW_ADD] THEN
 337   MATCH_MP_TAC(REAL_RING `d = &0 ==> a + c * d = a`) THEN
 338   MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 339   MAP_EVERY EXISTS_TAC [`i:num`; `j:num`] THEN
 340   ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ]);;
 341
 342 (* ------------------------------------------------------------------------- *)
 343 (* May as well do this, though it's a bit unsatisfactory since it ignores    *)
 344 (* exact duplicates by considering the rows/columns as a set.                *)
 345 (* ------------------------------------------------------------------------- *)
 346
 347 let DET_DEPENDENT_ROWS = prove
 348  (`!A:real^N^N. dependent(rows A) ==> det A = &0`,
 349   GEN_TAC THEN
 350   REWRITE_TAC[dependent; rows; IN_ELIM_THM; LEFT_AND_EXISTS_THM] THEN
 351   REWRITE_TAC[LEFT_IMP_EXISTS_THM] THEN GEN_TAC THEN X_GEN_TAC `i:num` THEN
 352   STRIP_TAC THEN FIRST_X_ASSUM SUBST_ALL_TAC THEN
 353   ASM_CASES_TAC
 354    `?i j. 1 <= i /\ i <= dimindex(:N) /\
 355           1 <= j /\ j <= dimindex(:N) /\ ~(i = j) /\
 356           row i (A:real^N^N) = row j A`
 357   THENL [ASM_MESON_TAC[DET_IDENTICAL_ROWS]; ALL_TAC] THEN
 358   MP_TAC(SPECL [`A:real^N^N`; `i:num`; `--(row i (A:real^N^N))`]
 359     DET_ROW_SPAN) THEN
 360   ANTS_TAC THENL
 361    [ASM_REWRITE_TAC[] THEN MATCH_MP_TAC SPAN_NEG THEN
 362     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [IN]) THEN
 363     MATCH_MP_TAC(TAUT `a = b ==> a ==> b`) THEN
 364     REWRITE_TAC[IN] THEN AP_THM_TAC THEN AP_TERM_TAC THEN
 365     REWRITE_TAC[EXTENSION; IN_DELETE; IN_ELIM_THM] THEN ASM_MESON_TAC[];
 366     DISCH_THEN(SUBST1_TAC o SYM) THEN MATCH_MP_TAC DET_ZERO_ROW THEN
 367     EXISTS_TAC `i:num` THEN
 368     ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ; VECTOR_ADD_COMPONENT;
 369                  VECTOR_NEG_COMPONENT; VEC_COMPONENT] THEN
 370     REAL_ARITH_TAC]);;
 371
 372 let DET_DEPENDENT_COLUMNS = prove
 373  (`!A:real^N^N. dependent(columns A) ==> det A = &0`,
 374   MESON_TAC[DET_DEPENDENT_ROWS; ROWS_TRANSP; DET_TRANSP]);;
 375
 376 (* ------------------------------------------------------------------------- *)
 377 (* Multilinearity and the multiplication formula.                            *)
 378 (* ------------------------------------------------------------------------- *)
 379
 380 let DET_LINEAR_ROW_VSUM = prove
 381  (`!a c s k.
 382          FINITE s /\ 1 <= k /\ k <= dimindex(:N)
 383          ==> det((lambda i. if i = k then vsum s a else c i):real^N^N) =
 384              sum s
 385                (\j. det((lambda i. if i = k then a(j) else c i):real^N^N))`,
 386   GEN_TAC THEN GEN_TAC THEN ONCE_REWRITE_TAC[IMP_CONJ] THEN
 387   REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN
 388   MATCH_MP_TAC FINITE_INDUCT_STRONG THEN
 389   SIMP_TAC[VSUM_CLAUSES; SUM_CLAUSES; DET_ROW_ADD] THEN
 390   REPEAT STRIP_TAC THEN MATCH_MP_TAC DET_ZERO_ROW THEN EXISTS_TAC `k:num` THEN
 391   ASM_SIMP_TAC[row; LAMBDA_BETA; CART_EQ; VEC_COMPONENT]);;
 392
 393 let BOUNDED_FUNCTIONS_BIJECTIONS_1 = prove
 394  (`!p. p IN {(y,g) | y IN s /\
 395                      g IN {f | (!i. 1 <= i /\ i <= k ==> f i IN s) /\
 396                                (!i. ~(1 <= i /\ i <= k) ==> f i = i)}}
 397        ==> (\(y,g) i. if i = SUC k then y else g(i)) p IN
 398              {f | (!i. 1 <= i /\ i <= SUC k ==> f i IN s) /\
 399                   (!i. ~(1 <= i /\ i <= SUC k) ==> f i = i)} /\
 400            (\h. h(SUC k),(\i. if i = SUC k then i else h(i)))
 401             ((\(y,g) i. if i = SUC k then y else g(i)) p) = p`,
 402   REWRITE_TAC[FORALL_PAIR_THM; IN_ELIM_PAIR_THM] THEN
 403   CONV_TAC(REDEPTH_CONV GEN_BETA_CONV) THEN REWRITE_TAC[IN_ELIM_THM] THEN
 404   MAP_EVERY X_GEN_TAC [`y:num`; `h:num->num`] THEN REPEAT STRIP_TAC THENL
 405    [ASM_MESON_TAC[LE];
 406     ASM_MESON_TAC[LE; ARITH_RULE `~(1 <= i /\ i <= SUC k) ==> ~(i = SUC k)`];
 407     REWRITE_TAC[PAIR_EQ; FUN_EQ_THM] THEN
 408     ASM_MESON_TAC[ARITH_RULE `~(SUC k <= k)`]]);;
 409
 410 let BOUNDED_FUNCTIONS_BIJECTIONS_2 = prove
 411  (`!h. h IN {f | (!i. 1 <= i /\ i <= SUC k ==> f i IN s) /\
 412                  (!i. ~(1 <= i /\ i <= SUC k) ==> f i = i)}
 413        ==> (\h. h(SUC k),(\i. if i = SUC k then i else h(i))) h IN
 414            {(y,g) | y IN s /\
 415                      g IN {f | (!i. 1 <= i /\ i <= k ==> f i IN s) /\
 416                                (!i. ~(1 <= i /\ i <= k) ==> f i = i)}} /\
 417            (\(y,g) i. if i = SUC k then y else g(i))
 418               ((\h. h(SUC k),(\i. if i = SUC k then i else h(i))) h) = h`,
 419   REWRITE_TAC[IN_ELIM_PAIR_THM] THEN
 420   CONV_TAC(REDEPTH_CONV GEN_BETA_CONV) THEN REWRITE_TAC[IN_ELIM_THM] THEN
 421   X_GEN_TAC `h:num->num` THEN REPEAT STRIP_TAC THENL
 422    [FIRST_X_ASSUM MATCH_MP_TAC THEN ARITH_TAC;
 423     ASM_MESON_TAC[ARITH_RULE `i <= k ==> i <= SUC k /\ ~(i = SUC k)`];
 424     ASM_MESON_TAC[ARITH_RULE `i <= SUC k /\ ~(i = SUC k) ==> i <= k`];
 425     REWRITE_TAC[FUN_EQ_THM] THEN ASM_MESON_TAC[LE_REFL]]);;
 426
 427 let FINITE_BOUNDED_FUNCTIONS = prove
 428  (`!s k. FINITE s
 429          ==> FINITE {f | (!i. 1 <= i /\ i <= k ==> f(i) IN s) /\
 430                          (!i. ~(1 <= i /\ i <= k) ==> f(i) = i)}`,
 431   REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN GEN_TAC THEN DISCH_TAC THEN
 432   INDUCT_TAC THENL
 433    [REWRITE_TAC[ARITH_RULE `~(1 <= i /\ i <= 0)`] THEN
 434     SIMP_TAC[GSYM FUN_EQ_THM; SET_RULE `{x | x = y} = {y}`; FINITE_RULES];
 435     ALL_TAC] THEN
 436   UNDISCH_TAC `FINITE(s:num->bool)` THEN POP_ASSUM MP_TAC THEN
 437   REWRITE_TAC[TAUT `a ==> b ==> c <=> b /\ a ==> c`] THEN
 438   DISCH_THEN(MP_TAC o MATCH_MP FINITE_PRODUCT) THEN
 439   DISCH_THEN(MP_TAC o ISPEC `\(y:num,g) i. if i = SUC k then y else g(i)` o
 440                       MATCH_MP FINITE_IMAGE) THEN
 441   MATCH_MP_TAC(TAUT `a = b ==> a ==> b`) THEN AP_TERM_TAC THEN
 442   REWRITE_TAC[EXTENSION; IN_IMAGE] THEN
 443   X_GEN_TAC `h:num->num` THEN EQ_TAC THENL
 444    [STRIP_TAC THEN ASM_SIMP_TAC[BOUNDED_FUNCTIONS_BIJECTIONS_1]; ALL_TAC] THEN
 445   DISCH_TAC THEN EXISTS_TAC
 446     `(\h. h(SUC k),(\i. if i = SUC k then i else h(i))) h` THEN
 447   PURE_ONCE_REWRITE_TAC[CONJ_SYM] THEN CONV_TAC (RAND_CONV SYM_CONV) THEN
 448   MATCH_MP_TAC BOUNDED_FUNCTIONS_BIJECTIONS_2 THEN ASM_REWRITE_TAC[]);;
 449
 450 let DET_LINEAR_ROWS_VSUM_LEMMA = prove
 451  (`!s k a c.
 452          FINITE s /\ k <= dimindex(:N)
 453          ==> det((lambda i. if i <= k then vsum s (a i) else c i):real^N^N) =
 454              sum {f | (!i. 1 <= i /\ i <= k ==> f(i) IN s) /\
 455                       !i. ~(1 <= i /\ i <= k) ==> f(i) = i}
 456                  (\f. det((lambda i. if i <= k then a i (f i) else c i)
 457                           :real^N^N))`,
 458   let lemma = prove
 459    (`(lambda i. if i <= 0 then x(i) else y(i)) = (lambda i. y i)`,
 460     SIMP_TAC[CART_EQ; ARITH; LAMBDA_BETA; ARITH_RULE
 461                  `1 <= k ==> ~(k <= 0)`]) in
 462   ONCE_REWRITE_TAC[IMP_CONJ] THEN
 463   REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN GEN_TAC THEN DISCH_TAC THEN
 464   INDUCT_TAC THENL
 465    [REWRITE_TAC[lemma; LE_0] THEN GEN_TAC THEN
 466     REWRITE_TAC[ARITH_RULE `~(1 <= i /\ i <= 0)`] THEN
 467     REWRITE_TAC[GSYM FUN_EQ_THM; SET_RULE `{x | x = y} = {y}`] THEN
 468     REWRITE_TAC[SUM_SING];
 469     ALL_TAC] THEN
 470   DISCH_TAC THEN FIRST_X_ASSUM(MP_TAC o check (is_imp o concl)) THEN
 471   ASM_SIMP_TAC[ARITH_RULE `SUC k <= n ==> k <= n`] THEN REPEAT STRIP_TAC THEN
 472   GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [LE] THEN
 473   REWRITE_TAC[TAUT
 474    `(if a \/ b then c else d) = (if a then c else if b then c else d)`] THEN
 475   ASM_SIMP_TAC[DET_LINEAR_ROW_VSUM; ARITH_RULE `1 <= SUC k`] THEN
 476   ONCE_REWRITE_TAC[TAUT
 477     `(if a then b else if c then d else e) =
 478      (if c then (if a then b else d) else (if a then b else e))`] THEN
 479   ASM_SIMP_TAC[ARITH_RULE `i <= k ==> ~(i = SUC k)`] THEN
 480   ASM_SIMP_TAC[SUM_SUM_PRODUCT; FINITE_BOUNDED_FUNCTIONS] THEN
 481   MATCH_MP_TAC SUM_EQ_GENERAL_INVERSES THEN
 482   EXISTS_TAC `\(y:num,g) i. if i = SUC k then y else g(i)` THEN
 483   EXISTS_TAC `\h. h(SUC k),(\i. if i = SUC k then i else h(i))` THEN
 484   CONJ_TAC THENL [ACCEPT_TAC BOUNDED_FUNCTIONS_BIJECTIONS_2; ALL_TAC] THEN
 485   X_GEN_TAC `p:num#(num->num)` THEN
 486   DISCH_THEN(STRIP_ASSUME_TAC o MATCH_MP BOUNDED_FUNCTIONS_BIJECTIONS_1) THEN
 487   ASM_REWRITE_TAC[] THEN
 488   SPEC_TAC(`p:num#(num->num)`,`q:num#(num->num)`) THEN
 489   REWRITE_TAC[FORALL_PAIR_THM] THEN
 490   CONV_TAC(ONCE_DEPTH_CONV GEN_BETA_CONV) THEN
 491   MAP_EVERY X_GEN_TAC [`y:num`; `g:num->num`] THEN AP_TERM_TAC THEN
 492   SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 493   REPEAT STRIP_TAC THEN REPEAT(COND_CASES_TAC THEN ASM_REWRITE_TAC[]) THEN
 494   ASM_MESON_TAC[LE; ARITH_RULE `~(SUC k <= k)`]);;
 495
 496 let DET_LINEAR_ROWS_VSUM = prove
 497  (`!s a.
 498          FINITE s
 499          ==> det((lambda i. vsum s (a i)):real^N^N) =
 500              sum {f | (!i. 1 <= i /\ i <= dimindex(:N) ==> f(i) IN s) /\
 501                       !i. ~(1 <= i /\ i <= dimindex(:N)) ==> f(i) = i}
 502                  (\f. det((lambda i. a i (f i)):real^N^N))`,
 503   let lemma = prove
 504    (`(lambda i. if i <= dimindex(:N) then x(i) else y(i)):real^N^N =
 505      (lambda i. x(i))`,
 506     SIMP_TAC[CART_EQ; LAMBDA_BETA]) in
 507   REPEAT STRIP_TAC THEN
 508   MP_TAC(SPECL [`s:num->bool`; `dimindex(:N)`] DET_LINEAR_ROWS_VSUM_LEMMA) THEN
 509   ASM_REWRITE_TAC[LE_REFL; lemma] THEN SIMP_TAC[]);;
 510
 511 let MATRIX_MUL_VSUM_ALT = prove
 512  (`!A:real^N^N B:real^N^N. A ** B =
 513                   lambda i. vsum (1..dimindex(:N)) (\k. A$i$k % B$k)`,
 514   SIMP_TAC[matrix_mul; CART_EQ; LAMBDA_BETA; VECTOR_MUL_COMPONENT;
 515            VSUM_COMPONENT]);;
 516
 517 let DET_ROWS_MUL = prove
 518  (`!a c. det((lambda i. c(i) % a(i)):real^N^N) =
 519          product(1..dimindex(:N)) (\i. c(i)) *
 520          det((lambda i. a(i)):real^N^N)`,
 521   REPEAT GEN_TAC THEN SIMP_TAC[det; LAMBDA_BETA] THEN
 522   SIMP_TAC[GSYM SUM_LMUL; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 523   MATCH_MP_TAC SUM_EQ THEN SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 524   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 525   MATCH_MP_TAC(REAL_RING `b = c * d ==> s * b = c * s * d`) THEN
 526   SIMP_TAC[GSYM PRODUCT_MUL_NUMSEG] THEN
 527   MATCH_MP_TAC PRODUCT_EQ_NUMSEG THEN
 528   ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; VECTOR_MUL_COMPONENT]);;
 529
 530 let DET_MUL = prove
 531  (`!A B:real^N^N. det(A ** B) = det(A) * det(B)`,
 532   REPEAT GEN_TAC THEN REWRITE_TAC[MATRIX_MUL_VSUM_ALT] THEN
 533   SIMP_TAC[DET_LINEAR_ROWS_VSUM; FINITE_NUMSEG] THEN
 534   MATCH_MP_TAC EQ_TRANS THEN
 535   EXISTS_TAC `sum {p | p permutes 1..dimindex(:N)}
 536             (\f. det (lambda i. (A:real^N^N)$i$f i % (B:real^N^N)$f i))` THEN
 537   CONJ_TAC THENL
 538    [REWRITE_TAC[DET_ROWS_MUL] THEN
 539     MATCH_MP_TAC SUM_SUPERSET THEN
 540     SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 541     REWRITE_TAC[SUBSET; IN_ELIM_THM] THEN CONJ_TAC THENL
 542      [MESON_TAC[permutes; IN_NUMSEG]; ALL_TAC] THEN
 543     X_GEN_TAC `f:num->num` THEN REWRITE_TAC[permutes; IN_NUMSEG] THEN
 544     DISCH_THEN(CONJUNCTS_THEN2 STRIP_ASSUME_TAC MP_TAC) THEN
 545     ASM_REWRITE_TAC[] THEN DISCH_TAC THEN
 546     REWRITE_TAC[REAL_ENTIRE] THEN DISJ2_TAC THEN
 547     MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 548     MP_TAC(ISPECL [`1..dimindex(:N)`; `f:num->num`]
 549        SURJECTIVE_IFF_INJECTIVE) THEN
 550     ASM_REWRITE_TAC[SUBSET; IN_NUMSEG; FINITE_NUMSEG; FORALL_IN_IMAGE] THEN
 551     MATCH_MP_TAC(TAUT `(~b ==> c) /\ (b ==> ~a) ==> (a <=> b) ==> c`) THEN
 552     CONJ_TAC THENL
 553      [REWRITE_TAC[NOT_FORALL_THM] THEN
 554       REPEAT(MATCH_MP_TAC MONO_EXISTS THEN GEN_TAC) THEN
 555       SIMP_TAC[CART_EQ; LAMBDA_BETA; row; NOT_IMP];
 556       ALL_TAC] THEN
 557     DISCH_TAC THEN
 558     SUBGOAL_THEN `!x y. (f:num->num)(x) = f(y) ==> x = y` ASSUME_TAC THENL
 559      [REPEAT GEN_TAC THEN
 560       ASM_CASES_TAC `1 <= x /\ x <= dimindex(:N)` THEN
 561       ASM_CASES_TAC `1 <= y /\ y <= dimindex(:N)` THEN
 562       ASM_MESON_TAC[];
 563       ALL_TAC] THEN
 564     ASM_MESON_TAC[];
 565     ALL_TAC] THEN
 566   SIMP_TAC[det; REAL_MUL_SUM; FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 567   MATCH_MP_TAC SUM_EQ THEN SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 568   X_GEN_TAC `p:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 569   FIRST_ASSUM(fun th -> GEN_REWRITE_TAC RAND_CONV
 570     [MATCH_MP SUM_PERMUTATIONS_COMPOSE_R (MATCH_MP PERMUTES_INVERSE th)]) THEN
 571   MATCH_MP_TAC SUM_EQ THEN SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG] THEN
 572   X_GEN_TAC `q:num->num` THEN REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
 573   REWRITE_TAC[o_THM] THEN ONCE_REWRITE_TAC[AC REAL_MUL_AC
 574    `(p * x) * (q * y) = (p * q) * (x * y)`] THEN
 575   BINOP_TAC THENL
 576    [SUBGOAL_THEN `sign(q o inverse p) = sign(p:num->num) * sign(q:num->num)`
 577      (fun t -> SIMP_TAC[REAL_MUL_ASSOC; SIGN_IDEMPOTENT; REAL_MUL_LID; t]) THEN
 578     ASM_MESON_TAC[SIGN_COMPOSE; PERMUTES_INVERSE; PERMUTATION_PERMUTES;
 579                   FINITE_NUMSEG; SIGN_INVERSE; REAL_MUL_SYM];
 580     ALL_TAC] THEN
 581   GEN_REWRITE_TAC (RAND_CONV o RAND_CONV)
 582    [MATCH_MP PRODUCT_PERMUTE_NUMSEG (ASSUME `p permutes 1..dimindex(:N)`)] THEN
 583   SIMP_TAC[GSYM PRODUCT_MUL; FINITE_NUMSEG] THEN
 584   MATCH_MP_TAC PRODUCT_EQ_NUMSEG THEN
 585   ASM_SIMP_TAC[LAMBDA_BETA; LAMBDA_BETA_PERM; o_THM] THEN
 586   X_GEN_TAC `i:num` THEN STRIP_TAC THEN MATCH_MP_TAC EQ_TRANS THEN
 587   EXISTS_TAC `(A:real^N^N)$i$p(i) * (B:real^N^N)$p(i)$q(i)` THEN CONJ_TAC THENL
 588    [ASM_MESON_TAC[VECTOR_MUL_COMPONENT; PERMUTES_IN_IMAGE; IN_NUMSEG];
 589     ASM_MESON_TAC[PERMUTES_INVERSES]]);;
 590
 591 (* ------------------------------------------------------------------------- *)
 592 (* Relation to invertibility.                                                *)
 593 (* ------------------------------------------------------------------------- *)
 594
 595 let INVERTIBLE_DET_NZ = prove
 596  (`!A:real^N^N. invertible(A) <=> ~(det A = &0)`,
 597   GEN_TAC THEN EQ_TAC THENL
 598    [REWRITE_TAC[INVERTIBLE_RIGHT_INVERSE; LEFT_IMP_EXISTS_THM] THEN
 599     GEN_TAC THEN DISCH_THEN(MP_TAC o AP_TERM `det:real^N^N->real`) THEN
 600     REWRITE_TAC[DET_MUL; DET_I] THEN CONV_TAC REAL_RING;
 601     ALL_TAC] THEN
 602   ONCE_REWRITE_TAC[GSYM CONTRAPOS_THM] THEN
 603   REWRITE_TAC[INVERTIBLE_RIGHT_INVERSE] THEN
 604   REWRITE_TAC[MATRIX_RIGHT_INVERTIBLE_INDEPENDENT_ROWS] THEN
 605   REWRITE_TAC[NOT_FORALL_THM; NOT_IMP] THEN
 606   REWRITE_TAC[RIGHT_AND_EXISTS_THM; LEFT_IMP_EXISTS_THM] THEN
 607   MAP_EVERY X_GEN_TAC [`c:num->real`; `i:num`] THEN STRIP_TAC THEN
 608   MP_TAC(SPECL [`A:real^N^N`; `i:num`; `--(row i (A:real^N^N))`]
 609     DET_ROW_SPAN) THEN
 610   ANTS_TAC THENL
 611    [ASM_REWRITE_TAC[] THEN
 612     SUBGOAL_THEN
 613       `--(row i (A:real^N^N)) =
 614        vsum ((1..dimindex(:N)) DELETE i) (\j. inv(c i) % c j % row j A)`
 615     SUBST1_TAC THENL
 616      [ASM_SIMP_TAC[VSUM_DELETE_CASES; FINITE_NUMSEG; IN_NUMSEG; VSUM_LMUL] THEN
 617       ASM_SIMP_TAC[VECTOR_MUL_ASSOC; REAL_MUL_LINV] THEN VECTOR_ARITH_TAC;
 618       ALL_TAC] THEN
 619     MATCH_MP_TAC SPAN_VSUM THEN
 620     REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG; FINITE_DELETE; IN_DELETE] THEN
 621     X_GEN_TAC `j:num` THEN STRIP_TAC THEN REPEAT(MATCH_MP_TAC SPAN_MUL) THEN
 622     MATCH_MP_TAC(CONJUNCT1 SPAN_CLAUSES) THEN
 623     REWRITE_TAC[IN_ELIM_THM] THEN ASM_MESON_TAC[];
 624     ALL_TAC] THEN
 625   DISCH_THEN(SUBST1_TAC o SYM) THEN MATCH_MP_TAC DET_ZERO_ROW THEN
 626   EXISTS_TAC `i:num` THEN
 627   ASM_SIMP_TAC[row; CART_EQ; LAMBDA_BETA; VEC_COMPONENT;
 628                VECTOR_ARITH `x + --x:real^N = vec 0`]);;
 629
 630 let DET_EQ_0 = prove
 631  (`!A:real^N^N. det(A) = &0 <=> ~invertible(A)`,
 632   REWRITE_TAC[INVERTIBLE_DET_NZ]);;
 633
 634 let MATRIX_MUL_LINV = prove
 635  (`!A:real^N^N. ~(det A = &0) ==> matrix_inv A ** A = mat 1`,
 636   SIMP_TAC[MATRIX_INV; DET_EQ_0]);;
 637
 638 let MATRIX_MUL_RINV = prove
 639  (`!A:real^N^N. ~(det A = &0) ==> A ** matrix_inv A = mat 1`,
 640   SIMP_TAC[MATRIX_INV; DET_EQ_0]);;
 641
 642 let DET_MATRIX_EQ_0 = prove
 643  (`!f:real^N->real^N.
 644         linear f
 645         ==> (det(matrix f) = &0 <=>
 646              ~(?g. linear g /\ f o g = I /\ g o f = I))`,
 647   SIMP_TAC[DET_EQ_0; MATRIX_INVERTIBLE]);;
 648
 649 let DET_MATRIX_EQ_0_LEFT = prove
 650  (`!f:real^N->real^N.
 651         linear f
 652         ==> (det(matrix f) = &0 <=>
 653              ~(?g. linear g /\ g o f = I))`,
 654    SIMP_TAC[DET_MATRIX_EQ_0] THEN MESON_TAC[LINEAR_INVERSE_LEFT]);;
 655
 656 let DET_MATRIX_EQ_0_RIGHT = prove
 657  (`!f:real^N->real^N.
 658         linear f
 659         ==> (det(matrix f) = &0 <=>
 660              ~(?g. linear g /\ f o g = I))`,
 661    SIMP_TAC[DET_MATRIX_EQ_0] THEN MESON_TAC[LINEAR_INVERSE_LEFT]);;
 662
 663 let DET_EQ_0_RANK = prove
 664  (`!A:real^N^N. det A = &0 <=> rank A < dimindex(:N)`,
 665   REWRITE_TAC[DET_EQ_0; INVERTIBLE_LEFT_INVERSE; GSYM FULL_RANK_INJECTIVE;
 666               MATRIX_LEFT_INVERTIBLE_INJECTIVE] THEN
 667   GEN_TAC THEN MP_TAC(ISPEC `A:real^N^N` RANK_BOUND) THEN
 668   ARITH_TAC);;
 669
 670 let HOMOGENEOUS_LINEAR_EQUATIONS_DET = prove
 671  (`!A:real^N^N. (?x. ~(x = vec 0) /\ A ** x = vec 0) <=> det A = &0`,
 672   GEN_TAC THEN
 673   REWRITE_TAC[MATRIX_NONFULL_LINEAR_EQUATIONS_EQ; DET_EQ_0_RANK] THEN
 674   MATCH_MP_TAC(ARITH_RULE `r <= MIN N N ==> (~(r = N) <=> r < N)`) THEN
 675   REWRITE_TAC[RANK_BOUND]);;
 676
 677 (* ------------------------------------------------------------------------- *)
 678 (* Cramer's rule.                                                            *)
 679 (* ------------------------------------------------------------------------- *)
 680
 681 let CRAMER_LEMMA_TRANSP = prove
 682  (`!A:real^N^N x:real^N.
 683         1 <= k /\ k <= dimindex(:N)
 684         ==> det((lambda i. if i = k
 685                            then vsum(1..dimindex(:N)) (\i. x$i % row i A)
 686                            else row i A):real^N^N) =
 687             x$k * det A`,
 688   REPEAT STRIP_TAC THEN
 689   SUBGOAL_THEN `1..dimindex(:N) = k INSERT ((1..dimindex(:N)) DELETE k)`
 690   SUBST1_TAC THENL [ASM_MESON_TAC[INSERT_DELETE; IN_NUMSEG]; ALL_TAC] THEN
 691   SIMP_TAC[VSUM_CLAUSES; FINITE_NUMSEG; FINITE_DELETE; IN_DELETE] THEN
 692   REWRITE_TAC[VECTOR_ARITH
 693    `(x:real^N)$k % row k (A:real^N^N) + s =
 694     (x$k - &1) % row k A + row k A + s`] THEN
 695   W(MP_TAC o PART_MATCH (lhs o rand) DET_ROW_ADD o lhand o snd) THEN
 696   ASM_SIMP_TAC[DET_ROW_MUL] THEN DISCH_THEN(K ALL_TAC) THEN
 697   MATCH_MP_TAC(REAL_RING `d = d' /\ e = d' ==> (c - &1) * d + e = c * d'`) THEN
 698   CONJ_TAC THENL
 699    [AP_TERM_TAC THEN ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 700     REPEAT STRIP_TAC THEN COND_CASES_TAC THEN ASM_SIMP_TAC[LAMBDA_BETA; row];
 701     MATCH_MP_TAC DET_ROW_SPAN THEN ASM_REWRITE_TAC[] THEN
 702     ASM_REWRITE_TAC[] THEN MATCH_MP_TAC SPAN_VSUM THEN
 703     REWRITE_TAC[FINITE_NUMSEG; IN_NUMSEG; FINITE_DELETE; IN_DELETE] THEN
 704     REPEAT STRIP_TAC THEN MATCH_MP_TAC SPAN_MUL THEN
 705     MATCH_MP_TAC(CONJUNCT1 SPAN_CLAUSES) THEN
 706     REWRITE_TAC[IN_ELIM_THM] THEN ASM_MESON_TAC[]]);;
 707
 708 let CRAMER_LEMMA = prove
 709  (`!A:real^N^N x:real^N.
 710         1 <= k /\ k <= dimindex(:N)
 711         ==> det((lambda i j. if j = k then (A**x)$i else A$i$j):real^N^N) =
 712             x$k * det(A)`,
 713   REPEAT GEN_TAC THEN DISCH_TAC THEN REWRITE_TAC[MATRIX_MUL_VSUM] THEN
 714   FIRST_ASSUM(MP_TAC o SYM o SPECL [`transp(A:real^N^N)`; `x:real^N`] o
 715               MATCH_MP CRAMER_LEMMA_TRANSP) THEN
 716   REWRITE_TAC[DET_TRANSP] THEN DISCH_THEN SUBST1_TAC THEN
 717   GEN_REWRITE_TAC LAND_CONV [GSYM DET_TRANSP] THEN AP_TERM_TAC THEN
 718   ASM_SIMP_TAC[CART_EQ; transp; LAMBDA_BETA; MATRIX_MUL_VSUM; row; column;
 719         COND_COMPONENT; VECTOR_MUL_COMPONENT; VSUM_COMPONENT]);;
 720
 721 let CRAMER = prove
 722  (`!A:real^N^N x b.
 723         ~(det(A) = &0)
 724         ==> (A ** x = b <=>
 725              x = lambda k.
 726                    det((lambda i j. if j = k then b$i else A$i$j):real^N^N) /
 727                    det(A))`,
 728   GEN_TAC THEN REWRITE_TAC[RIGHT_FORALL_IMP_THM] THEN DISCH_TAC THEN
 729   ONCE_REWRITE_TAC[SWAP_FORALL_THM] THEN GEN_TAC THEN MATCH_MP_TAC(MESON[]
 730    `(?x. p(x)) /\ (!x. p(x) ==> x = a) ==> !x. p(x) <=> x = a`) THEN
 731   CONJ_TAC THENL
 732    [MP_TAC(SPEC `A:real^N^N` INVERTIBLE_DET_NZ) THEN
 733     ASM_MESON_TAC[invertible; MATRIX_VECTOR_MUL_ASSOC; MATRIX_VECTOR_MUL_LID];
 734     GEN_TAC THEN DISCH_THEN(SUBST1_TAC o SYM) THEN
 735     ASM_SIMP_TAC[CART_EQ; CRAMER_LEMMA; LAMBDA_BETA; REAL_FIELD
 736     `~(z = &0) ==> (x = y / z <=> x * z = y)`]]);;
 737
 738 (* ------------------------------------------------------------------------- *)
 739 (* Variants of Cramer's rule for matrix-matrix multiplication.               *)
 740 (* ------------------------------------------------------------------------- *)
 741
 742 let CRAMER_MATRIX_LEFT = prove
 743  (`!A:real^N^N X:real^N^N B:real^N^N.
 744         ~(det A = &0)
 745         ==> (X ** A = B <=>
 746              X = lambda k l.
 747                    det((lambda i j. if j = l then B$k$i else A$j$i):real^N^N) /
 748                    det A)`,
 749   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[CART_EQ] THEN
 750   ASM_SIMP_TAC[MATRIX_MUL_COMPONENT; CRAMER; DET_TRANSP] THEN
 751   SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 752   REPLICATE_TAC 2 (AP_TERM_TAC THEN ABS_TAC THEN AP_TERM_TAC) THEN
 753   AP_TERM_TAC THEN AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 754   SIMP_TAC[CART_EQ; LAMBDA_BETA; transp]);;
 755
 756 let CRAMER_MATRIX_RIGHT = prove
 757  (`!A:real^N^N X:real^N^N B:real^N^N.
 758         ~(det A = &0)
 759         ==> (A ** X = B <=>
 760              X = lambda k l.
 761                    det((lambda i j. if j = k then B$i$l else A$i$j):real^N^N) /
 762                    det A)`,
 763   REPEAT STRIP_TAC THEN
 764   GEN_REWRITE_TAC LAND_CONV [GSYM TRANSP_EQ] THEN
 765   REWRITE_TAC[MATRIX_TRANSP_MUL] THEN
 766   ASM_SIMP_TAC[CRAMER_MATRIX_LEFT; DET_TRANSP] THEN
 767   GEN_REWRITE_TAC LAND_CONV [GSYM TRANSP_EQ] THEN
 768   REWRITE_TAC[TRANSP_TRANSP] THEN AP_TERM_TAC THEN
 769   SIMP_TAC[CART_EQ; LAMBDA_BETA; transp] THEN
 770   REPEAT(GEN_TAC THEN STRIP_TAC) THEN
 771   AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 772   SIMP_TAC[CART_EQ; LAMBDA_BETA; transp]);;
 773
 774 let CRAMER_MATRIX_RIGHT_INVERSE = prove
 775  (`!A:real^N^N A':real^N^N.
 776         A ** A' = mat 1 <=>
 777         ~(det A = &0) /\
 778         A' = lambda k l.
 779                 det((lambda i j. if j = k then if i = l then &1 else &0
 780                                  else A$i$j):real^N^N) /
 781                 det A`,
 782   REPEAT GEN_TAC THEN ASM_CASES_TAC `det(A:real^N^N) = &0` THENL
 783    [ASM_REWRITE_TAC[] THEN
 784     DISCH_THEN(MP_TAC o AP_TERM `det:real^N^N->real`) THEN
 785     ASM_REWRITE_TAC[DET_MUL; DET_I] THEN REAL_ARITH_TAC;
 786     ASM_SIMP_TAC[CRAMER_MATRIX_RIGHT] THEN AP_TERM_TAC THEN
 787     SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 788     REPEAT(GEN_TAC THEN STRIP_TAC) THEN
 789     AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 790     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA; mat]]);;
 791
 792 let CRAMER_MATRIX_LEFT_INVERSE = prove
 793  (`!A:real^N^N A':real^N^N.
 794         A' ** A = mat 1 <=>
 795         ~(det A = &0) /\
 796         A' = lambda k l.
 797                 det((lambda i j. if j = l then if i = k then &1 else &0
 798                                  else A$j$i):real^N^N) /
 799                 det A`,
 800   REPEAT GEN_TAC THEN ASM_CASES_TAC `det(A:real^N^N) = &0` THENL
 801    [ASM_REWRITE_TAC[] THEN
 802     DISCH_THEN(MP_TAC o AP_TERM `det:real^N^N->real`) THEN
 803     ASM_REWRITE_TAC[DET_MUL; DET_I] THEN REAL_ARITH_TAC;
 804     ASM_SIMP_TAC[CRAMER_MATRIX_LEFT] THEN AP_TERM_TAC THEN
 805     SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
 806     REPEAT(GEN_TAC THEN STRIP_TAC) THEN
 807     AP_THM_TAC THEN AP_TERM_TAC THEN AP_TERM_TAC THEN
 808     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA; mat] THEN MESON_TAC[]]);;
 809
 810 (* ------------------------------------------------------------------------- *)
 811 (* Cofactors and their relationship to inverse matrices.                     *)
 812 (* ------------------------------------------------------------------------- *)
 813
 814 let cofactor = new_definition
 815   `(cofactor:real^N^N->real^N^N) A =
 816         lambda i j. det((lambda k l. if k = i /\ l = j then &1
 817                                      else if k = i \/ l = j then &0
 818                                      else A$k$l):real^N^N)`;;
 819
 820 let COFACTOR_TRANSP = prove
 821  (`!A:real^N^N. cofactor(transp A) = transp(cofactor A)`,
 822   SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA; transp] THEN REPEAT STRIP_TAC THEN
 823   GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN
 824   AP_TERM_TAC THEN SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA; transp] THEN
 825   MESON_TAC[]);;
 826
 827 let COFACTOR_COLUMN = prove
 828  (`!A:real^N^N.
 829         cofactor A =
 830         lambda i j. det((lambda k l. if l = j then if k = i then &1 else &0
 831                                      else A$k$l):real^N^N)`,
 832   GEN_TAC THEN CONV_TAC SYM_CONV THEN
 833   SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA] THEN
 834   X_GEN_TAC `i:num` THEN STRIP_TAC THEN
 835   X_GEN_TAC `j:num` THEN STRIP_TAC THEN
 836   REWRITE_TAC[det] THEN MATCH_MP_TAC SUM_EQ THEN
 837   REWRITE_TAC[FORALL_IN_GSPEC] THEN GEN_TAC THEN
 838   DISCH_TAC THEN AP_TERM_TAC THEN
 839   ASM_CASES_TAC `(p:num->num) i = j` THENL
 840    [MATCH_MP_TAC PRODUCT_EQ THEN
 841     X_GEN_TAC `k:num` THEN SIMP_TAC[IN_NUMSEG; LAMBDA_BETA] THEN STRIP_TAC THEN
 842     SUBGOAL_THEN `(p:num->num) k IN 1..dimindex(:N)` MP_TAC THENL
 843      [ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG];
 844       SIMP_TAC[LAMBDA_BETA; IN_NUMSEG] THEN STRIP_TAC] THEN
 845     ASM_CASES_TAC `(p:num->num) k = j` THEN ASM_REWRITE_TAC[] THEN
 846     COND_CASES_TAC THEN ASM_REWRITE_TAC[] THEN ASM_MESON_TAC[];
 847     MATCH_MP_TAC(REAL_ARITH `s = &0 /\ t = &0 ==> s = t`) THEN
 848     ASM_SIMP_TAC[PRODUCT_EQ_0; FINITE_NUMSEG] THEN CONJ_TAC THEN
 849     EXISTS_TAC `inverse (p:num->num) j` THEN
 850     ASM_SIMP_TAC[IN_NUMSEG; LAMBDA_BETA] THEN
 851     (SUBGOAL_THEN `inverse(p:num->num) j IN 1..dimindex(:N)` MP_TAC THENL
 852       [ASM_MESON_TAC[PERMUTES_IN_IMAGE; PERMUTES_INVERSE; IN_NUMSEG];
 853        SIMP_TAC[LAMBDA_BETA; IN_NUMSEG] THEN STRIP_TAC] THEN
 854      SUBGOAL_THEN `(p:num->num)(inverse p j) = j` SUBST1_TAC THENL
 855       [ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG];
 856        ASM_SIMP_TAC[LAMBDA_BETA] THEN
 857         ASM_MESON_TAC[PERMUTES_INVERSE_EQ]])]);;
 858
 859 let COFACTOR_ROW = prove
 860  (`!A:real^N^N.
 861         cofactor A =
 862         lambda i j. det((lambda k l. if k = i then if l = j then &1 else &0
 863                                      else A$k$l):real^N^N)`,
 864   GEN_TAC THEN ONCE_REWRITE_TAC[GSYM TRANSP_EQ] THEN
 865   REWRITE_TAC[GSYM COFACTOR_TRANSP] THEN
 866   SIMP_TAC[COFACTOR_COLUMN; CART_EQ; LAMBDA_BETA; transp] THEN
 867   REPEAT STRIP_TAC THEN
 868   GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN
 869   AP_TERM_TAC THEN SIMP_TAC[cofactor; CART_EQ; LAMBDA_BETA; transp]);;
 870
 871 let MATRIX_RIGHT_INVERSE_COFACTOR = prove
 872  (`!A:real^N^N A':real^N^N.
 873         A ** A' = mat 1 <=>
 874         ~(det A = &0) /\ A' = inv(det A) %% transp(cofactor A)`,
 875   REPEAT GEN_TAC THEN REWRITE_TAC[CRAMER_MATRIX_RIGHT_INVERSE] THEN
 876   ASM_CASES_TAC `det(A:real^N^N) = &0` THEN ASM_REWRITE_TAC[] THEN
 877   AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA; MATRIX_CMUL_COMPONENT] THEN
 878   X_GEN_TAC `k:num` THEN STRIP_TAC THEN
 879   X_GEN_TAC `l:num` THEN STRIP_TAC THEN
 880   REWRITE_TAC[ONCE_REWRITE_RULE[REAL_MUL_SYM] real_div] THEN AP_TERM_TAC THEN
 881   ASM_SIMP_TAC[transp; COFACTOR_COLUMN; LAMBDA_BETA] THEN
 882   AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA]);;
 883
 884 let MATRIX_LEFT_INVERSE_COFACTOR = prove
 885  (`!A:real^N^N A':real^N^N.
 886         A' ** A = mat 1 <=>
 887         ~(det A = &0) /\ A' = inv(det A) %% transp(cofactor A)`,
 888   REPEAT GEN_TAC THEN
 889   ONCE_REWRITE_TAC[MATRIX_LEFT_RIGHT_INVERSE] THEN
 890   REWRITE_TAC[MATRIX_RIGHT_INVERSE_COFACTOR]);;
 891
 892 let MATRIX_INV_COFACTOR = prove
 893  (`!A. ~(det A = &0) ==> matrix_inv A = inv(det A) %% transp(cofactor A)`,
 894   GEN_TAC THEN DISCH_THEN(MP_TAC o MATCH_MP MATRIX_MUL_LINV) THEN
 895   SIMP_TAC[MATRIX_LEFT_INVERSE_COFACTOR]);;
 896
 897 let COFACTOR_MATRIX_INV = prove
 898  (`!A:real^N^N. ~(det A = &0) ==> cofactor A = det(A) %% transp(matrix_inv A)`,
 899   SIMP_TAC[MATRIX_INV_COFACTOR; TRANSP_MATRIX_CMUL; TRANSP_TRANSP] THEN
 900   SIMP_TAC[MATRIX_CMUL_ASSOC; REAL_MUL_RINV; MATRIX_CMUL_LID]);;
 901
 902 let DET_COFACTOR_EXPANSION = prove
 903  (`!A:real^N^N i.
 904         1 <= i /\ i <= dimindex(:N)
 905         ==> det A = sum (1..dimindex(:N))
 906                         (\j. A$i$j * (cofactor A)$i$j)`,
 907   REPEAT STRIP_TAC THEN ASM_SIMP_TAC[COFACTOR_COLUMN; LAMBDA_BETA; det] THEN
 908   REWRITE_TAC[GSYM SUM_LMUL] THEN
 909   W(MP_TAC o PART_MATCH (lhand o rand) SUM_SWAP o rand o snd) THEN
 910   ANTS_TAC THENL [SIMP_TAC[FINITE_PERMUTATIONS; FINITE_NUMSEG]; ALL_TAC] THEN
 911   DISCH_THEN SUBST1_TAC THEN
 912   MATCH_MP_TAC SUM_EQ THEN REWRITE_TAC[FORALL_IN_GSPEC] THEN
 913   GEN_TAC THEN DISCH_TAC THEN
 914   ONCE_REWRITE_TAC[REAL_ARITH `a * s * p:real = s * a * p`] THEN
 915   REWRITE_TAC[SUM_LMUL] THEN AP_TERM_TAC THEN MATCH_MP_TAC EQ_TRANS THEN
 916   EXISTS_TAC
 917    `sum (1..dimindex (:N))
 918         (\j. (A:real^N^N)$i$j *
 919              product
 920               (inverse p j INSERT ((1..dimindex(:N)) DELETE (inverse p j)))
 921               (\k. if k = inverse p j then if k = i then &1 else &0
 922                    else A$k$(p k)))` THEN
 923   CONJ_TAC THENL
 924    [SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_PERMUTATIONS;
 925              FINITE_NUMSEG; IN_DELETE] THEN
 926     SUBGOAL_THEN `!j. inverse (p:num->num) j = i <=> j = p i`
 927      (fun th -> REWRITE_TAC[th])
 928     THENL [ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG]; ALL_TAC] THEN
 929     REWRITE_TAC[REAL_ARITH
 930      `x * (if p then &1 else &0) * y = if p then x * y else &0`] THEN
 931     SIMP_TAC[SUM_DELTA] THEN COND_CASES_TAC THENL
 932      [ALL_TAC; ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]] THEN
 933     SUBGOAL_THEN
 934      `1..dimindex(:N) = i INSERT ((1..dimindex(:N)) DELETE i)`
 935      (fun th -> GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) [th])
 936     THENL
 937      [ASM_SIMP_TAC[IN_NUMSEG; SET_RULE `s = x INSERT (s DELETE x) <=> x IN s`];
 938       SIMP_TAC[PRODUCT_CLAUSES; FINITE_DELETE; FINITE_NUMSEG; IN_DELETE] THEN
 939       AP_TERM_TAC THEN MATCH_MP_TAC(MESON[PRODUCT_EQ]
 940        `s = t /\ (!x. x IN t ==> f x = g x) ==> product s f = product t g`) THEN
 941       SIMP_TAC[IN_DELETE] THEN ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG]];
 942     MATCH_MP_TAC SUM_EQ_NUMSEG THEN X_GEN_TAC `j:num` THEN STRIP_TAC THEN
 943     REWRITE_TAC[] THEN AP_TERM_TAC THEN MATCH_MP_TAC(MESON[PRODUCT_EQ]
 944      `s = t /\ (!x. x IN t ==> f x = g x) ==> product s f = product t g`) THEN
 945     CONJ_TAC THENL
 946      [REWRITE_TAC[SET_RULE `x INSERT (s DELETE x) = s <=> x IN s`] THEN
 947       ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG; PERMUTES_INVERSE];
 948       X_GEN_TAC `k:num` THEN REWRITE_TAC[IN_NUMSEG] THEN STRIP_TAC THEN
 949       SUBGOAL_THEN `(p:num->num) k IN 1..dimindex(:N)` MP_TAC THENL
 950        [ASM_MESON_TAC[PERMUTES_IN_IMAGE; IN_NUMSEG]; ALL_TAC] THEN
 951       SIMP_TAC[LAMBDA_BETA; IN_NUMSEG] THEN
 952       ASM_MESON_TAC[PERMUTES_INVERSES; IN_NUMSEG]]]);;
 953
 954 let MATRIX_MUL_RIGHT_COFACTOR = prove
 955  (`!A:real^N^N. A ** transp(cofactor A) = det(A) %% mat 1`,
 956   GEN_TAC THEN
 957   SIMP_TAC[CART_EQ; MATRIX_CMUL_COMPONENT; mat;
 958            matrix_mul; LAMBDA_BETA; transp] THEN
 959   X_GEN_TAC `i:num` THEN STRIP_TAC THEN
 960   X_GEN_TAC `i':num` THEN STRIP_TAC THEN
 961   COND_CASES_TAC THEN
 962   ASM_SIMP_TAC[GSYM DET_COFACTOR_EXPANSION; REAL_MUL_RID] THEN
 963   MATCH_MP_TAC EQ_TRANS THEN
 964   EXISTS_TAC `det((lambda k l. if k = i' then (A:real^N^N)$i$l
 965                                else A$k$l):real^N^N)` THEN
 966   CONJ_TAC THENL
 967    [MP_TAC(GEN `A:real^N^N`
 968      (ISPECL [`A:real^N^N`; `i':num`] DET_COFACTOR_EXPANSION)) THEN
 969     ASM_REWRITE_TAC[] THEN ASM_SIMP_TAC[] THEN DISCH_THEN(K ALL_TAC) THEN
 970     MATCH_MP_TAC SUM_EQ THEN X_GEN_TAC `j:num` THEN
 971     REWRITE_TAC[IN_NUMSEG] THEN STRIP_TAC THEN
 972     ASM_SIMP_TAC[LAMBDA_BETA] THEN AP_TERM_TAC THEN
 973     ASM_SIMP_TAC[cofactor; LAMBDA_BETA] THEN AP_TERM_TAC THEN
 974     SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN ASM_MESON_TAC[];
 975     REWRITE_TAC[REAL_MUL_RZERO] THEN MATCH_MP_TAC DET_IDENTICAL_ROWS THEN
 976     MAP_EVERY EXISTS_TAC [`i:num`;` i':num`] THEN
 977     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA; row]]);;
 978
 979 let MATRIX_MUL_LEFT_COFACTOR = prove
 980  (`!A:real^N^N. transp(cofactor A) ** A = det(A) %% mat 1`,
 981   GEN_TAC THEN ONCE_REWRITE_TAC[GSYM TRANSP_EQ] THEN
 982   REWRITE_TAC[MATRIX_TRANSP_MUL] THEN
 983   ONCE_REWRITE_TAC[GSYM COFACTOR_TRANSP] THEN
 984   REWRITE_TAC[MATRIX_MUL_RIGHT_COFACTOR; TRANSP_MATRIX_CMUL] THEN
 985   REWRITE_TAC[DET_TRANSP; TRANSP_MAT]);;
 986
 987 (* ------------------------------------------------------------------------- *)
 988 (* Explicit formulas for low dimensions.                                     *)
 989 (* ------------------------------------------------------------------------- *)
 990
 991 let PRODUCT_1 = prove
 992  (`product(1..1) f = f(1)`,
 993   REWRITE_TAC[PRODUCT_SING_NUMSEG]);;
 994
 995 let PRODUCT_2 = prove
 996  (`!t. product(1..2) t = t(1) * t(2)`,
 997   REWRITE_TAC[num_CONV `2`; PRODUCT_CLAUSES_NUMSEG] THEN
 998   REWRITE_TAC[PRODUCT_SING_NUMSEG; ARITH; REAL_MUL_ASSOC]);;
 999
1000 let PRODUCT_3 = prove
1001  (`!t. product(1..3) t = t(1) * t(2) * t(3)`,
1002   REWRITE_TAC[num_CONV `3`; num_CONV `2`; PRODUCT_CLAUSES_NUMSEG] THEN
1003   REWRITE_TAC[PRODUCT_SING_NUMSEG; ARITH; REAL_MUL_ASSOC]);;
1004
1005 let DET_1 = prove
1006  (`!A:real^1^1. det A = A$1$1`,
1007   REWRITE_TAC[det; DIMINDEX_1; PERMUTES_SING; NUMSEG_SING] THEN
1008   REWRITE_TAC[SUM_SING; SET_RULE `{x | x = a} = {a}`; PRODUCT_SING] THEN
1009   REWRITE_TAC[SIGN_I; I_THM] THEN REAL_ARITH_TAC);;
1010
1011 let DET_2 = prove
1012  (`!A:real^2^2. det A = A$1$1 * A$2$2 - A$1$2 * A$2$1`,
1013   GEN_TAC THEN REWRITE_TAC[det; DIMINDEX_2] THEN
1014   CONV_TAC(LAND_CONV(RATOR_CONV(ONCE_DEPTH_CONV NUMSEG_CONV))) THEN
1015   SIMP_TAC[SUM_OVER_PERMUTATIONS_INSERT; FINITE_INSERT; FINITE_EMPTY;
1016            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1017   REWRITE_TAC[PERMUTES_EMPTY; SUM_SING; SET_RULE `{x | x = a} = {a}`] THEN
1018   REWRITE_TAC[SWAP_REFL; I_O_ID] THEN
1019   REWRITE_TAC[GSYM(NUMSEG_CONV `1..2`); SUM_2] THEN
1020   SIMP_TAC[SUM_CLAUSES; FINITE_INSERT; FINITE_EMPTY;
1021            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1022   SIMP_TAC[SIGN_COMPOSE; PERMUTATION_SWAP] THEN
1023   REWRITE_TAC[SIGN_SWAP; ARITH] THEN REWRITE_TAC[PRODUCT_2] THEN
1024   REWRITE_TAC[o_THM; swap; ARITH] THEN REAL_ARITH_TAC);;
1025
1026 let DET_3 = prove
1027  (`!A:real^3^3.
1028         det(A) = A$1$1 * A$2$2 * A$3$3 +
1029                  A$1$2 * A$2$3 * A$3$1 +
1030                  A$1$3 * A$2$1 * A$3$2 -
1031                  A$1$1 * A$2$3 * A$3$2 -
1032                  A$1$2 * A$2$1 * A$3$3 -
1033                  A$1$3 * A$2$2 * A$3$1`,
1034   GEN_TAC THEN REWRITE_TAC[det; DIMINDEX_3] THEN
1035   CONV_TAC(LAND_CONV(RATOR_CONV(ONCE_DEPTH_CONV NUMSEG_CONV))) THEN
1036   SIMP_TAC[SUM_OVER_PERMUTATIONS_INSERT; FINITE_INSERT; FINITE_EMPTY;
1037            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1038   REWRITE_TAC[PERMUTES_EMPTY; SUM_SING; SET_RULE `{x | x = a} = {a}`] THEN
1039   REWRITE_TAC[SWAP_REFL; I_O_ID] THEN
1040   REWRITE_TAC[GSYM(NUMSEG_CONV `1..3`); SUM_3] THEN
1041   SIMP_TAC[SUM_CLAUSES; FINITE_INSERT; FINITE_EMPTY;
1042            ARITH_EQ; IN_INSERT; NOT_IN_EMPTY] THEN
1043   SIMP_TAC[SIGN_COMPOSE; PERMUTATION_SWAP] THEN
1044   REWRITE_TAC[SIGN_SWAP; ARITH] THEN REWRITE_TAC[PRODUCT_3] THEN
1045   REWRITE_TAC[o_THM; swap; ARITH] THEN REAL_ARITH_TAC);;
1046
1047 (* ------------------------------------------------------------------------- *)
1048 (* Grassmann-Plucker relations for n = 2 and n = 3.                          *)
1049 (* I have a proof of the general n case but the proof is a bit long and the  *)
1050 (* result doesn't seem generally useful enough to go in the main theories.   *)
1051 (* ------------------------------------------------------------------------- *)
1052
1053 let GRASSMANN_PLUCKER_2 = prove
1054  (`!x1 x2 y1 y2:real^2.
1055         det(vector[x1;x2]) * det(vector[y1;y2]) =
1056           det(vector[y1;x2]) * det(vector[x1;y2]) +
1057           det(vector[y2;x2]) * det(vector[y1;x1])`,
1058   REWRITE_TAC[DET_2; VECTOR_2] THEN REAL_ARITH_TAC);;
1059
1060 let GRASSMANN_PLUCKER_3 = prove
1061  (`!x1 x2 x3 y1 y2 y3:real^3.
1062         det(vector[x1;x2;x3]) * det(vector[y1;y2;y3]) =
1063           det(vector[y1;x2;x3]) * det(vector[x1;y2;y3]) +
1064           det(vector[y2;x2;x3]) * det(vector[y1;x1;y3]) +
1065           det(vector[y3;x2;x3]) * det(vector[y1;y2;x1])`,
1066   REWRITE_TAC[DET_3; VECTOR_3] THEN REAL_ARITH_TAC);;
1067
1068 (* ------------------------------------------------------------------------- *)
1069 (* Determinants of integer matrices.                                         *)
1070 (* ------------------------------------------------------------------------- *)
1071
1072 let INTEGER_PRODUCT = prove
1073  (`!f s. (!x. x IN s ==> integer(f x)) ==> integer(product s f)`,
1074   REPEAT STRIP_TAC THEN MATCH_MP_TAC PRODUCT_CLOSED THEN
1075   ASM_REWRITE_TAC[INTEGER_CLOSED]);;
1076
1077 let INTEGER_SIGN = prove
1078  (`!p. integer(sign p)`,
1079   SIMP_TAC[sign; COND_RAND; INTEGER_CLOSED; COND_ID]);;
1080
1081 let INTEGER_DET = prove
1082  (`!M:real^N^N.
1083         (!i j. 1 <= i /\ i <= dimindex(:N) /\
1084                1 <= j /\ j <= dimindex(:N)
1085                ==> integer(M$i$j))
1086         ==> integer(det M)`,
1087   REPEAT STRIP_TAC THEN REWRITE_TAC[det] THEN
1088   MATCH_MP_TAC INTEGER_SUM THEN X_GEN_TAC `p:num->num` THEN
1089   REWRITE_TAC[IN_ELIM_THM] THEN DISCH_TAC THEN
1090   MATCH_MP_TAC INTEGER_MUL THEN REWRITE_TAC[INTEGER_SIGN] THEN
1091   MATCH_MP_TAC INTEGER_PRODUCT THEN REWRITE_TAC[IN_NUMSEG] THEN
1092   REPEAT STRIP_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN
1093   ASM_MESON_TAC[IN_NUMSEG; permutes]);;
1094
1095 (* ------------------------------------------------------------------------- *)
1096 (* Orthogonality of a transformation and matrix.                             *)
1097 (* ------------------------------------------------------------------------- *)
1098
1099 let orthogonal_transformation = new_definition
1100  `orthogonal_transformation(f:real^N->real^N) <=>
1101         linear f /\ !v w. f(v) dot f(w) = v dot w`;;
1102
1103 let ORTHOGONAL_TRANSFORMATION = prove
1104  (`!f. orthogonal_transformation f <=> linear f /\ !v. norm(f v) = norm(v)`,
1105   GEN_TAC THEN REWRITE_TAC[orthogonal_transformation] THEN EQ_TAC THENL
1106    [MESON_TAC[vector_norm]; SIMP_TAC[DOT_NORM] THEN MESON_TAC[LINEAR_ADD]]);;
1107
1108 let orthogonal_matrix = new_definition
1109  `orthogonal_matrix(Q:real^N^N) <=>
1110       transp(Q) ** Q = mat 1 /\ Q ** transp(Q) = mat 1`;;
1111
1112 let ORTHOGONAL_MATRIX = prove
1113  (`orthogonal_matrix(Q:real^N^N) <=> transp(Q) ** Q = mat 1`,
1114   MESON_TAC[MATRIX_LEFT_RIGHT_INVERSE; orthogonal_matrix]);;
1115
1116 let ORTHOGONAL_MATRIX_ID = prove
1117  (`orthogonal_matrix(mat 1)`,
1118   REWRITE_TAC[orthogonal_matrix; TRANSP_MAT; MATRIX_MUL_LID]);;
1119
1120 let ORTHOGONAL_MATRIX_MUL = prove
1121  (`!A B. orthogonal_matrix A /\ orthogonal_matrix B
1122          ==> orthogonal_matrix(A ** B)`,
1123   REWRITE_TAC[orthogonal_matrix; MATRIX_TRANSP_MUL] THEN
1124   REPEAT STRIP_TAC THEN ONCE_REWRITE_TAC[GSYM MATRIX_MUL_ASSOC] THEN
1125   GEN_REWRITE_TAC (LAND_CONV o RAND_CONV) [MATRIX_MUL_ASSOC] THEN
1126   ASM_REWRITE_TAC[MATRIX_MUL_LID; MATRIX_MUL_RID]);;
1127
1128 let ORTHOGONAL_TRANSFORMATION_MATRIX = prove
1129  (`!f:real^N->real^N.
1130      orthogonal_transformation f <=> linear f /\ orthogonal_matrix(matrix f)`,
1131   REPEAT STRIP_TAC THEN EQ_TAC THENL
1132    [REWRITE_TAC[orthogonal_transformation; ORTHOGONAL_MATRIX] THEN
1133     STRIP_TAC THEN ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN
1134     X_GEN_TAC `i:num` THEN STRIP_TAC THEN
1135     X_GEN_TAC `j:num` THEN STRIP_TAC THEN
1136     FIRST_X_ASSUM(MP_TAC o SPECL [`basis i:real^N`; `basis j:real^N`]) THEN
1137     FIRST_ASSUM(fun th -> REWRITE_TAC[GSYM(MATCH_MP MATRIX_WORKS th)]) THEN
1138     REWRITE_TAC[DOT_MATRIX_VECTOR_MUL] THEN
1139     ABBREV_TAC `A = transp (matrix f) ** matrix(f:real^N->real^N)` THEN
1140     ASM_SIMP_TAC[matrix_mul; columnvector; rowvector; basis; LAMBDA_BETA;
1141              SUM_DELTA; DIMINDEX_1; LE_REFL; dot; IN_NUMSEG; mat;
1142              MESON[REAL_MUL_LID; REAL_MUL_LZERO; REAL_MUL_RID; REAL_MUL_RZERO]
1143               `(if b then &1 else &0) * x = (if b then x else &0) /\
1144                x * (if b then &1 else &0) = (if b then x else &0)`];
1145     REWRITE_TAC[orthogonal_matrix; ORTHOGONAL_TRANSFORMATION; NORM_EQ] THEN
1146     DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
1147     FIRST_ASSUM(fun th -> REWRITE_TAC[GSYM(MATCH_MP MATRIX_WORKS th)]) THEN
1148     ASM_REWRITE_TAC[DOT_MATRIX_VECTOR_MUL] THEN
1149     SIMP_TAC[DOT_MATRIX_PRODUCT; MATRIX_MUL_LID]]);;
1150
1151 let DET_ORTHOGONAL_MATRIX = prove
1152  (`!Q. orthogonal_matrix Q ==> det(Q) = &1 \/ det(Q) = -- &1`,
1153   GEN_TAC THEN REWRITE_TAC[REAL_RING `x = &1 \/ x = -- &1 <=> x * x = &1`] THEN
1154   GEN_REWRITE_TAC (RAND_CONV o LAND_CONV o RAND_CONV) [GSYM DET_TRANSP] THEN
1155   SIMP_TAC[GSYM DET_MUL; orthogonal_matrix; DET_I]);;
1156
1157 let ORTHOGONAL_MATRIX_TRANSP = prove
1158  (`!A:real^N^N. orthogonal_matrix(transp A) <=> orthogonal_matrix A`,
1159   REWRITE_TAC[orthogonal_matrix; TRANSP_TRANSP; CONJ_ACI]);;
1160
1161 let MATRIX_MUL_LTRANSP_DOT_COLUMN = prove
1162  (`!A:real^N^N. transp A ** A = (lambda i j. (column i A) dot (column j A))`,
1163   SIMP_TAC[matrix_mul; CART_EQ; LAMBDA_BETA; transp; dot; column]);;
1164
1165 let MATRIX_MUL_RTRANSP_DOT_ROW = prove
1166  (`!A:real^N^N. A ** transp A = (lambda i j. (row i A) dot (row j A))`,
1167   SIMP_TAC[matrix_mul; CART_EQ; LAMBDA_BETA; transp; dot; row]);;
1168
1169 let ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS = prove
1170  (`!A:real^N^N.
1171         orthogonal_matrix A <=>
1172         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(column i A) = &1) /\
1173         (!i j. 1 <= i /\ i <= dimindex(:N) /\
1174                1 <= j /\ j <= dimindex(:N) /\ ~(i = j)
1175                ==> orthogonal (column i A) (column j A))`,
1176   REWRITE_TAC[ORTHOGONAL_MATRIX] THEN
1177   SIMP_TAC[MATRIX_MUL_LTRANSP_DOT_COLUMN; CART_EQ; mat; LAMBDA_BETA] THEN
1178   REWRITE_TAC[orthogonal; NORM_EQ_1] THEN MESON_TAC[]);;
1179
1180 let ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS = prove
1181  (`!A:real^N^N.
1182         orthogonal_matrix A <=>
1183         (!i. 1 <= i /\ i <= dimindex(:N) ==> norm(row i A) = &1) /\
1184         (!i j. 1 <= i /\ i <= dimindex(:N) /\
1185                1 <= j /\ j <= dimindex(:N) /\ ~(i = j)
1186                ==> orthogonal (row i A) (row j A))`,
1187   ONCE_REWRITE_TAC[GSYM ORTHOGONAL_MATRIX_TRANSP] THEN
1188   SIMP_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS; COLUMN_TRANSP]);;
1189
1190 let ORTHOGONAL_MATRIX_2 = prove
1191  (`!A:real^2^2. orthogonal_matrix A <=>
1192                 A$1$1 pow 2 + A$2$1 pow 2 = &1 /\
1193                 A$1$2 pow 2 + A$2$2 pow 2 = &1 /\
1194                 A$1$1 * A$1$2 + A$2$1 * A$2$2 = &0`,
1195   SIMP_TAC[orthogonal_matrix; CART_EQ; matrix_mul; LAMBDA_BETA;
1196            TRANSP_COMPONENT; MAT_COMPONENT] THEN
1197   REWRITE_TAC[DIMINDEX_2; FORALL_2; SUM_2] THEN
1198   CONV_TAC NUM_REDUCE_CONV THEN CONV_TAC REAL_RING);;
1199
1200 let ORTHOGONAL_MATRIX_2_ALT = prove
1201  (`!A:real^2^2. orthogonal_matrix A <=>
1202                 A$1$1 pow 2 + A$2$1 pow 2 = &1 /\
1203                 (A$1$1 = A$2$2 /\ A$1$2 = --(A$2$1) \/
1204                  A$1$1 = --(A$2$2) /\ A$1$2 = A$2$1)`,
1205   REWRITE_TAC[ORTHOGONAL_MATRIX_2] THEN CONV_TAC REAL_RING);;
1206
1207 (* ------------------------------------------------------------------------- *)
1208 (* Linearity of scaling, and hence isometry, that preserves origin.          *)
1209 (* ------------------------------------------------------------------------- *)
1210
1211 let SCALING_LINEAR = prove
1212  (`!f:real^M->real^N c.
1213         (f(vec 0) = vec 0) /\ (!x y. dist(f x,f y) = c * dist(x,y))
1214         ==> linear(f)`,
1215   REPEAT STRIP_TAC THEN
1216   SUBGOAL_THEN `!v w. ((f:real^M->real^N) v) dot (f w) = c pow 2 * (v dot w)`
1217   ASSUME_TAC THENL
1218    [FIRST_ASSUM(MP_TAC o GEN `v:real^M` o
1219       SPECL [`v:real^M`; `vec 0 :real^M`]) THEN
1220     REWRITE_TAC[dist] THEN ASM_REWRITE_TAC[VECTOR_SUB_RZERO] THEN
1221     DISCH_TAC THEN ASM_REWRITE_TAC[DOT_NORM_NEG; GSYM dist] THEN
1222     REAL_ARITH_TAC;
1223     ALL_TAC] THEN
1224   REWRITE_TAC[linear; VECTOR_EQ] THEN
1225   ASM_REWRITE_TAC[DOT_LADD; DOT_RADD; DOT_LMUL; DOT_RMUL] THEN
1226   REAL_ARITH_TAC);;
1227
1228 let ISOMETRY_LINEAR = prove
1229  (`!f:real^M->real^N.
1230         (f(vec 0) = vec 0) /\ (!x y. dist(f x,f y) = dist(x,y))
1231         ==> linear(f)`,
1232   MESON_TAC[SCALING_LINEAR; REAL_MUL_LID]);;
1233
1234 (* ------------------------------------------------------------------------- *)
1235 (* Hence another formulation of orthogonal transformation.                   *)
1236 (* ------------------------------------------------------------------------- *)
1237
1238 let ORTHOGONAL_TRANSFORMATION_ISOMETRY = prove
1239  (`!f:real^N->real^N.
1240         orthogonal_transformation f <=>
1241         (f(vec 0) = vec 0) /\ (!x y. dist(f x,f y) = dist(x,y))`,
1242   GEN_TAC THEN REWRITE_TAC[ORTHOGONAL_TRANSFORMATION] THEN EQ_TAC THENL
1243    [MESON_TAC[LINEAR_0; LINEAR_SUB; dist]; STRIP_TAC] THEN
1244   ASM_SIMP_TAC[ISOMETRY_LINEAR] THEN X_GEN_TAC `x:real^N` THEN
1245   FIRST_X_ASSUM(MP_TAC o SPECL [`x:real^N`; `vec 0:real^N`]) THEN
1246   ASM_REWRITE_TAC[dist; VECTOR_SUB_RZERO]);;
1247
1248 (* ------------------------------------------------------------------------- *)
1249 (* Can extend an isometry from unit sphere.                                  *)
1250 (* ------------------------------------------------------------------------- *)
1251
1252 let ISOMETRY_SPHERE_EXTEND = prove
1253  (`!f:real^N->real^N.
1254         (!x. norm(x) = &1 ==> norm(f x) = &1) /\
1255         (!x y. norm(x) = &1 /\ norm(y) = &1 ==> dist(f x,f y) = dist(x,y))
1256         ==> ?g. orthogonal_transformation g /\
1257                 (!x. norm(x) = &1 ==> g(x) = f(x))`,
1258   let lemma = prove
1259    (`!x:real^N y:real^N x':real^N y':real^N x0 y0 x0' y0'.
1260           x = norm(x) % x0 /\ y = norm(y) % y0 /\
1261           x' = norm(x) % x0' /\ y' = norm(y) % y0' /\
1262           norm(x0) = &1 /\ norm(x0') = &1 /\ norm(y0) = &1 /\ norm(y0') = &1 /\
1263           norm(x0' - y0') = norm(x0 - y0)
1264           ==> norm(x' - y') = norm(x - y)`,
1265     REPEAT GEN_TAC THEN
1266     MAP_EVERY ABBREV_TAC [`a = norm(x:real^N)`; `b = norm(y:real^N)`] THEN
1267     REPLICATE_TAC 4 (DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC)) THEN
1268     ASM_REWRITE_TAC[] THEN REWRITE_TAC[NORM_EQ; NORM_EQ_1] THEN
1269     REWRITE_TAC[DOT_LSUB; DOT_RSUB; DOT_LMUL; DOT_RMUL] THEN
1270     REWRITE_TAC[DOT_SYM] THEN CONV_TAC REAL_RING) in
1271   REPEAT STRIP_TAC THEN
1272   EXISTS_TAC `\x. if x = vec 0 then vec 0
1273                   else norm(x) % (f:real^N->real^N)(inv(norm x) % x)` THEN
1274   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_ISOMETRY] THEN
1275   SIMP_TAC[VECTOR_MUL_LID; REAL_INV_1] THEN CONJ_TAC THENL
1276    [ALL_TAC; MESON_TAC[NORM_0; REAL_ARITH `~(&1 = &0)`]] THEN
1277   REPEAT GEN_TAC THEN REPEAT(COND_CASES_TAC THEN ASM_REWRITE_TAC[]) THEN
1278   REWRITE_TAC[dist; VECTOR_SUB_LZERO; VECTOR_SUB_RZERO; NORM_NEG; NORM_MUL;
1279               REAL_ABS_NORM] THEN
1280   ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
1281   ASM_SIMP_TAC[GSYM REAL_EQ_RDIV_EQ; NORM_POS_LT] THEN
1282   ASM_SIMP_TAC[REAL_DIV_REFL; REAL_LT_IMP_NZ; NORM_EQ_0] THEN
1283   TRY(FIRST_X_ASSUM MATCH_MP_TAC) THEN
1284   REWRITE_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
1285   ASM_SIMP_TAC[REAL_MUL_LINV; NORM_EQ_0] THEN
1286   MATCH_MP_TAC lemma THEN MAP_EVERY EXISTS_TAC
1287    [`inv(norm x) % x:real^N`; `inv(norm y) % y:real^N`;
1288     `(f:real^N->real^N) (inv (norm x) % x)`;
1289     `(f:real^N->real^N) (inv (norm y) % y)`] THEN
1290   REWRITE_TAC[NORM_MUL; VECTOR_MUL_ASSOC; REAL_ABS_INV; REAL_ABS_NORM] THEN
1291   ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_RINV; NORM_EQ_0] THEN
1292   ASM_REWRITE_TAC[GSYM dist; VECTOR_MUL_LID] THEN
1293   REPEAT CONJ_TAC THEN FIRST_X_ASSUM MATCH_MP_TAC THEN
1294   REWRITE_TAC[NORM_MUL; VECTOR_MUL_ASSOC; REAL_ABS_INV; REAL_ABS_NORM] THEN
1295   ASM_SIMP_TAC[REAL_MUL_LINV; REAL_MUL_RINV; NORM_EQ_0]);;
1296
1297 let ORTHOGONAL_TRANSFORMATION_LINEAR = prove
1298  (`!f:real^N->real^N. orthogonal_transformation f ==> linear f`,
1299   SIMP_TAC[orthogonal_transformation]);;
1300
1301 let ORTHOGONAL_TRANSFORMATION_INJECTIVE = prove
1302  (`!f:real^N->real^N.
1303         orthogonal_transformation f ==> !x y. f x = f y ==> x = y`,
1304   SIMP_TAC[LINEAR_INJECTIVE_0; ORTHOGONAL_TRANSFORMATION; GSYM NORM_EQ_0]);;
1305
1306 let ORTHOGONAL_TRANSFORMATION_SURJECTIVE = prove
1307  (`!f:real^N->real^N.
1308         orthogonal_transformation f ==> !y. ?x. f x = y`,
1309   MESON_TAC[LINEAR_INJECTIVE_IMP_SURJECTIVE;
1310             ORTHOGONAL_TRANSFORMATION_INJECTIVE; orthogonal_transformation]);;
1311
1312 let ORTHOGONAL_TRANSFORMATION_INVERSE_o = prove
1313  (`!f:real^N->real^N.
1314         orthogonal_transformation f
1315         ==> ?g. orthogonal_transformation g /\ g o f = I /\ f o g = I`,
1316   REPEAT STRIP_TAC THEN
1317   FIRST_ASSUM(ASSUME_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
1318   FIRST_ASSUM(ASSUME_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_INJECTIVE) THEN
1319   MP_TAC(ISPEC `f:real^N->real^N` LINEAR_INJECTIVE_LEFT_INVERSE) THEN
1320   ASM_REWRITE_TAC[] THEN MATCH_MP_TAC MONO_EXISTS THEN
1321   X_GEN_TAC `g:real^N->real^N` THEN  STRIP_TAC THEN
1322   MP_TAC(ISPECL [`f:real^N->real^N`; `g:real^N->real^N`]
1323     LINEAR_INVERSE_LEFT) THEN
1324   ASM_REWRITE_TAC[] THEN DISCH_TAC THEN ASM_REWRITE_TAC[] THEN
1325   ASM_REWRITE_TAC[ORTHOGONAL_TRANSFORMATION] THEN X_GEN_TAC `v:real^N` THEN
1326   MATCH_MP_TAC EQ_TRANS THEN
1327   EXISTS_TAC `norm((f:real^N->real^N)((g:real^N->real^N) v))` THEN
1328   CONJ_TAC THENL [ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION]; ALL_TAC] THEN
1329   RULE_ASSUM_TAC(REWRITE_RULE[FUN_EQ_THM; o_THM; I_THM]) THEN
1330   ASM_REWRITE_TAC[]);;
1331
1332 let ORTHOGONAL_TRANSFORMATION_INVERSE = prove
1333  (`!f:real^N->real^N.
1334         orthogonal_transformation f
1335         ==> ?g. orthogonal_transformation g /\
1336                 (!x. g(f x) = x) /\ (!y. f(g y) = y)`,
1337   GEN_TAC THEN
1338   DISCH_THEN(MP_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_INVERSE_o) THEN
1339   REWRITE_TAC[FUN_EQ_THM; o_THM; I_THM]);;
1340
1341 let ORTHOGONAL_TRANSFORMATION_ID = prove
1342  (`orthogonal_transformation(\x. x)`,
1343   REWRITE_TAC[orthogonal_transformation; LINEAR_ID]);;
1344
1345 let ORTHOGONAL_TRANSFORMATION_I = prove
1346  (`orthogonal_transformation I`,
1347   REWRITE_TAC[I_DEF; ORTHOGONAL_TRANSFORMATION_ID]);;
1348
1349 (* ------------------------------------------------------------------------- *)
1350 (* We can find an orthogonal matrix taking any unit vector to any other.     *)
1351 (* ------------------------------------------------------------------------- *)
1352
1353 let FINITE_INDEX_NUMSEG_SPECIAL = prove
1354  (`!s a:A.
1355         FINITE s /\ a IN s
1356         ==> ?f. (!i j. i IN 1..CARD s /\ j IN 1..CARD s /\ f i = f j
1357                        ==> i = j) /\
1358                 s = IMAGE f (1..CARD s) /\
1359                 f 1 = a`,
1360   REPEAT STRIP_TAC THEN
1361   FIRST_ASSUM(MP_TAC o GEN_REWRITE_RULE I [FINITE_INDEX_NUMSEG]) THEN
1362   DISCH_THEN(X_CHOOSE_THEN `f:num->A` STRIP_ASSUME_TAC) THEN
1363   SUBGOAL_THEN `?k. k IN 1..CARD(s:A->bool) /\ (a:A) = f k`
1364   STRIP_ASSUME_TAC THENL[ASM SET_TAC[]; ALL_TAC] THEN
1365   EXISTS_TAC `(f:num->A) o swap(1,k)` THEN
1366   SUBGOAL_THEN `1 IN 1..CARD(s:A->bool)` ASSUME_TAC THENL
1367    [REWRITE_TAC[IN_NUMSEG; LE_REFL; ARITH_RULE `1 <= x <=> ~(x = 0)`] THEN
1368     ASM_SIMP_TAC[CARD_EQ_0; ARITH_EQ] THEN ASM SET_TAC[];
1369     ALL_TAC] THEN
1370   ASM_REWRITE_TAC[o_THM; swap] THEN
1371   CONJ_TAC THENL [ASM SET_TAC[]; ALL_TAC] THEN
1372   UNDISCH_THEN `s = IMAGE (f:num->A) (1..CARD(s:A->bool))`
1373    (fun th -> GEN_REWRITE_TAC LAND_CONV [th]) THEN
1374   REWRITE_TAC[EXTENSION; IN_IMAGE; o_THM] THEN
1375   X_GEN_TAC `b:A` THEN EQ_TAC THEN
1376   DISCH_THEN(X_CHOOSE_THEN `i:num` STRIP_ASSUME_TAC) THEN
1377   EXISTS_TAC `swap(1,k) i` THEN
1378   REWRITE_TAC[swap] THEN ASM_MESON_TAC[swap]);;
1379
1380 let ORTHOGONAL_MATRIX_EXISTS_BASIS = prove
1381  (`!a:real^N.
1382         norm(a) = &1
1383         ==> ?A. orthogonal_matrix A /\ A**(basis 1) = a`,
1384   REPEAT STRIP_TAC THEN
1385   FIRST_ASSUM(MP_TAC o MATCH_MP VECTOR_IN_ORTHONORMAL_BASIS) THEN
1386   REWRITE_TAC[HAS_SIZE] THEN
1387   DISCH_THEN(X_CHOOSE_THEN `s:real^N->bool` STRIP_ASSUME_TAC) THEN
1388   MP_TAC(ISPECL [`s:real^N->bool`; `a:real^N`]
1389    FINITE_INDEX_NUMSEG_SPECIAL) THEN ASM_REWRITE_TAC[IN_NUMSEG] THEN
1390   REWRITE_TAC[TAUT `a /\ b ==> c <=> c \/ ~a \/ ~b`] THEN
1391   DISCH_THEN(X_CHOOSE_THEN `f:num->real^N`
1392    (CONJUNCTS_THEN2 ASSUME_TAC (CONJUNCTS_THEN2 (ASSUME_TAC o SYM)
1393      ASSUME_TAC))) THEN
1394   EXISTS_TAC `(lambda i j. ((f j):real^N)$i):real^N^N` THEN
1395   SIMP_TAC[CART_EQ; LAMBDA_BETA; matrix_vector_mul; BASIS_COMPONENT;
1396            IN_NUMSEG] THEN
1397   ONCE_REWRITE_TAC[COND_RAND] THEN SIMP_TAC[REAL_MUL_RZERO; SUM_DELTA] THEN
1398   ASM_REWRITE_TAC[IN_NUMSEG; REAL_MUL_RID; LE_REFL; DIMINDEX_GE_1] THEN
1399   REWRITE_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_COLUMNS] THEN
1400   SIMP_TAC[column; LAMBDA_BETA] THEN CONJ_TAC THENL
1401    [X_GEN_TAC `i:num` THEN REPEAT STRIP_TAC THEN
1402     MATCH_MP_TAC EQ_TRANS THEN
1403     EXISTS_TAC `norm((f:num->real^N) i)` THEN CONJ_TAC THENL
1404      [AP_TERM_TAC THEN ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA];
1405       ASM_MESON_TAC[IN_IMAGE; IN_NUMSEG]];
1406     MAP_EVERY X_GEN_TAC [`i:num`; `j:num`] THEN STRIP_TAC THEN
1407     SUBGOAL_THEN `orthogonal ((f:num->real^N) i) (f j)` MP_TAC THENL
1408      [ASM_MESON_TAC[pairwise; IN_IMAGE; IN_NUMSEG]; ALL_TAC] THEN
1409     MATCH_MP_TAC EQ_IMP THEN BINOP_TAC THEN
1410     ASM_SIMP_TAC[CART_EQ; LAMBDA_BETA]]);;
1411
1412 let ORTHOGONAL_TRANSFORMATION_EXISTS_1 = prove
1413  (`!a b:real^N.
1414         norm(a) = &1 /\ norm(b) = &1
1415         ==> ?f. orthogonal_transformation f /\ f a = b`,
1416   REPEAT STRIP_TAC THEN
1417   MP_TAC(ISPEC `b:real^N` ORTHOGONAL_MATRIX_EXISTS_BASIS) THEN
1418   MP_TAC(ISPEC `a:real^N` ORTHOGONAL_MATRIX_EXISTS_BASIS) THEN
1419   ASM_REWRITE_TAC[] THEN
1420   DISCH_THEN(X_CHOOSE_THEN `A:real^N^N` (STRIP_ASSUME_TAC o GSYM)) THEN
1421   DISCH_THEN(X_CHOOSE_THEN `B:real^N^N` (STRIP_ASSUME_TAC o GSYM)) THEN
1422   EXISTS_TAC `\x:real^N. ((B:real^N^N) ** transp(A:real^N^N)) ** x` THEN
1423   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX; MATRIX_VECTOR_MUL_LINEAR;
1424               MATRIX_OF_MATRIX_VECTOR_MUL] THEN
1425   ASM_SIMP_TAC[ORTHOGONAL_MATRIX_MUL; ORTHOGONAL_MATRIX_TRANSP] THEN
1426   REWRITE_TAC[GSYM MATRIX_VECTOR_MUL_ASSOC] THEN AP_TERM_TAC THEN
1427   RULE_ASSUM_TAC(REWRITE_RULE[ORTHOGONAL_MATRIX]) THEN
1428   ASM_REWRITE_TAC[MATRIX_VECTOR_MUL_ASSOC; MATRIX_VECTOR_MUL_LID]);;
1429
1430 let ORTHOGONAL_TRANSFORMATION_EXISTS = prove
1431  (`!a b:real^N.
1432         norm(a) = norm(b) ==> ?f. orthogonal_transformation f /\ f a = b`,
1433   REPEAT GEN_TAC THEN ASM_CASES_TAC `b:real^N = vec 0` THEN
1434   ASM_SIMP_TAC[NORM_0; NORM_EQ_0] THENL
1435    [MESON_TAC[ORTHOGONAL_TRANSFORMATION_ID]; ALL_TAC] THEN
1436   ASM_CASES_TAC `a:real^N = vec 0` THENL
1437    [ASM_MESON_TAC[NORM_0; NORM_EQ_0]; ALL_TAC] THEN
1438   DISCH_TAC THEN
1439   MP_TAC(ISPECL [`inv(norm a) % a:real^N`; `inv(norm b) % b:real^N`]
1440                 ORTHOGONAL_TRANSFORMATION_EXISTS_1) THEN
1441   REWRITE_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
1442   ASM_SIMP_TAC[NORM_EQ_0; REAL_MUL_LINV] THEN
1443   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
1444   DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
1445   ASM_REWRITE_TAC[] THEN
1446   FIRST_ASSUM(ASSUME_TAC o MATCH_MP LINEAR_CMUL o
1447               MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
1448   ASM_REWRITE_TAC[VECTOR_ARITH
1449     `a % x:real^N = a % y <=> a % (x - y) = vec 0`] THEN
1450   ASM_REWRITE_TAC[VECTOR_MUL_EQ_0; REAL_INV_EQ_0; NORM_EQ_0; VECTOR_SUB_EQ]);;
1451
1452 (* ------------------------------------------------------------------------- *)
1453 (* Rotation, reflection, rotoinversion.                                      *)
1454 (* ------------------------------------------------------------------------- *)
1455
1456 let rotation_matrix = new_definition
1457  `rotation_matrix Q <=> orthogonal_matrix Q /\ det(Q) = &1`;;
1458
1459 let rotoinversion_matrix = new_definition
1460  `rotoinversion_matrix Q <=> orthogonal_matrix Q /\ det(Q) = -- &1`;;
1461
1462 let ORTHOGONAL_ROTATION_OR_ROTOINVERSION = prove
1463  (`!Q. orthogonal_matrix Q <=> rotation_matrix Q \/ rotoinversion_matrix Q`,
1464   MESON_TAC[rotation_matrix; rotoinversion_matrix; DET_ORTHOGONAL_MATRIX]);;
1465
1466 let ROTATION_MATRIX_2 = prove
1467  (`!A:real^2^2. rotation_matrix A <=>
1468                 A$1$1 pow 2 + A$2$1 pow 2 = &1 /\
1469                 A$1$1 = A$2$2 /\ A$1$2 = --(A$2$1)`,
1470   REWRITE_TAC[rotation_matrix; ORTHOGONAL_MATRIX_2; DET_2] THEN
1471   CONV_TAC REAL_RING);;
1472
1473 (* ------------------------------------------------------------------------- *)
1474 (* Slightly stronger results giving rotation, but only in >= 2 dimensions.   *)
1475 (* ------------------------------------------------------------------------- *)
1476
1477 let ROTATION_MATRIX_EXISTS_BASIS = prove
1478  (`!a:real^N.
1479         2 <= dimindex(:N) /\ norm(a) = &1
1480         ==> ?A. rotation_matrix A /\ A**(basis 1) = a`,
1481   REPEAT STRIP_TAC THEN
1482   FIRST_X_ASSUM(X_CHOOSE_THEN `A:real^N^N` STRIP_ASSUME_TAC o
1483    MATCH_MP ORTHOGONAL_MATRIX_EXISTS_BASIS) THEN
1484   FIRST_ASSUM(DISJ_CASES_TAC o GEN_REWRITE_RULE I
1485    [ORTHOGONAL_ROTATION_OR_ROTOINVERSION])
1486   THENL [ASM_MESON_TAC[]; ALL_TAC] THEN
1487   EXISTS_TAC `transp(lambda i. if i = dimindex(:N) then -- &1 % transp A$i
1488                                else (transp A:real^N^N)$i):real^N^N` THEN
1489   REWRITE_TAC[rotation_matrix; DET_TRANSP] THEN REPEAT CONJ_TAC THENL
1490    [REWRITE_TAC[ORTHOGONAL_MATRIX_TRANSP];
1491     SIMP_TAC[DET_ROW_MUL; DIMINDEX_GE_1; LE_REFL] THEN
1492     MATCH_MP_TAC(REAL_ARITH `x = -- &1 ==> -- &1 * x = &1`) THEN
1493     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [rotoinversion_matrix]) THEN
1494     DISCH_THEN(SUBST1_TAC o SYM o CONJUNCT2) THEN
1495     GEN_REWRITE_TAC RAND_CONV [GSYM DET_TRANSP] THEN
1496     AP_TERM_TAC THEN SIMP_TAC[CART_EQ; LAMBDA_BETA] THEN MESON_TAC[];
1497     FIRST_X_ASSUM(SUBST1_TAC o SYM) THEN
1498     SIMP_TAC[matrix_vector_mul; LAMBDA_BETA; CART_EQ; transp;
1499              BASIS_COMPONENT] THEN
1500     ONCE_REWRITE_TAC[REAL_ARITH
1501       `x * (if p then &1 else &0) = if p then x else &0`] THEN
1502     ASM_SIMP_TAC[ARITH_RULE `2 <= n ==> ~(1 = n)`; LAMBDA_BETA]] THEN
1503   FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I
1504    [GSYM ORTHOGONAL_MATRIX_TRANSP]) THEN
1505   SPEC_TAC(`transp(A:real^N^N)`,`B:real^N^N`) THEN GEN_TAC THEN
1506   SUBGOAL_THEN `!i. 1 <= i /\ i <= dimindex(:N)
1507                     ==> row i ((lambda i. if i = dimindex(:N) then -- &1 % B$i
1508                                 else (B:real^N^N)$i):real^N^N) =
1509                         if i = dimindex(:N) then --(row i B) else row i B`
1510   ASSUME_TAC THENL
1511    [SIMP_TAC[row; LAMBDA_BETA; LAMBDA_ETA; VECTOR_MUL_LID; VECTOR_MUL_LNEG];
1512     ASM_SIMP_TAC[ORTHOGONAL_MATRIX_ORTHONORMAL_ROWS] THEN
1513     ASM_MESON_TAC[ORTHOGONAL_LNEG; ORTHOGONAL_RNEG; NORM_NEG]]);;
1514
1515 let ROTATION_EXISTS_1 = prove
1516  (`!a b:real^N.
1517         2 <= dimindex(:N) /\ norm(a) = &1 /\ norm(b) = &1
1518         ==> ?f. orthogonal_transformation f /\ det(matrix f) = &1 /\ f a = b`,
1519   REPEAT STRIP_TAC THEN
1520   MP_TAC(ISPEC `b:real^N` ROTATION_MATRIX_EXISTS_BASIS) THEN
1521   MP_TAC(ISPEC `a:real^N` ROTATION_MATRIX_EXISTS_BASIS) THEN
1522   ASM_REWRITE_TAC[rotation_matrix] THEN
1523   DISCH_THEN(X_CHOOSE_THEN `A:real^N^N`
1524    (CONJUNCTS_THEN2 STRIP_ASSUME_TAC (ASSUME_TAC o SYM))) THEN
1525   DISCH_THEN(X_CHOOSE_THEN `B:real^N^N`
1526    (CONJUNCTS_THEN2 STRIP_ASSUME_TAC (ASSUME_TAC o SYM))) THEN
1527   EXISTS_TAC `\x:real^N. ((B:real^N^N) ** transp(A:real^N^N)) ** x` THEN
1528   REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_MATRIX; MATRIX_VECTOR_MUL_LINEAR;
1529               MATRIX_OF_MATRIX_VECTOR_MUL; DET_MUL; DET_TRANSP] THEN
1530   ASM_SIMP_TAC[ORTHOGONAL_MATRIX_MUL; ORTHOGONAL_MATRIX_TRANSP] THEN
1531   REWRITE_TAC[GSYM MATRIX_VECTOR_MUL_ASSOC; REAL_MUL_LID] THEN AP_TERM_TAC THEN
1532   RULE_ASSUM_TAC(REWRITE_RULE[ORTHOGONAL_MATRIX]) THEN
1533   ASM_REWRITE_TAC[MATRIX_VECTOR_MUL_ASSOC; MATRIX_VECTOR_MUL_LID]);;
1534
1535 let ROTATION_EXISTS = prove
1536  (`!a b:real^N.
1537         2 <= dimindex(:N) /\ norm(a) = norm(b)
1538         ==> ?f. orthogonal_transformation f /\ det(matrix f) = &1 /\ f a = b`,
1539   REPEAT GEN_TAC THEN ASM_CASES_TAC `b:real^N = vec 0` THEN
1540   ASM_SIMP_TAC[NORM_0; NORM_EQ_0] THENL
1541    [MESON_TAC[ORTHOGONAL_TRANSFORMATION_ID; MATRIX_ID; DET_I]; ALL_TAC] THEN
1542   ASM_CASES_TAC `a:real^N = vec 0` THENL
1543    [ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_ID; MATRIX_ID; DET_I; NORM_0;
1544                   NORM_EQ_0]; ALL_TAC] THEN
1545   DISCH_TAC THEN
1546   MP_TAC(ISPECL [`inv(norm a) % a:real^N`; `inv(norm b) % b:real^N`]
1547                 ROTATION_EXISTS_1) THEN
1548   REWRITE_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM] THEN
1549   ASM_SIMP_TAC[NORM_EQ_0; REAL_MUL_LINV] THEN
1550   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
1551   DISCH_THEN(CONJUNCTS_THEN2 ASSUME_TAC MP_TAC) THEN
1552   ASM_REWRITE_TAC[] THEN
1553   FIRST_ASSUM(ASSUME_TAC o MATCH_MP LINEAR_CMUL o
1554               MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
1555   ASM_REWRITE_TAC[VECTOR_ARITH
1556    `a % x:real^N = a % y <=> a % (x - y) = vec 0`] THEN
1557   ASM_REWRITE_TAC[VECTOR_MUL_EQ_0; REAL_INV_EQ_0; NORM_EQ_0; VECTOR_SUB_EQ]);;
1558
1559 let ROTATION_RIGHTWARD_LINE = prove
1560  (`!a:real^N k.
1561         1 <= k /\ k <= dimindex(:N)
1562         ==> ?b f. orthogonal_transformation f /\
1563                   (2 <= dimindex(:N) ==> det(matrix f) = &1) /\
1564                   f(b % basis k) = a /\
1565                   &0 <= b`,
1566   REPEAT STRIP_TAC THEN EXISTS_TAC `norm(a:real^N)` THEN
1567   ASM_SIMP_TAC[VECTOR_MUL_COMPONENT; BASIS_COMPONENT; LE_REFL; DIMINDEX_GE_1;
1568                REAL_MUL_RID; NORM_POS_LE; LT_IMP_LE; LTE_ANTISYM] THEN
1569   REWRITE_TAC[ARITH_RULE `2 <= n <=> 1 <= n /\ ~(n = 1)`; DIMINDEX_GE_1] THEN
1570   ASM_CASES_TAC `dimindex(:N) = 1` THEN ASM_REWRITE_TAC[] THENL
1571    [MATCH_MP_TAC ORTHOGONAL_TRANSFORMATION_EXISTS;
1572     MATCH_MP_TAC ROTATION_EXISTS] THEN
1573    ASM_SIMP_TAC[NORM_MUL; NORM_BASIS; LE_REFL; DIMINDEX_GE_1] THEN
1574    REWRITE_TAC[REAL_ABS_NORM; REAL_MUL_RID] THEN
1575    MATCH_MP_TAC(ARITH_RULE `~(n = 1) /\ 1 <= n ==> 2 <= n`) THEN
1576    ASM_REWRITE_TAC[DIMINDEX_GE_1]);;
1577
1578 (* ------------------------------------------------------------------------- *)
1579 (* We can always rotate so that a hyperplane is "horizontal".                *)
1580 (* ------------------------------------------------------------------------- *)
1581
1582 let ROTATION_LOWDIM_HORIZONTAL = prove
1583  (`!s:real^N->bool.
1584         dim s < dimindex(:N)
1585         ==> ?f. orthogonal_transformation f /\ det(matrix f) = &1 /\
1586                (IMAGE f s) SUBSET {z | z$(dimindex(:N)) = &0}`,
1587   GEN_TAC THEN ASM_CASES_TAC `dim(s:real^N->bool) = 0` THENL
1588    [RULE_ASSUM_TAC(REWRITE_RULE[DIM_EQ_0]) THEN DISCH_TAC THEN
1589     EXISTS_TAC `\x:real^N. x` THEN
1590     REWRITE_TAC[ORTHOGONAL_TRANSFORMATION_ID; MATRIX_ID; DET_I] THEN
1591     FIRST_X_ASSUM(MATCH_MP_TAC o MATCH_MP (SET_RULE
1592      `s SUBSET {a} ==> a IN t ==> IMAGE (\x. x) s SUBSET t`)) THEN
1593     SIMP_TAC[IN_ELIM_THM; VEC_COMPONENT; LE_REFL; DIMINDEX_GE_1];
1594     DISCH_TAC] THEN
1595   SUBGOAL_THEN `2 <= dimindex(:N)` ASSUME_TAC THENL
1596    [ASM_ARITH_TAC; ALL_TAC] THEN
1597   FIRST_X_ASSUM(X_CHOOSE_THEN `a:real^N` STRIP_ASSUME_TAC o MATCH_MP
1598     LOWDIM_SUBSET_HYPERPLANE) THEN
1599   MP_TAC(ISPECL [`a:real^N`; `norm(a:real^N) % basis(dimindex(:N)):real^N`]
1600         ROTATION_EXISTS) THEN
1601   ASM_SIMP_TAC[NORM_MUL; NORM_BASIS; LE_REFL; DIMINDEX_GE_1] THEN
1602   REWRITE_TAC[REAL_ABS_NORM; REAL_MUL_RID] THEN
1603   MATCH_MP_TAC MONO_EXISTS THEN X_GEN_TAC `f:real^N->real^N` THEN
1604   STRIP_TAC THEN ASM_REWRITE_TAC[FORALL_IN_IMAGE; SUBSET; IN_ELIM_THM] THEN
1605   X_GEN_TAC `x:real^N` THEN DISCH_TAC THEN
1606   SUBGOAL_THEN `(f:real^N->real^N) x dot (f a) = &0` MP_TAC THENL
1607    [RULE_ASSUM_TAC(REWRITE_RULE[orthogonal_transformation]) THEN
1608     ASM_REWRITE_TAC[] THEN ONCE_REWRITE_TAC[DOT_SYM] THEN
1609     FIRST_X_ASSUM(MP_TAC o GEN_REWRITE_RULE I [SUBSET]) THEN
1610     ASM_SIMP_TAC[SPAN_SUPERSET; IN_ELIM_THM];
1611     ASM_SIMP_TAC[DOT_BASIS; LE_REFL; DIMINDEX_GE_1; DOT_RMUL] THEN
1612     ASM_REWRITE_TAC[REAL_ENTIRE; NORM_EQ_0]]);;
1613
1614 let ORTHOGONAL_TRANSFORMATION_LOWDIM_HORIZONTAL = prove
1615  (`!s:real^N->bool.
1616         dim s < dimindex(:N)
1617         ==> ?f. orthogonal_transformation f /\
1618                (IMAGE f s) SUBSET {z | z$(dimindex(:N)) = &0}`,
1619   GEN_TAC THEN DISCH_THEN(MP_TAC o MATCH_MP ROTATION_LOWDIM_HORIZONTAL) THEN
1620   MATCH_MP_TAC MONO_EXISTS THEN SIMP_TAC[]);;
1621
1622 (* ------------------------------------------------------------------------- *)
1623 (* Extract scaling, translation and linear invariance theorems.              *)
1624 (* For the linear case, chain through some basic consequences automatically, *)
1625 (* e.g. norm-preserving and linear implies injective.                        *)
1626 (* ------------------------------------------------------------------------- *)
1627
1628 let SCALING_THEOREMS v =
1629   let th1 = UNDISCH(snd(EQ_IMP_RULE(ISPEC v NORM_POS_LT))) in
1630   let t = rand(concl th1) in
1631   end_itlist CONJ (map (C MP th1 o SPEC t) (!scaling_theorems));;
1632
1633 let TRANSLATION_INVARIANTS x =
1634   end_itlist CONJ (mapfilter (ISPEC x) (!invariant_under_translation));;
1635
1636 let USABLE_CONCLUSION f ths th =
1637   let ith = PURE_REWRITE_RULE[RIGHT_FORALL_IMP_THM] (ISPEC f th) in
1638   let bod = concl ith in
1639   let cjs = conjuncts(fst(dest_imp bod)) in
1640   let ths = map (fun t -> find(fun th -> aconv (concl th) t) ths) cjs in
1641   GEN_ALL(MP ith (end_itlist CONJ ths));;
1642
1643 let LINEAR_INVARIANTS =
1644   let sths = (CONJUNCTS o prove)
1645    (`(!f:real^M->real^N.
1646          linear f /\ (!x. norm(f x) = norm x)
1647          ==> (!x y. f x = f y ==> x = y)) /\
1648      (!f:real^N->real^N.
1649          linear f /\ (!x. norm(f x) = norm x) ==> (!y. ?x. f x = y)) /\
1650      (!f:real^N->real^N. linear f /\ (!x y. f x = f y ==> x = y)
1651                          ==> (!y. ?x. f x = y)) /\
1652      (!f:real^N->real^N. linear f /\ (!y. ?x. f x = y)
1653                          ==> (!x y. f x = f y ==> x = y))`,
1654     CONJ_TAC THENL
1655      [ONCE_REWRITE_TAC[GSYM VECTOR_SUB_EQ] THEN
1656       SIMP_TAC[GSYM LINEAR_SUB; GSYM NORM_EQ_0];
1657       MESON_TAC[ORTHOGONAL_TRANSFORMATION_SURJECTIVE;
1658                 ORTHOGONAL_TRANSFORMATION_INJECTIVE; ORTHOGONAL_TRANSFORMATION;
1659                 LINEAR_SURJECTIVE_IFF_INJECTIVE]]) in
1660   fun f ths ->
1661     let ths' = ths @ mapfilter (USABLE_CONCLUSION f ths) sths in
1662     end_itlist CONJ
1663      (mapfilter (USABLE_CONCLUSION f ths') (!invariant_under_linear));;
1664
1665 (* ------------------------------------------------------------------------- *)
1666 (* Tactic to pick WLOG a particular point as the origin. The conversion form *)
1667 (* assumes it's the outermost universal variable; the tactic is more general *)
1668 (* and allows any free or outer universally quantified variable. The list    *)
1669 (* "avoid" is the points not to translate. There is also a tactic to help in *)
1670 (* proving new translation theorems, which uses similar machinery.           *)
1671 (* ------------------------------------------------------------------------- *)
1672
1673 let GEOM_ORIGIN_CONV,GEOM_TRANSLATE_CONV =
1674   let pth = prove
1675    (`!a:real^N. a = a + vec 0 /\
1676                 {} = IMAGE (\x. a + x) {} /\
1677                 {} = IMAGE (IMAGE (\x. a + x)) {} /\
1678                 (:real^N) = IMAGE (\x. a + x) (:real^N) /\
1679                 (:real^N->bool) = IMAGE (IMAGE (\x. a + x)) (:real^N->bool) /\
1680                 [] = MAP (\x. a + x) []`,
1681     REWRITE_TAC[IMAGE_CLAUSES; VECTOR_ADD_RID; MAP] THEN
1682     REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
1683     REWRITE_TAC[SURJECTIVE_IMAGE] THEN
1684     REWRITE_TAC[VECTOR_ARITH `a + y:real^N = x <=> y = x - a`; EXISTS_REFL])
1685   and qth = prove
1686    (`!a:real^N.
1687         ((!P. (!x. P x) <=> (!x. P (a + x))) /\
1688          (!P. (?x. P x) <=> (?x. P (a + x))) /\
1689          (!Q. (!s. Q s) <=> (!s. Q(IMAGE (\x. a + x) s))) /\
1690          (!Q. (?s. Q s) <=> (?s. Q(IMAGE (\x. a + x) s))) /\
1691          (!Q. (!s. Q s) <=> (!s. Q(IMAGE (IMAGE (\x. a + x)) s))) /\
1692          (!Q. (?s. Q s) <=> (?s. Q(IMAGE (IMAGE (\x. a + x)) s))) /\
1693          (!P. (!g:real^1->real^N. P g) <=> (!g. P ((\x. a + x) o g))) /\
1694          (!P. (?g:real^1->real^N. P g) <=> (?g. P ((\x. a + x) o g))) /\
1695          (!P. (!g:num->real^N. P g) <=> (!g. P ((\x. a + x) o g))) /\
1696          (!P. (?g:num->real^N. P g) <=> (?g. P ((\x. a + x) o g))) /\
1697          (!Q. (!l. Q l) <=> (!l. Q(MAP (\x. a + x) l))) /\
1698          (!Q. (?l. Q l) <=> (?l. Q(MAP (\x. a + x) l)))) /\
1699         ((!P. {x | P x} = IMAGE (\x. a + x) {x | P(a + x)}) /\
1700          (!Q. {s | Q s} =
1701               IMAGE (IMAGE (\x. a + x)) {s | Q(IMAGE (\x. a + x) s)}) /\
1702          (!R. {l | R l} = IMAGE (MAP (\x. a + x)) {l | R(MAP (\x. a + x) l)}))`,
1703     GEN_TAC THEN MATCH_MP_TAC QUANTIFY_SURJECTION_HIGHER_THM THEN
1704     X_GEN_TAC `y:real^N` THEN EXISTS_TAC `y - a:real^N` THEN
1705     VECTOR_ARITH_TAC) in
1706   let GEOM_ORIGIN_CONV avoid tm =
1707     let x,tm0 = dest_forall tm in
1708     let th0 = ISPEC x pth in
1709     let x' = genvar(type_of x) in
1710     let ith = ISPEC x' qth in
1711     let th1 = PARTIAL_EXPAND_QUANTS_CONV avoid (ASSUME(concl ith)) tm0 in
1712     let th2 = CONV_RULE(RAND_CONV(SUBS_CONV(CONJUNCTS th0))) th1 in
1713     let th3 = INST[x,x'] (PROVE_HYP ith th2) in
1714     let ths = TRANSLATION_INVARIANTS x in
1715     let thr = REFL x in
1716     let th4 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV)
1717       [BETA_THM;ADD_ASSUM(concl thr) ths] th3 in
1718     let th5 = MK_FORALL x (PROVE_HYP thr th4) in
1719     GEN_REWRITE_RULE (RAND_CONV o TRY_CONV) [FORALL_SIMP] th5
1720   and GEOM_TRANSLATE_CONV avoid a tm =
1721     let cth = CONJUNCT2(ISPEC a pth)
1722     and vth = ISPEC a qth in
1723     let th1 = PARTIAL_EXPAND_QUANTS_CONV avoid (ASSUME(concl vth)) tm in
1724     let th2 = CONV_RULE(RAND_CONV(SUBS_CONV(CONJUNCTS cth))) th1 in
1725     let th3 = PROVE_HYP vth th2 in
1726     let ths = TRANSLATION_INVARIANTS a in
1727     let thr = REFL a in
1728     let th4 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV)
1729         [BETA_THM;ADD_ASSUM(concl thr) ths] th3 in
1730     PROVE_HYP thr th4 in
1731   GEOM_ORIGIN_CONV,GEOM_TRANSLATE_CONV;;
1732
1733 let GEN_GEOM_ORIGIN_TAC x avoid (asl,w as gl) =
1734   let avs,bod = strip_forall w
1735   and avs' = subtract (frees w) (freesl(map (concl o snd) asl)) in
1736   (MAP_EVERY X_GEN_TAC avs THEN
1737    MAP_EVERY (fun t -> SPEC_TAC(t,t)) (rev(subtract (avs@avs') [x])) THEN
1738    SPEC_TAC(x,x) THEN CONV_TAC(GEOM_ORIGIN_CONV avoid)) gl;;
1739
1740 let GEOM_ORIGIN_TAC x = GEN_GEOM_ORIGIN_TAC x [];;
1741
1742 let GEOM_TRANSLATE_TAC avoid (asl,w) =
1743   let a,bod = dest_forall w in
1744   let n = length(fst(strip_forall bod)) in
1745   (X_GEN_TAC a THEN
1746    CONV_TAC(funpow n BINDER_CONV (LAND_CONV(GEOM_TRANSLATE_CONV avoid a))) THEN
1747    REWRITE_TAC[]) (asl,w);;
1748
1749 (* ------------------------------------------------------------------------- *)
1750 (* Rename existential variables in conclusion to fresh genvars.              *)
1751 (* ------------------------------------------------------------------------- *)
1752
1753 let EXISTS_GENVAR_RULE =
1754   let rec rule vs th =
1755     match vs with
1756       [] -> th
1757     | v::ovs -> let x,bod = dest_exists(concl th) in
1758                 let th1 = rule ovs (ASSUME bod) in
1759                 let th2 = SIMPLE_CHOOSE x (SIMPLE_EXISTS x th1) in
1760                 PROVE_HYP th (CONV_RULE (GEN_ALPHA_CONV v) th2) in
1761   fun th -> rule (map (genvar o type_of) (fst(strip_exists(concl th)))) th;;
1762
1763 (* ------------------------------------------------------------------------- *)
1764 (* Rotate so that WLOG some point is a +ve multiple of basis vector k.       *)
1765 (* For general N, it's better to use k = 1 so the side-condition can be      *)
1766 (* discharged. For dimensions 1, 2 and 3 anything will work automatically.   *)
1767 (* Could generalize by asking the user to prove theorem 1 <= k <= N.         *)
1768 (* ------------------------------------------------------------------------- *)
1769
1770 let GEOM_BASIS_MULTIPLE_RULE =
1771   let pth = prove
1772    (`!f. orthogonal_transformation (f:real^N->real^N)
1773          ==> (vec 0 = f(vec 0) /\
1774               {} = IMAGE f {} /\
1775               {} = IMAGE (IMAGE f) {} /\
1776               (:real^N) = IMAGE f (:real^N) /\
1777               (:real^N->bool) = IMAGE (IMAGE f) (:real^N->bool) /\
1778               [] = MAP f []) /\
1779              ((!P. (!x. P x) <=> (!x. P (f x))) /\
1780               (!P. (?x. P x) <=> (?x. P (f x))) /\
1781               (!Q. (!s. Q s) <=> (!s. Q (IMAGE f s))) /\
1782               (!Q. (?s. Q s) <=> (?s. Q (IMAGE f s))) /\
1783               (!Q. (!s. Q s) <=> (!s. Q (IMAGE (IMAGE f) s))) /\
1784               (!Q. (?s. Q s) <=> (?s. Q (IMAGE (IMAGE f) s))) /\
1785               (!P. (!g:real^1->real^N. P g) <=> (!g. P (f o g))) /\
1786               (!P. (?g:real^1->real^N. P g) <=> (?g. P (f o g))) /\
1787               (!P. (!g:num->real^N. P g) <=> (!g. P (f o g))) /\
1788               (!P. (?g:num->real^N. P g) <=> (?g. P (f o g))) /\
1789               (!Q. (!l. Q l) <=> (!l. Q(MAP f l))) /\
1790               (!Q. (?l. Q l) <=> (?l. Q(MAP f l)))) /\
1791              ((!P. {x | P x} = IMAGE f {x | P(f x)}) /\
1792               (!Q. {s | Q s} = IMAGE (IMAGE f) {s | Q(IMAGE f s)}) /\
1793               (!R. {l | R l} = IMAGE (MAP f) {l | R(MAP f l)}))`,
1794     REPEAT GEN_TAC THEN DISCH_TAC THEN
1795     FIRST_ASSUM(ASSUME_TAC o
1796           MATCH_MP ORTHOGONAL_TRANSFORMATION_SURJECTIVE) THEN
1797     CONJ_TAC THENL
1798      [REWRITE_TAC[IMAGE_CLAUSES; MAP] THEN
1799       FIRST_ASSUM(ASSUME_TAC o MATCH_MP ORTHOGONAL_TRANSFORMATION_LINEAR) THEN
1800       CONJ_TAC THENL [ASM_MESON_TAC[LINEAR_0]; ALL_TAC] THEN
1801       REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
1802       ASM_REWRITE_TAC[SURJECTIVE_IMAGE];
1803       MATCH_MP_TAC QUANTIFY_SURJECTION_HIGHER_THM THEN ASM_REWRITE_TAC[]])
1804   and oth = prove
1805    (`!f:real^N->real^N.
1806         orthogonal_transformation f /\
1807         (2 <= dimindex(:N) ==> det(matrix f) = &1)
1808         ==> linear f /\
1809             (!x y. f x = f y ==> x = y) /\
1810             (!y. ?x. f x = y) /\
1811             (!x. norm(f x) = norm x) /\
1812             (2 <= dimindex(:N) ==> det(matrix f) = &1)`,
1813     GEN_TAC THEN STRIP_TAC THEN ASM_REWRITE_TAC[] THEN REPEAT CONJ_TAC THENL
1814      [ASM_SIMP_TAC[ORTHOGONAL_TRANSFORMATION_LINEAR];
1815       ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_INJECTIVE];
1816       ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION_SURJECTIVE];
1817       ASM_MESON_TAC[ORTHOGONAL_TRANSFORMATION]])
1818   and arithconv = REWRITE_CONV[DIMINDEX_1; DIMINDEX_2; DIMINDEX_3;
1819                                ARITH_RULE `1 <= 1`; DIMINDEX_GE_1] THENC
1820                   NUM_REDUCE_CONV in
1821   fun k tm ->
1822     let x,bod = dest_forall tm in
1823     let th0 = ISPECL [x; mk_small_numeral k] ROTATION_RIGHTWARD_LINE in
1824     let th1 = EXISTS_GENVAR_RULE
1825      (MP th0 (EQT_ELIM(arithconv(lhand(concl th0))))) in
1826     let [a;f],tm1 = strip_exists(concl th1) in
1827     let th_orth,th2 = CONJ_PAIR(ASSUME tm1) in
1828     let th_det,th2a = CONJ_PAIR th2 in
1829     let th_works,th_zero = CONJ_PAIR th2a in
1830     let thc,thq = CONJ_PAIR(PROVE_HYP th2 (UNDISCH(ISPEC f pth))) in
1831     let th3 = CONV_RULE(RAND_CONV(SUBS_CONV(GSYM th_works::CONJUNCTS thc)))
1832                (EXPAND_QUANTS_CONV(ASSUME(concl thq)) bod) in
1833     let th4 = PROVE_HYP thq th3 in
1834     let thps = CONJUNCTS(MATCH_MP oth (CONJ th_orth th_det)) in
1835     let th5 = LINEAR_INVARIANTS f thps in
1836     let th6 = PROVE_HYP th_orth
1837      (GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV) [BETA_THM; th5] th4) in
1838     let ntm = mk_forall(a,mk_imp(concl th_zero,rand(concl th6))) in
1839     let th7 = MP(SPEC a (ASSUME ntm)) th_zero in
1840     let th8 = DISCH ntm (EQ_MP (SYM th6) th7) in
1841     if intersect (frees(concl th8)) [a;f] = [] then
1842       let th9 = PROVE_HYP th1 (itlist SIMPLE_CHOOSE [a;f] th8) in
1843       let th10 = DISCH ntm (GEN x (UNDISCH th9)) in
1844       let a' = variant (frees(concl th10))
1845                 (mk_var(fst(dest_var x),snd(dest_var a))) in
1846       CONV_RULE(LAND_CONV (GEN_ALPHA_CONV a')) th10
1847     else
1848       let mtm = list_mk_forall([a;f],mk_imp(hd(hyp th8),rand(concl th6))) in
1849       let th9 = EQ_MP (SYM th6) (UNDISCH(SPECL [a;f] (ASSUME mtm))) in
1850       let th10 = itlist SIMPLE_CHOOSE [a;f] (DISCH mtm th9) in
1851       let th11 = GEN x (PROVE_HYP th1 th10) in
1852       MATCH_MP MONO_FORALL th11;;
1853
1854 let GEOM_BASIS_MULTIPLE_TAC k l (asl,w as gl) =
1855   let avs,bod = strip_forall w
1856   and avs' = subtract (frees w) (freesl(map (concl o snd) asl)) in
1857   (MAP_EVERY X_GEN_TAC avs THEN
1858    MAP_EVERY (fun t -> SPEC_TAC(t,t)) (rev(subtract (avs@avs') [l])) THEN
1859    SPEC_TAC(l,l) THEN
1860    W(MATCH_MP_TAC o GEOM_BASIS_MULTIPLE_RULE k o snd)) gl;;
1861
1862 (* ------------------------------------------------------------------------- *)
1863 (* Create invariance theorems automatically, in simple cases. If there are   *)
1864 (* any nested quantifiers, this will need surjectivity. It's often possible  *)
1865 (* to prove a stronger theorem by more delicate manual reasoning, so this    *)
1866 (* isn't used nearly as often as GEOM_TRANSLATE_CONV / GEOM_TRANSLATE_TAC.   *)
1867 (* As a small step, some ad-hoc rewrites analogous to FORALL_IN_IMAGE are    *)
1868 (* tried if the first step doesn't finish the goal, but it's very ad hoc.    *)
1869 (* ------------------------------------------------------------------------- *)
1870
1871 let GEOM_TRANSFORM_TAC =
1872   let cth0 = prove
1873    (`!f:real^M->real^N.
1874           linear f
1875           ==> vec 0 = f(vec 0) /\
1876               {} = IMAGE f {} /\
1877               {} = IMAGE (IMAGE f) {}`,
1878     REWRITE_TAC[IMAGE_CLAUSES] THEN MESON_TAC[LINEAR_0])
1879   and cth1 = prove
1880    (`!f:real^M->real^N.
1881           (!y. ?x. f x = y)
1882           ==> (:real^N) = IMAGE f (:real^M) /\
1883               (:real^N->bool) = IMAGE (IMAGE f) (:real^M->bool)`,
1884     REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
1885     REWRITE_TAC[SURJECTIVE_IMAGE])
1886   and sths = (CONJUNCTS o prove)
1887    (`(!f:real^M->real^N.
1888          linear f /\ (!x. norm(f x) = norm x)
1889          ==> (!x y. f x = f y ==> x = y)) /\
1890      (!f:real^N->real^N.
1891          linear f /\ (!x. norm(f x) = norm x) ==> (!y. ?x. f x = y)) /\
1892      (!f:real^N->real^N. linear f /\ (!x y. f x = f y ==> x = y)
1893                          ==> (!y. ?x. f x = y)) /\
1894      (!f:real^N->real^N. linear f /\ (!y. ?x. f x = y)
1895                          ==> (!x y. f x = f y ==> x = y))`,
1896     CONJ_TAC THENL
1897      [ONCE_REWRITE_TAC[GSYM VECTOR_SUB_EQ] THEN
1898       SIMP_TAC[GSYM LINEAR_SUB; GSYM NORM_EQ_0];
1899       MESON_TAC[ORTHOGONAL_TRANSFORMATION_SURJECTIVE;
1900                 ORTHOGONAL_TRANSFORMATION_INJECTIVE; ORTHOGONAL_TRANSFORMATION;
1901                 LINEAR_SURJECTIVE_IFF_INJECTIVE]])
1902   and aths = (CONJUNCTS o prove)
1903    (`(!f s P. (!y. y IN IMAGE f s ==> P y) <=> (!x. x IN s ==> P(f x))) /\
1904      (!f s P. (!u. u IN IMAGE (IMAGE f) s ==> P u) <=>
1905               (!t. t IN s ==> P(IMAGE f t))) /\
1906      (!f s P. (?y. y IN IMAGE f s /\ P y) <=> (?x. x IN s /\ P(f x))) /\
1907      (!f s P. (?u. u IN IMAGE (IMAGE f) s /\ P u) <=>
1908               (?t. t IN s /\ P(IMAGE f t)))`,
1909     SET_TAC[]) in
1910   fun avoid (asl,w as gl) ->
1911     let f,wff = dest_forall w in
1912     let vs,bod = strip_forall wff in
1913     let ant,cons = dest_imp bod in
1914     let hths = CONJUNCTS(ASSUME ant) in
1915     let fths = hths @ mapfilter (USABLE_CONCLUSION f hths) sths in
1916     let cths = mapfilter (USABLE_CONCLUSION f fths) [cth0; cth1]
1917     and vconv =
1918       try let vth = USABLE_CONCLUSION f fths QUANTIFY_SURJECTION_HIGHER_THM in
1919           PROVE_HYP vth o PARTIAL_EXPAND_QUANTS_CONV avoid (ASSUME(concl vth))
1920       with Failure _ -> ALL_CONV
1921     and bths = LINEAR_INVARIANTS f fths in
1922     (MAP_EVERY X_GEN_TAC (f::vs) THEN DISCH_TAC THEN
1923      GEN_REWRITE_TAC (LAND_CONV o ONCE_DEPTH_CONV) cths THEN
1924      CONV_TAC(LAND_CONV vconv) THEN
1925      GEN_REWRITE_TAC (LAND_CONV o REDEPTH_CONV) [bths] THEN
1926      REWRITE_TAC[] THEN
1927      REWRITE_TAC(mapfilter (ADD_ASSUM ant o ISPEC f) aths) THEN
1928      GEN_REWRITE_TAC (LAND_CONV o REDEPTH_CONV) [bths] THEN
1929      REWRITE_TAC[]) gl;;
1930
1931 (* ------------------------------------------------------------------------- *)
1932 (* Scale so that a chosen vector has size 1. Generates a conjunction of      *)
1933 (* two formulas, one for the zero case (which one hopes is trivial) and      *)
1934 (* one just like the original goal but with a norm(...) = 1 assumption.      *)
1935 (* ------------------------------------------------------------------------- *)
1936
1937 let GEOM_NORMALIZE_RULE =
1938   let pth = prove
1939    (`!a:real^N. ~(a = vec 0)
1940                 ==> vec 0 = norm(a) % vec 0 /\
1941                     a = norm(a) % inv(norm a) % a /\
1942                     {} = IMAGE (\x. norm(a) % x) {} /\
1943                     {} = IMAGE (IMAGE (\x. norm(a) % x)) {} /\
1944                     (:real^N) = IMAGE (\x. norm(a) % x) (:real^N) /\
1945                     (:real^N->bool) =
1946                     IMAGE (IMAGE (\x. norm(a) % x)) (:real^N->bool) /\
1947                     [] = MAP (\x. norm(a) % x) []`,
1948     REWRITE_TAC[IMAGE_CLAUSES; VECTOR_MUL_ASSOC; VECTOR_MUL_RZERO; MAP] THEN
1949     SIMP_TAC[NORM_EQ_0; REAL_MUL_RINV; VECTOR_MUL_LID] THEN
1950     GEN_TAC THEN DISCH_TAC THEN
1951     REWRITE_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> !y. ?x. f x = y`] THEN
1952     ASM_REWRITE_TAC[SURJECTIVE_IMAGE] THEN
1953     X_GEN_TAC `y:real^N` THEN EXISTS_TAC `inv(norm(a:real^N)) % y:real^N` THEN
1954     ASM_SIMP_TAC[VECTOR_MUL_ASSOC; NORM_EQ_0; REAL_MUL_RINV; VECTOR_MUL_LID])
1955   and qth = prove
1956    (`!a:real^N.
1957         ~(a = vec 0)
1958         ==> ((!P. (!r:real. P r) <=> (!r. P(norm a * r))) /\
1959              (!P. (?r:real. P r) <=> (?r. P(norm a * r))) /\
1960              (!P. (!x:real^N. P x) <=> (!x. P (norm(a) % x))) /\
1961              (!P. (?x:real^N. P x) <=> (?x. P (norm(a) % x))) /\
1962              (!Q. (!s:real^N->bool. Q s) <=>
1963                   (!s. Q(IMAGE (\x. norm(a) % x) s))) /\
1964              (!Q. (?s:real^N->bool. Q s) <=>
1965                   (?s. Q(IMAGE (\x. norm(a) % x) s))) /\
1966              (!Q. (!s:(real^N->bool)->bool. Q s) <=>
1967                   (!s. Q(IMAGE (IMAGE (\x. norm(a) % x)) s))) /\
1968              (!Q. (?s:(real^N->bool)->bool. Q s) <=>
1969                   (?s. Q(IMAGE (IMAGE (\x. norm(a) % x)) s))) /\
1970              (!P. (!g:real^1->real^N. P g) <=>
1971                   (!g. P ((\x. norm(a) % x) o g))) /\
1972              (!P. (?g:real^1->real^N. P g) <=>
1973                   (?g. P ((\x. norm(a) % x) o g))) /\
1974              (!P. (!g:num->real^N. P g) <=>
1975                   (!g. P ((\x. norm(a) % x) o g))) /\
1976              (!P. (?g:num->real^N. P g) <=>
1977                   (?g. P ((\x. norm(a) % x) o g))) /\
1978              (!Q. (!l. Q l) <=> (!l. Q(MAP (\x:real^N. norm(a) % x) l))) /\
1979              (!Q. (?l. Q l) <=> (?l. Q(MAP (\x:real^N. norm(a) % x) l)))) /\
1980             ((!P. {x:real^N | P x} =
1981                   IMAGE (\x. norm(a) % x) {x | P(norm(a) % x)}) /\
1982              (!Q. {s:real^N->bool | Q s} =
1983                   IMAGE (IMAGE (\x. norm(a) % x))
1984                        {s | Q(IMAGE (\x. norm(a) % x) s)}) /\
1985              (!R. {l:(real^N)list | R l} =
1986                   IMAGE (MAP (\x:real^N. norm(a) % x))
1987                         {l | R(MAP (\x:real^N. norm(a) % x) l)}))`,
1988     GEN_TAC THEN DISCH_TAC THEN MATCH_MP_TAC(TAUT
1989      `(a /\ b) /\ c /\ d ==> (a /\ b /\ c) /\ d`) THEN
1990     CONJ_TAC THENL
1991      [ASM_MESON_TAC[NORM_EQ_0; REAL_FIELD `~(x = &0) ==> x * inv x * a = a`];
1992       MP_TAC(ISPEC `\x:real^N. norm(a:real^N) % x`
1993         (INST_TYPE [`:real^1`,`:C`] QUANTIFY_SURJECTION_HIGHER_THM)) THEN
1994       ASM_REWRITE_TAC[] THEN DISCH_THEN MATCH_MP_TAC THEN
1995       ASM_SIMP_TAC[SURJECTIVE_SCALING; NORM_EQ_0]])
1996   and lth = prove
1997    (`(!b:real^N. ~(b = vec 0) ==> (P(b) <=> Q(inv(norm b) % b)))
1998      ==> P(vec 0) /\ (!b. norm(b) = &1 ==> Q b) ==> (!b. P b)`,
1999     REPEAT STRIP_TAC THEN
2000     ASM_CASES_TAC `b:real^N = vec 0` THEN ASM_SIMP_TAC[] THEN
2001     FIRST_X_ASSUM MATCH_MP_TAC THEN
2002     ASM_SIMP_TAC[NORM_MUL; REAL_ABS_INV; REAL_ABS_NORM;
2003                  REAL_MUL_LINV; NORM_EQ_0]) in
2004   fun avoid tm ->
2005     let x,tm0 = dest_forall tm in
2006     let cth = UNDISCH(ISPEC x pth)
2007     and vth = UNDISCH(ISPEC x qth) in
2008     let th1 = ONCE_REWRITE_CONV[cth] tm0 in
2009     let th2 = CONV_RULE (RAND_CONV
2010      (PARTIAL_EXPAND_QUANTS_CONV avoid vth)) th1 in
2011     let th3 = SCALING_THEOREMS x in
2012     let th3' = (end_itlist CONJ (map
2013        (fun th -> let avs,_ = strip_forall(concl th) in
2014                   let gvs = map (genvar o type_of) avs in
2015                   GENL gvs (SPECL gvs th))
2016        (CONJUNCTS th3))) in
2017     let th4 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV)
2018                [BETA_THM; th3'] th2 in
2019     MATCH_MP lth (GEN x (DISCH_ALL th4));;
2020
2021 let GEN_GEOM_NORMALIZE_TAC x avoid (asl,w as gl) =
2022   let avs,bod = strip_forall w
2023   and avs' = subtract (frees w) (freesl(map (concl o snd) asl)) in
2024   (MAP_EVERY X_GEN_TAC avs THEN
2025    MAP_EVERY (fun t -> SPEC_TAC(t,t)) (rev(subtract (avs@avs') [x])) THEN
2026    SPEC_TAC(x,x) THEN
2027    W(MATCH_MP_TAC o GEOM_NORMALIZE_RULE avoid o snd)) gl;;
2028
2029 let GEOM_NORMALIZE_TAC x = GEN_GEOM_NORMALIZE_TAC x [];;
2030
2031 (* ------------------------------------------------------------------------- *)
2032 (* Add invariance theorems for collinearity.                                 *)
2033 (* ------------------------------------------------------------------------- *)
2034
2035 let COLLINEAR_TRANSLATION_EQ = prove
2036  (`!a s. collinear (IMAGE (\x. a + x) s) <=> collinear s`,
2037   REWRITE_TAC[collinear] THEN GEOM_TRANSLATE_TAC["u"]);;
2038
2039 add_translation_invariants [COLLINEAR_TRANSLATION_EQ];;
2040
2041 let COLLINEAR_TRANSLATION = prove
2042  (`!s a. collinear s ==> collinear (IMAGE (\x. a + x) s)`,
2043   REWRITE_TAC[COLLINEAR_TRANSLATION_EQ]);;
2044
2045 let COLLINEAR_LINEAR_IMAGE = prove
2046  (`!f s. collinear s /\ linear f ==> collinear(IMAGE f s)`,
2047   REPEAT GEN_TAC THEN DISCH_THEN(CONJUNCTS_THEN2 MP_TAC ASSUME_TAC) THEN
2048   REWRITE_TAC[collinear; IMP_CONJ; RIGHT_FORALL_IMP_THM; FORALL_IN_IMAGE] THEN
2049   ASM_MESON_TAC[LINEAR_SUB; LINEAR_CMUL]);;
2050
2051 let COLLINEAR_LINEAR_IMAGE_EQ = prove
2052  (`!f s. linear f /\ (!x y. f x = f y ==> x = y)
2053          ==> (collinear (IMAGE f s) <=> collinear s)`,
2054   MATCH_ACCEPT_TAC(LINEAR_INVARIANT_RULE COLLINEAR_LINEAR_IMAGE));;
2055
2056 add_linear_invariants [COLLINEAR_LINEAR_IMAGE_EQ];;
2057
2058 (* ------------------------------------------------------------------------- *)
2059 (* Take a theorem "th" with outer universal quantifiers involving real^N     *)
2060 (* and a theorem "dth" asserting |- dimindex(:M) <= dimindex(:N) and         *)
2061 (* return a theorem replacing type :N by :M in th. Neither N or M need be a  *)
2062 (* type variable.                                                            *)
2063 (* ------------------------------------------------------------------------- *)
2064
2065 let GEOM_DROP_DIMENSION_RULE =
2066   let oth = prove
2067    (`!f:real^M->real^N.
2068           linear f /\ (!x. norm(f x) = norm x)
2069           ==> linear f /\
2070               (!x y. f x = f y ==> x = y) /\
2071               (!x. norm(f x) = norm x)`,
2072     MESON_TAC[PRESERVES_NORM_INJECTIVE])
2073   and cth = prove
2074    (`linear(f:real^M->real^N)
2075      ==> vec 0 = f(vec 0) /\
2076          {} = IMAGE f {} /\
2077          {} = IMAGE (IMAGE f) {} /\
2078          [] = MAP f []`,
2079     REWRITE_TAC[IMAGE_CLAUSES; MAP; GSYM LINEAR_0]) in
2080   fun dth th ->
2081     let ath = GEN_ALL th
2082     and eth = MATCH_MP ISOMETRY_UNIV_UNIV dth
2083     and avoid = variables(concl th) in
2084     let f,bod = dest_exists(concl eth) in
2085     let fimage = list_mk_icomb "IMAGE" [f]
2086     and fmap = list_mk_icomb "MAP" [f]
2087     and fcompose = list_mk_icomb "o" [f] in
2088     let fimage2 = list_mk_icomb "IMAGE" [fimage] in
2089     let lin,iso = CONJ_PAIR(ASSUME bod) in
2090     let olduniv = rand(rand(concl dth))
2091     and newuniv = rand(lhand(concl dth)) in
2092     let oldty = fst(dest_fun_ty(type_of olduniv))
2093     and newty = fst(dest_fun_ty(type_of newuniv)) in
2094     let newvar v =
2095        let n,t = dest_var v in
2096        variant avoid (mk_var(n,tysubst[newty,oldty] t)) in
2097     let newterm v =
2098       try let v' = newvar v in
2099           tryfind (fun f -> mk_comb(f,v')) [f;fimage;fmap;fcompose;fimage2]
2100       with Failure _ -> v in
2101     let specrule th =
2102       let v = fst(dest_forall(concl th)) in SPEC (newterm v) th in
2103     let sth = SUBS(CONJUNCTS(MATCH_MP cth lin)) ath in
2104     let fth = SUBS[SYM(MATCH_MP LINEAR_0 lin)] (repeat specrule sth) in
2105     let thps = CONJUNCTS(MATCH_MP oth (ASSUME bod)) in
2106     let th5 = LINEAR_INVARIANTS f thps in
2107     let th6 = GEN_REWRITE_RULE REDEPTH_CONV [th5] fth in
2108     let th7 = PROVE_HYP eth (SIMPLE_CHOOSE f th6) in
2109     GENL (map newvar (fst(strip_forall(concl ath)))) th7;;
2110
2111 (* ------------------------------------------------------------------------- *)
2112 (* Transfer theorems automatically between same-dimension spaces.            *)
2113 (* Given dth = A |- dimindex(:M) = dimindex(:N)                              *)
2114 (* and a theorem th involving variables of type real^N                       *)
2115 (* returns a corresponding theorem mapped to type real^M with assumptions A. *)
2116 (* ------------------------------------------------------------------------- *)
2117
2118 let GEOM_EQUAL_DIMENSION_RULE =
2119   let bth = prove
2120    (`dimindex(:M) = dimindex(:N)
2121      ==> ?f:real^M->real^N.
2122              (linear f /\ (!y. ?x. f x = y)) /\
2123              (!x. norm(f x) = norm x)`,
2124     REWRITE_TAC[SET_RULE `(!y. ?x. f x = y) <=> IMAGE f UNIV = UNIV`] THEN
2125     DISCH_TAC THEN REWRITE_TAC[GSYM CONJ_ASSOC] THEN
2126     MATCH_MP_TAC ISOMETRY_UNIV_SUBSPACE THEN
2127     REWRITE_TAC[SUBSPACE_UNIV; DIM_UNIV] THEN FIRST_ASSUM ACCEPT_TAC)
2128   and pth = prove
2129    (`!f:real^M->real^N.
2130         linear f /\ (!y. ?x. f x = y)
2131          ==> (vec 0 = f(vec 0) /\
2132               {} = IMAGE f {} /\
2133               {} = IMAGE (IMAGE f) {} /\
2134               (:real^N) = IMAGE f (:real^M) /\
2135               (:real^N->bool) = IMAGE (IMAGE f) (:real^M->bool) /\
2136               [] = MAP f []) /\
2137              ((!P. (!x. P x) <=> (!x. P (f x))) /\
2138               (!P. (?x. P x) <=> (?x. P (f x))) /\
2139               (!Q. (!s. Q s) <=> (!s. Q (IMAGE f s))) /\
2140               (!Q. (?s. Q s) <=> (?s. Q (IMAGE f s))) /\
2141               (!Q. (!s. Q s) <=> (!s. Q (IMAGE (IMAGE f) s))) /\
2142               (!Q. (?s. Q s) <=> (?s. Q (IMAGE (IMAGE f) s))) /\
2143               (!P. (!g:real^1->real^N. P g) <=> (!g. P (f o g))) /\
2144               (!P. (?g:real^1->real^N. P g) <=> (?g. P (f o g))) /\
2145               (!P. (!g:num->real^N. P g) <=> (!g. P (f o g))) /\
2146               (!P. (?g:num->real^N. P g) <=> (?g. P (f o g))) /\
2147               (!Q. (!l. Q l) <=> (!l. Q(MAP f l))) /\
2148               (!Q. (?l. Q l) <=> (?l. Q(MAP f l)))) /\
2149              ((!P. {x | P x} = IMAGE f {x | P(f x)}) /\
2150               (!Q. {s | Q s} = IMAGE (IMAGE f) {s | Q(IMAGE f s)}) /\
2151               (!R. {l | R l} = IMAGE (MAP f) {l | R(MAP f l)}))`,
2152     GEN_TAC THEN
2153     SIMP_TAC[SET_RULE `UNIV = IMAGE f UNIV <=> (!y. ?x. f x = y)`;
2154              SURJECTIVE_IMAGE] THEN
2155     MATCH_MP_TAC MONO_AND THEN
2156     REWRITE_TAC[QUANTIFY_SURJECTION_HIGHER_THM] THEN
2157     REWRITE_TAC[IMAGE_CLAUSES; MAP] THEN MESON_TAC[LINEAR_0]) in
2158   fun dth th ->
2159     let eth = EXISTS_GENVAR_RULE (MATCH_MP bth dth) in
2160     let f,bod = dest_exists(concl eth) in
2161     let lsth,neth = CONJ_PAIR(ASSUME bod) in
2162     let cth,qth = CONJ_PAIR(MATCH_MP pth lsth) in
2163     let th1 = CONV_RULE (EXPAND_QUANTS_CONV qth) th in
2164     let ith = LINEAR_INVARIANTS f (neth::CONJUNCTS lsth) in
2165     let th2 = GEN_REWRITE_RULE (RAND_CONV o REDEPTH_CONV) [BETA_THM;ith] th1 in
2166     PROVE_HYP eth (SIMPLE_CHOOSE f th2);;