1(* ------------------------------------------------------------------------- *) 2(* Binomial coefficients and expansion. *) 3(* ------------------------------------------------------------------------- *) 4 5(*===========================================================================*) 6 7(* add all dependent libraries for script *) 8open HolKernel boolLib bossLib Parse; 9 10(* declare new theory at start *) 11val _ = new_theory "binomial"; 12 13(* ------------------------------------------------------------------------- *) 14 15 16(* val _ = load "jcLib"; *) 17open jcLib; 18 19(* Get dependent theories in lib *) 20(* val _ = load "helperFunctionTheory"; *) 21(* (* val _ = load "helperNumTheory"; -- in helperFunctionTheory *) *) 22(* (* val _ = load "helperSetTheory"; -- in helperFunctionTheory *) *) 23open helperNumTheory helperSetTheory helperFunctionTheory; 24 25(* val _ = load "helperListTheory"; *) 26open helperListTheory; 27 28(* open dependent theories *) 29open pred_setTheory listTheory; 30(* (* val _ = load "dividesTheory"; -- in helperNumTheory *) *) 31(* (* val _ = load "gcdTheory"; -- in helperNumTheory *) *) 32open arithmeticTheory dividesTheory gcdTheory; 33 34 35(* ------------------------------------------------------------------------- *) 36(* Binomial scripts in HOL: 37C:\jc\www\ml\hol\info\Hol\examples\miller\RSA\summationScript.sml 38C:\jc\www\ml\hol\info\Hol\examples\miller\RSA\powerScript.sml 39C:\jc\www\ml\hol\info\Hol\examples\miller\RSA\binomialScript.sml 40*) 41(* ------------------------------------------------------------------------- *) 42 43(* ------------------------------------------------------------------------- *) 44(* Binomial Documentation *) 45(* ------------------------------------------------------------------------- *) 46(* Overloading: 47*) 48(* Definitions and Theorems (# are exported): 49 50 Binomial Coefficients: 51 binomial_def |- (binomial 0 0 = 1) /\ (!n. binomial (SUC n) 0 = 1) /\ 52 (!k. binomial 0 (SUC k) = 0) /\ 53 !n k. binomial (SUC n) (SUC k) = binomial n k + binomial n (SUC k) 54 binomial_alt |- !n k. binomial n 0 = 1 /\ binomial 0 (k + 1) = 0 /\ 55 binomial (n + 1) (k + 1) = binomial n k + binomial n (k + 1) 56 binomial_less_0 |- !n k. n < k ==> (binomial n k = 0) 57 binomial_n_0 |- !n. binomial n 0 = 1 58 binomial_n_n |- !n. binomial n n = 1 59 binomial_0_n |- !n. binomial 0 n = if n = 0 then 1 else 0 60 binomial_recurrence |- !n k. binomial (SUC n) (SUC k) = binomial n k + binomial n (SUC k) 61 binomial_formula |- !n k. binomial (n + k) k * (FACT n * FACT k) = FACT (n + k) 62 binomial_formula2 |- !n k. k <= n ==> (FACT n = binomial n k * (FACT (n - k) * FACT k)) 63 binomial_formula3 |- !n k. k <= n ==> (binomial n k = FACT n DIV (FACT k * FACT (n - k))) 64 binomial_fact |- !n k. k <= n ==> (binomial n k = FACT n DIV (FACT k * FACT (n - k))) 65 binomial_n_k |- !n k. k <= n ==> (binomial n k = FACT n DIV FACT k DIV FACT (n - k) 66 binomial_n_1 |- !n. binomial n 1 = n 67 binomial_sym |- !n k. k <= n ==> (binomial n k = binomial n (n - k)) 68 binomial_is_integer |- !n k. k <= n ==> (FACT k * FACT (n - k)) divides (FACT n) 69 binomial_pos |- !n k. k <= n ==> 0 < binomial n k 70 binomial_eq_0 |- !n k. (binomial n k = 0) <=> n < k 71 binomial_up_eqn |- !n. 0 < n ==> !k. n * binomial (n - 1) k = (n - k) * binomial n k 72 binomial_up |- !n. 0 < n ==> !k. binomial (n - 1) k = (n - k) * binomial n k DIV n 73 binomial_right_eqn |- !n. 0 < n ==> !k. (k + 1) * binomial n (k + 1) = (n - k) * binomial n k 74 binomial_right |- !n. 0 < n ==> !k. binomial n (k + 1) = (n - k) * binomial n k DIV (k + 1) 75 binomial_monotone |- !n k. k < HALF n ==> binomial n k < binomial n (k + 1) 76 binomial_max |- !n k. binomial n k <= binomial n (HALF n) 77 78 Primes and Binomial Coefficients: 79 prime_divides_binomials |- !n. prime n ==> 1 < n /\ !k. 0 < k /\ k < n ==> n divides (binomial n k) 80 prime_divides_binomials_alt |- !n k. prime n /\ 0 < k /\ k < n ==> n divides binomial n k 81 prime_divisor_property |- !n p. 1 < n /\ p < n /\ prime p /\ p divides n ==> ~(p divides (FACT (n - 1) DIV FACT (n - p))) 82 divides_binomials_imp_prime |- !n. 1 < n /\ (!k. 0 < k /\ k < n ==> n divides (binomial n k)) ==> prime n 83 prime_iff_divides_binomials |- !n. prime n <=> 1 < n /\ !k. 0 < k /\ k < n ==> n divides (binomial n k) 84 prime_iff_divides_binomials_alt 85 |- !n. prime n <=> 1 < n /\ !k. 0 < k /\ k < n ==> binomial n k MOD n = 0 86 87 Binomial Theorem: 88 GENLIST_binomial_index_shift |- !n x y. GENLIST ((\k. binomial n k * x ** SUC (n - k) * y ** k) o SUC) n = 89 GENLIST (\k. binomial n (SUC k) * x ** (n - k) * y ** SUC k) n 90 binomial_index_shift |- !n x y. (\k. binomial (SUC n) k * x ** (SUC n - k) * y ** k) o SUC = 91 (\k. binomial (SUC n) (SUC k) * x ** (n - k) * y ** SUC k) 92 binomial_term_merge_x |- !n x y. (\k. x * k) o (\k. binomial n k * x ** (n - k) * y ** k) = 93 (\k. binomial n k * x ** SUC (n - k) * y ** k) 94 binomial_term_merge_y |- !n x y. (\k. y * k) o (\k. binomial n k * x ** (n - k) * y ** k) = 95 (\k. binomial n k * x ** (n - k) * y ** SUC k) 96 binomial_thm |- !n x y. (x + y) ** n = SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** k) (SUC n)) 97 binomial_thm_alt |- !n x y. (x + y) ** n = SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** k) (n + 1)) 98 binomial_sum |- !n. SUM (GENLIST (binomial n) (SUC n)) = 2 ** n 99 binomial_sum_alt |- !n. SUM (GENLIST (binomial n) (n + 1)) = 2 ** n 100 101 Binomial Horizontal List: 102 binomial_horizontal_0 |- binomial_horizontal 0 = [1] 103 binomial_horizontal_len |- !n. LENGTH (binomial_horizontal n) = n + 1 104 binomial_horizontal_mem |- !n k. k < n + 1 ==> MEM (binomial n k) (binomial_horizontal n) 105 binomial_horizontal_mem_iff |- !n k. MEM (binomial n k) (binomial_horizontal n) <=> k <= n 106 binomial_horizontal_member |- !n x. MEM x (binomial_horizontal n) <=> ?k. k <= n /\ (x = binomial n k) 107 binomial_horizontal_element |- !n k. k <= n ==> (EL k (binomial_horizontal n) = binomial n k) 108 binomial_horizontal_pos |- !n. EVERY (\x. 0 < x) (binomial_horizontal n) 109 binomial_horizontal_pos_alt |- !n x. MEM x (binomial_horizontal n) ==> 0 < x 110 binomial_horizontal_sum |- !n. SUM (binomial_horizontal n) = 2 ** n 111 binomial_horizontal_max |- !n. MAX_LIST (binomial_horizontal n) = binomial n (HALF n) 112 binomial_row_max |- !n. MAX_SET (IMAGE (binomial n) (count (n + 1))) = binomial n (HALF n) 113 binomial_product_identity |- !m n k. k <= m /\ m <= n ==> 114 (binomial m k * binomial n m = binomial n k * binomial (n - k) (m - k)) 115 binomial_middle_upper_bound |- !n. binomial n (HALF n) <= 4 ** HALF n 116 117 Stirling's Approximation: 118 Stirling = (!n. FACT n = (SQRT (2 * pi * n)) * (n DIV e) ** n) /\ 119 (!n. SQRT n = n ** h) /\ (2 * h = 1) /\ (0 < pi) /\ (0 < e) /\ 120 (!a b x y. (a * b) DIV (x * y) = (a DIV x) * (b DIV y)) /\ 121 (!a b c. (a DIV c) DIV (b DIV c) = a DIV b) 122 binomial_middle_by_stirling |- Stirling ==> 123 !n. 0 < n /\ EVEN n ==> (binomial n (HALF n) = 2 ** (n + 1) DIV SQRT (2 * pi * n)) 124 125 Useful theorems for Binomial: 126 binomial_range_shift |- !n . 0 < n ==> ((!k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0)) <=> 127 (!h. h < PRE n ==> ((binomial n (SUC h)) MOD n = 0))) 128 binomial_mod_zero |- !n. 0 < n ==> !k. (binomial n k MOD n = 0) <=> 129 (!x y. (binomial n k * x ** (n-k) * y ** k) MOD n = 0) 130 binomial_range_shift_alt |- !n . 0 < n ==> ((!k. 0 < k /\ k < n ==> 131 (!x y. ((binomial n k * x ** (n - k) * y ** k) MOD n = 0))) <=> 132 (!h. h < PRE n ==> (!x y. ((binomial n (SUC h) * x ** (n - (SUC h)) * y ** (SUC h)) MOD n = 0)))) 133 binomial_mod_zero_alt |- !n. 0 < n ==> ((!k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0)) <=> 134 !x y. SUM (GENLIST ((\k. (binomial n k * x ** (n - k) * y ** k) MOD n) o SUC) (PRE n)) = 0) 135 136 Binomial Theorem with prime exponent: 137 binomial_thm_prime |- !p. prime p ==> (!x y. (x + y) ** p MOD p = (x ** p + y ** p) MOD p) 138*) 139 140(* ------------------------------------------------------------------------- *) 141(* Helper Theorems *) 142(* ------------------------------------------------------------------------- *) 143 144(* ------------------------------------------------------------------------- *) 145(* Binomial Coefficients *) 146(* ------------------------------------------------------------------------- *) 147 148(* Define Binomials: 149 C(n,0) = 1 150 C(0,k) = 0 if k > 0 151 C(n+1,k+1) = C(n,k) + C(n,k+1) 152*) 153val binomial_def = Define` 154 (binomial 0 0 = 1) /\ 155 (binomial (SUC n) 0 = 1) /\ 156 (binomial 0 (SUC k) = 0) /\ 157 (binomial (SUC n) (SUC k) = binomial n k + binomial n (SUC k)) 158`; 159 160(* Theorem: alternative definition of C(n,k). *) 161(* Proof: by binomial_def. *) 162Theorem binomial_alt: 163 !n k. (binomial n 0 = 1) /\ 164 (binomial 0 (k + 1) = 0) /\ 165 (binomial (n + 1) (k + 1) = binomial n k + binomial n (k + 1)) 166Proof 167 rewrite_tac[binomial_def, GSYM ADD1] >> 168 (Cases_on `n` >> simp[binomial_def]) 169QED 170 171(* Basic properties *) 172 173(* Theorem: C(n,k) = 0 if n < k *) 174(* Proof: 175 By induction on n. 176 Base case: C(0,k) = 0 if 0 < k, by definition. 177 Step case: assume C(n,k) = 0 if n < k. 178 then for SUC n < k, 179 C(SUC n, k) 180 = C(SUC n, SUC h) where k = SUC h 181 = C(n,h) + C(n,SUC h) h < SUC h = k 182 = 0 + 0 by induction hypothesis 183 = 0 184*) 185val binomial_less_0 = store_thm( 186 "binomial_less_0", 187 ``!n k. n < k ==> (binomial n k = 0)``, 188 Induct_on `n` >- 189 metis_tac[binomial_def, num_CASES, NOT_ZERO_LT_ZERO] >> 190 rw[binomial_def] >> 191 `?h. k = SUC h` by metis_tac[SUC_NOT, NOT_ZERO_LT_ZERO, LESS_EQ_SUC, LESS_TRANS] >> 192 metis_tac[binomial_def, LESS_MONO_EQ, LESS_TRANS, LESS_SUC, ADD_0]); 193 194(* Theorem: C(n,0) = 1 *) 195(* Proof: 196 If n = 0, C(n, 0) = C(0, 0) = 1 by binomial_def 197 If n <> 0, n = SUC m, and C(SUC m, 0) = 1 by binomial_def 198*) 199val binomial_n_0 = store_thm( 200 "binomial_n_0", 201 ``!n. binomial n 0 = 1``, 202 metis_tac[binomial_def, num_CASES]); 203 204(* Theorem: C(n,n) = 1 *) 205(* Proof: 206 By induction on n. 207 Base case: C(0,0) = 1, true by binomial_def. 208 Step case: assume C(n,n) = 1 209 C(SUC n, SUC n) 210 = C(n,n) + C(n,SUC n) 211 = 1 + C(n,SUC n) by induction hypothesis 212 = 1 + 0 by binomial_less_0 213 = 1 214*) 215val binomial_n_n = store_thm( 216 "binomial_n_n", 217 ``!n. binomial n n = 1``, 218 Induct_on `n` >- 219 metis_tac[binomial_def] >> 220 metis_tac[binomial_def, LESS_SUC, binomial_less_0, ADD_0]); 221 222(* Theorem: binomial 0 n = if n = 0 then 1 else 0 *) 223(* Proof: 224 If n = 0, 225 binomial 0 0 = 1 by binomial_n_0 226 If n <> 0, then 0 < n. 227 binomial 0 n = 0 by binomial_less_0 228*) 229val binomial_0_n = store_thm( 230 "binomial_0_n", 231 ``!n. binomial 0 n = if n = 0 then 1 else 0``, 232 rw[binomial_n_0, binomial_less_0]); 233 234(* Theorem: C(n+1,k+1) = C(n,k) + C(n,k+1) *) 235(* Proof: by definition. *) 236val binomial_recurrence = store_thm( 237 "binomial_recurrence", 238 ``!n k. binomial (SUC n) (SUC k) = binomial n k + binomial n (SUC k)``, 239 rw[binomial_def]); 240 241(* Theorem: C(n+k,k) = (n+k)!/n!k! *) 242(* Proof: 243 By induction on k. 244 Base case: C(n,0) = n!n! = 1 by binomial_n_0 245 Step case: assume C(n+k,k) = (n+k)!/n!k! 246 To prove C(n+SUC k, SUC k) = (n+SUC k)!/n!(SUC k)! 247 By induction on n. 248 Base case: C(SUC k, SUC k) = (SUC k)!/(SUC k)! = 1 by binomial_n_n 249 Step case: assume C(n+SUC k, SUC k) = (n +SUC k)!/n!(SUC k)! 250 To prove C(SUC n + SUC k, SUC k) = (SUC n + SUC k)!/(SUC n)!(SUC k)! 251 C(SUC n + SUC k, SUC k) 252 = C(SUC SUC (n+k), SUC k) 253 = C(SUC (n+k),k) + C(SUC (n+k), SUC k) 254 = C(SUC n + k, k) + C(n + SUC k, SUC k) 255 = (SUC n + k)!/(SUC n)!k! + (n + SUC k)!/n!(SUC k)! by two induction hypothesis 256 = ((SUC n + k)!(SUC k) + (n + SUC k)(SUC n))/(SUC n)!(SUC k)! 257 = (SUC n + SUC k)!/(SUC n)!(SUC k)! 258*) 259val binomial_formula = store_thm( 260 "binomial_formula", 261 ``!n k. binomial (n+k) k * (FACT n * FACT k) = FACT (n+k)``, 262 Induct_on `k` >- 263 metis_tac[binomial_n_0, FACT, MULT_CLAUSES, ADD_0] >> 264 Induct_on `n` >- 265 metis_tac[binomial_n_n, FACT, MULT_CLAUSES, ADD_CLAUSES] >> 266 `SUC n + SUC k = SUC (SUC (n+k))` by decide_tac >> 267 `SUC (n + k) = SUC n + k` by decide_tac >> 268 `binomial (SUC n + SUC k) (SUC k) * (FACT (SUC n) * FACT (SUC k)) = 269 (binomial (SUC (n + k)) k + 270 binomial (SUC (n + k)) (SUC k)) * (FACT (SUC n) * FACT (SUC k))` 271 by metis_tac[binomial_recurrence] >> 272 `_ = binomial (SUC (n + k)) k * (FACT (SUC n) * FACT (SUC k)) + 273 binomial (SUC (n + k)) (SUC k) * (FACT (SUC n) * FACT (SUC k))` 274 by metis_tac[RIGHT_ADD_DISTRIB] >> 275 `_ = binomial (SUC n + k) k * (FACT (SUC n) * ((SUC k) * FACT k)) + 276 binomial (n + SUC k) (SUC k) * ((SUC n) * FACT n * FACT (SUC k))` 277 by metis_tac[ADD_COMM, SUC_ADD_SYM, FACT] >> 278 `_ = binomial (SUC n + k) k * FACT (SUC n) * FACT k * (SUC k) + 279 binomial (n + SUC k) (SUC k) * FACT n * FACT (SUC k) * (SUC n)` 280 by metis_tac[MULT_COMM, MULT_ASSOC] >> 281 `_ = FACT (SUC n + k) * SUC k + FACT (n + SUC k) * SUC n` 282 by metis_tac[MULT_COMM, MULT_ASSOC] >> 283 `_ = FACT (SUC (n+k)) * SUC k + FACT (SUC (n+k)) * SUC n` 284 by metis_tac[ADD_COMM, SUC_ADD_SYM] >> 285 `_ = FACT (SUC (n+k)) * (SUC k + SUC n)` by metis_tac[LEFT_ADD_DISTRIB] >> 286 `_ = (SUC n + SUC k) * FACT (SUC (n+k))` by metis_tac[MULT_COMM, ADD_COMM] >> 287 metis_tac[FACT]); 288 289(* Theorem: C(n,k) = n!/k!(n-k)! for 0 <= k <= n *) 290(* Proof: 291 FACT n 292 = FACT ((n-k)+k) by SUB_ADD, k <= n. 293 = binomial ((n-k)+k) k * (FACT (n-k) * FACT k) by binomial_formula 294 = binomial n k * (FACT (n-k) * FACT k)) by SUB_ADD, k <= n. 295*) 296val binomial_formula2 = store_thm( 297 "binomial_formula2", 298 ``!n k. k <= n ==> (FACT n = binomial n k * (FACT (n-k) * FACT k))``, 299 metis_tac[binomial_formula, SUB_ADD]); 300 301(* Theorem: k <= n ==> binomial n k = (FACT n) DIV ((FACT k) * (FACT (n - k))) *) 302(* Proof: 303 binomial n k 304 = (binomial n k * (FACT (n - k) * FACT k)) DIV ((FACT (n - k) * FACT k)) by MULT_DIV 305 = (FACT n) DIV ((FACT (n - k) * FACT k)) by binomial_formula2 306 = (FACT n) DIV ((FACT k * FACT (n - k))) by MULT_COMM 307*) 308val binomial_formula3 = store_thm( 309 "binomial_formula3", 310 ``!n k. k <= n ==> (binomial n k = (FACT n) DIV ((FACT k) * (FACT (n - k))))``, 311 metis_tac[binomial_formula2, MULT_COMM, MULT_DIV, MULT_EQ_0, FACT_LESS, NOT_ZERO_LT_ZERO]); 312 313(* Theorem alias. *) 314val binomial_fact = save_thm("binomial_fact", binomial_formula3); 315(* val binomial_fact = |- !n k. k <= n ==> (binomial n k = FACT n DIV (FACT k * FACT (n - k))): thm *) 316 317(* Theorem: k <= n ==> binomial n k = (FACT n) DIV (FACT k) DIV (FACT (n - k)) *) 318(* Proof: 319 binomial n k 320 = (FACT n) DIV ((FACT k * FACT (n - k))) by binomial_formula3 321 = (FACT n) DIV (FACT k) DIV (FACT (n - k)) by DIV_DIV_DIV_MULT 322*) 323val binomial_n_k = store_thm( 324 "binomial_n_k", 325 ``!n k. k <= n ==> (binomial n k = (FACT n) DIV (FACT k) DIV (FACT (n - k)))``, 326 metis_tac[DIV_DIV_DIV_MULT, binomial_formula3, MULT_EQ_0, FACT_LESS, NOT_ZERO_LT_ZERO]); 327 328(* Theorem: binomial n 1 = n *) 329(* Proof: 330 If n = 0, 331 binomial 0 1 332 = if 1 = 0 then 1 else 0 by binomial_0_n 333 = 0 by 1 = 0 = F 334 If n <> 0, then 0 < n. 335 Thus 1 <= n, and n = SUC (n-1) by 0 < n 336 binomial n 1 337 = FACT n DIV FACT 1 DIV FACT (n - 1) by binomial_n_k, 1 <= n 338 = FACT n DIV 1 DIV (FACT (n-1)) by FACT, ONE 339 = FACT n DIV (FACT (n-1)) by DIV_1 340 = (n * FACT (n-1)) DIV (FACT (n-1)) by FACT 341 = n by MULT_DIV, FACT_LESS 342*) 343val binomial_n_1 = store_thm( 344 "binomial_n_1", 345 ``!n. binomial n 1 = n``, 346 rpt strip_tac >> 347 Cases_on `n = 0` >- 348 rw[binomial_0_n] >> 349 `1 <= n /\ (n = SUC (n-1))` by decide_tac >> 350 `binomial n 1 = FACT n DIV FACT 1 DIV FACT (n - 1)` by rw[binomial_n_k] >> 351 `_ = FACT n DIV 1 DIV (FACT (n-1))` by EVAL_TAC >> 352 `_ = FACT n DIV (FACT (n-1))` by rw[] >> 353 `_ = (n * FACT (n-1)) DIV (FACT (n-1))` by metis_tac[FACT] >> 354 `_ = n` by rw[MULT_DIV, FACT_LESS] >> 355 rw[]); 356 357(* Theorem: k <= n ==> (binomial n k = binomial n (n-k)) *) 358(* Proof: 359 Note (n-k) <= n always. 360 binomial n k 361 = (FACT n) DIV (FACT k * FACT (n - k)) by binomial_formula3, k <= n. 362 = (FACT n) DIV (FACT (n - k) * FACT k) by MULT_COMM 363 = (FACT n) DIV (FACT (n - k) * FACT (n-(n-k))) by n - (n-k) = k 364 = binomial n (n-k) by binomial_formula3, (n-k) <= n. 365*) 366val binomial_sym = store_thm( 367 "binomial_sym", 368 ``!n k. k <= n ==> (binomial n k = binomial n (n-k))``, 369 rpt strip_tac >> 370 `n - (n-k) = k` by decide_tac >> 371 `(n-k) <= n` by decide_tac >> 372 rw[binomial_formula3, MULT_COMM]); 373 374(* Theorem: k <= n ==> (FACT k * FACT (n-k)) divides (FACT n) *) 375(* Proof: 376 Since FACT n = binomial n k * (FACT (n - k) * FACT k) by binomial_formula2 377 = binomial n k * (FACT k * FACT (n - k)) by MULT_COMM 378 Hence (FACT k * FACT (n-k)) divides (FACT n) by divides_def 379*) 380val binomial_is_integer = store_thm( 381 "binomial_is_integer", 382 ``!n k. k <= n ==> (FACT k * FACT (n-k)) divides (FACT n)``, 383 metis_tac[binomial_formula2, MULT_COMM, divides_def]); 384 385(* Theorem: k <= n ==> 0 < binomial n k *) 386(* Proof: 387 Since FACT n = binomial n k * (FACT (n - k) * FACT k) by binomial_formula2 388 and 0 < FACT n, 0 < FACT (n-k), 0 < FACT k by FACT_LESS 389 Hence 0 < binomial n k by ZERO_LESS_MULT 390*) 391val binomial_pos = store_thm( 392 "binomial_pos", 393 ``!n k. k <= n ==> 0 < binomial n k``, 394 metis_tac[binomial_formula2, FACT_LESS, ZERO_LESS_MULT]); 395 396(* Theorem: (binomial n k = 0) <=> n < k *) 397(* Proof: 398 If part: (binomial n k = 0) ==> n < k 399 By contradiction, suppose k <= n. 400 Then 0 < binomial n k by binomial_pos 401 This contradicts binomial n k = 0 by NOT_ZERO_LT_ZERO 402 Only-if part: n < k ==> (binomial n k = 0) 403 This is true by binomial_less_0 404*) 405val binomial_eq_0 = store_thm( 406 "binomial_eq_0", 407 ``!n k. (binomial n k = 0) <=> n < k``, 408 rw[EQ_IMP_THM] >| [ 409 spose_not_then strip_assume_tac >> 410 `k <= n` by decide_tac >> 411 metis_tac[binomial_pos, NOT_ZERO_LT_ZERO], 412 rw[binomial_less_0] 413 ]); 414 415(* Relating Binomial to its up-entry: 416 417 binomial n k = (n, k, n-k) = n! / k! (n-k)! 418 binomial (n-1) k = (n-1, k, n-1-k) = (n-1)! / k! (n-1-k)! 419 = (n!/n) / k! ((n-k)!/(n-k)) 420 = (n-k) * binomial n k / n 421*) 422 423(* Theorem: 0 < n ==> !k. n * binomial (n-1) k = (n-k) * (binomial n k) *) 424(* Proof: 425 If n <= k, that is n-1 < k. 426 So binomial (n-1) k = 0 by binomial_less_0 427 and n - k = 0 by arithmetic 428 Hence true by MULT_EQ_0 429 Otherwise k < n, 430 or k <= n, 1 <= n-k, k <= n-1 431 Therefore, 432 FACT n = binomial n k * (FACT (n - k) * FACT k) by binomial_formula2, k <= n. 433 = binomial n k * ((n - k) * FACT (n-1-k) * FACT k) by FACT 434 = binomial n k * (n - k) * (FACT (n-1-k) * FACT k) by MULT_ASSOC 435 = (n - k) * binomial n k * (FACT (n-1-k) * FACT k) by MULT_COMM 436 FACT n = n * FACT (n-1) by FACT 437 = n * (binomial (n-1) k * (FACT (n-1-k) * FACT k)) by binomial_formula2, k <= n-1. 438 = (n * binomial (n-1) k) * (FACT (n-1-k) * FACT k) by MULT_ASSOC 439 Since 0 < FACT (n-1-k) * FACT k by FACT_LESS, MULT_EQ_0 440 n * binomial (n-1) k = (n-k) * (binomial n k) by MULT_RIGHT_CANCEL 441*) 442val binomial_up_eqn = store_thm( 443 "binomial_up_eqn", 444 ``!n. 0 < n ==> !k. n * binomial (n-1) k = (n-k) * (binomial n k)``, 445 rpt strip_tac >> 446 `!n. n <> 0 <=> 0 < n` by decide_tac >> 447 Cases_on `n <= k` >| [ 448 `n-1 < k /\ (n - k = 0)` by decide_tac >> 449 `binomial (n - 1) k = 0` by rw[binomial_less_0] >> 450 metis_tac[MULT_EQ_0], 451 `k < n /\ k <= n /\ 1 <= n-k /\ k <= n-1` by decide_tac >> 452 `SUC (n-1) = n` by decide_tac >> 453 `SUC (n-1-k) = n - k` by metis_tac[SUB_PLUS, ADD_COMM, ADD1, SUB_ADD] >> 454 `FACT n = binomial n k * (FACT (n - k) * FACT k)` by rw[binomial_formula2] >> 455 `_ = binomial n k * ((n - k) * FACT (n-1-k) * FACT k)` by metis_tac[FACT] >> 456 `_ = binomial n k * (n - k) * (FACT (n-1-k) * FACT k)` by rw[MULT_ASSOC] >> 457 `_ = (n - k) * binomial n k * (FACT (n-1-k) * FACT k)` by rw_tac std_ss[MULT_COMM] >> 458 `FACT n = n * FACT (n-1)` by metis_tac[FACT] >> 459 `_ = n * (binomial (n-1) k * (FACT (n-1-k) * FACT k))` by rw_tac std_ss[GSYM binomial_formula2] >> 460 `_ = (n * binomial (n-1) k) * (FACT (n-1-k) * FACT k)` by rw[MULT_ASSOC] >> 461 metis_tac[FACT_LESS, MULT_EQ_0, MULT_RIGHT_CANCEL] 462 ]); 463 464(* Theorem: 0 < n ==> !k. binomial (n-1) k = ((n-k) * (binomial n k)) DIV n *) 465(* Proof: 466 Since n * binomial (n-1) k = (n-k) * (binomial n k) by binomial_up_eqn 467 binomial (n-1) k = (n-k) * (binomial n k) DIV n by DIV_SOLVE, 0 < n. 468*) 469val binomial_up = store_thm( 470 "binomial_up", 471 ``!n. 0 < n ==> !k. binomial (n-1) k = ((n-k) * (binomial n k)) DIV n``, 472 rw[binomial_up_eqn, DIV_SOLVE]); 473 474(* Relating Binomial to its right-entry: 475 476 binomial n k = (n, k, n-k) = n! / k! (n-k)! 477 binomial n (k+1) = (n, k+1, n-k-1) = n! / (k+1)! (n-k-1)! 478 = n! / (k+1) * k! ((n-k)!/(n-k)) 479 = (n-k) * binomial n k / (k+1) 480*) 481 482(* Theorem: 0 < n ==> !k. (k + 1) * binomial n (k+1) = (n - k) * binomial n k *) 483(* Proof: 484 If n <= k, that is n < k+1. 485 So binomial n (k+1) = 0 by binomial_less_0 486 and n - k = 0 by arithmetic 487 Hence true by MULT_EQ_0 488 Otherwise k < n, 489 or k <= n, 1 <= n-k, k+1 <= n 490 Therefore, 491 FACT n = binomial n k * (FACT (n - k) * FACT k) by binomial_formula2, k <= n. 492 = binomial n k * ((n - k) * FACT (n-1-k) * FACT k) by FACT 493 = binomial n k * (n - k) * (FACT (n-1-k) * FACT k) by MULT_ASSOC 494 = (n - k) * binomial n k * (FACT (n-1-k) * FACT k) by MULT_COMM 495 FACT n = binomial n (k+1) * (FACT (n-(k+1)) * FACT (k+1)) by binomial_formula2, k+1 <= n. 496 = binomial n (k+1) * (FACT (n-1-k) * FACT (k+1)) by SUB_PLUS, ADD_COMM 497 = binomial n (k+1) * (FACT (n-1-k) * ((k+1) * FACT k)) by FACT 498 = binomial n (k+1) * ((k+1) * (FACT (n-1-k) * FACT k)) by MULT_ASSOC, MULT_COMM 499 = (k+1) * binomial n (k+1) * (FACT (n-1-k) * FACT k) by MULT_COMM, MULT_ASSOC 500 Since 0 < FACT (n-1-k) * FACT k by FACT_LESS, MULT_EQ_0 501 (k+1) * binomial n (k+1) = (n-k) * (binomial n k) by MULT_RIGHT_CANCEL 502*) 503val binomial_right_eqn = store_thm( 504 "binomial_right_eqn", 505 ``!n. 0 < n ==> !k. (k + 1) * binomial n (k+1) = (n - k) * binomial n k``, 506 rpt strip_tac >> 507 `!n. n <> 0 <=> 0 < n` by decide_tac >> 508 Cases_on `n <= k` >| [ 509 `n < k+1` by decide_tac >> 510 `binomial n (k+1) = 0` by rw[binomial_less_0] >> 511 `n - k = 0` by decide_tac >> 512 metis_tac[MULT_EQ_0], 513 `k < n /\ k <= n /\ 1 <= n-k /\ k+1 <= n` by decide_tac >> 514 `SUC k = k + 1` by decide_tac >> 515 `SUC (n-1-k) = n - k` by metis_tac[SUB_PLUS, ADD_COMM, ADD1, SUB_ADD] >> 516 `FACT n = binomial n k * (FACT (n - k) * FACT k)` by rw[binomial_formula2] >> 517 `_ = binomial n k * ((n - k) * FACT (n-1-k) * FACT k)` by metis_tac[FACT] >> 518 `_ = binomial n k * (n - k) * (FACT (n-1-k) * FACT k)` by rw[MULT_ASSOC] >> 519 `_ = (n - k) * binomial n k * (FACT (n-1-k) * FACT k)` by rw_tac std_ss[MULT_COMM] >> 520 `FACT n = binomial n (k+1) * (FACT (n-(k+1)) * FACT (k+1))` by rw[binomial_formula2] >> 521 `_ = binomial n (k+1) * (FACT (n-1-k) * FACT (k+1))` by metis_tac[SUB_PLUS, ADD_COMM] >> 522 `_ = binomial n (k+1) * (FACT (n-1-k) * ((k+1) * FACT k))` by metis_tac[FACT] >> 523 `_ = binomial n (k+1) * ((FACT (n-1-k) * (k+1)) * FACT k)` by rw[MULT_ASSOC] >> 524 `_ = binomial n (k+1) * ((k+1) * (FACT (n-1-k)) * FACT k)` by rw_tac std_ss[MULT_COMM] >> 525 `_ = (binomial n (k+1) * (k+1)) * (FACT (n-1-k) * FACT k)` by rw[MULT_ASSOC] >> 526 `_ = (k+1) * binomial n (k+1) * (FACT (n-1-k) * FACT k)` by rw_tac std_ss[MULT_COMM] >> 527 metis_tac[FACT_LESS, MULT_EQ_0, MULT_RIGHT_CANCEL] 528 ]); 529 530(* Theorem: 0 < n ==> !k. binomial n (k+1) = (n - k) * binomial n k DIV (k+1) *) 531(* Proof: 532 Since (k + 1) * binomial n (k+1) = (n - k) * binomial n k by binomial_right_eqn 533 binomial n (k+1) = (n - k) * binomial n k DIV (k+1) by DIV_SOLVE, 0 < k+1. 534*) 535val binomial_right = store_thm( 536 "binomial_right", 537 ``!n. 0 < n ==> !k. binomial n (k+1) = (n - k) * binomial n k DIV (k+1)``, 538 rw[binomial_right_eqn, DIV_SOLVE, DECIDE ``!k. 0 < k+1``]); 539 540(* 541 k < HALF n <=> k + 1 <= n - k 542n = 5, HALF n = 2, binomial 5 k: 1, 5, 10, 10, 5, 1 543 k= 0, 1, 2, 3, 4, 5 544 k < 2 <=> k + 1 <= 5 - k 545 k = 0 1 <= 5 binomial 5 1 >= binomial 5 0 546 k = 1 2 <= 4 binomial 5 2 >= binomial 5 1 547n = 6, HALF n = 3, binomial 6 k: 1, 6, 15, 20, 15, 6, 1 548 k= 0, 1, 2, 3, 4, 5, 6 549 k < 3 <=> k + 1 <= 6 - k 550 k = 0 1 <= 6 binomial 6 1 >= binomial 6 0 551 k = 1 2 <= 5 binomial 6 2 >= binomial 6 1 552 k = 2 3 <= 4 binomial 6 3 >= binomial 6 2 553*) 554 555(* Theorem: k < HALF n ==> binomial n k < binomial n (k + 1) *) 556(* Proof: 557 Note k < HALF n ==> 0 < n by ZERO_DIV, 0 < 2 558 also k < HALF n ==> k + 1 < n - k by LESS_HALF_IFF 559 so 0 < k + 1 /\ 0 < n - k by arithmetic 560 Now (k + 1) * binomial n (k + 1) = (n - k) * binomial n k by binomial_right_eqn, 0 < n 561 Note HALF n <= n by DIV_LESS_EQ, 0 < 2 562 so k < HALF n <= n by above 563 Thus 0 < binomial n k by binomial_pos, k <= n 564 and 0 < binomial n (k + 1) by MULT_0, MULT_EQ_0 565 Hence binomial n k < binomial n (k + 1) by MULT_EQ_LESS_TO_MORE 566*) 567val binomial_monotone = store_thm( 568 "binomial_monotone", 569 ``!n k. k < HALF n ==> binomial n k < binomial n (k + 1)``, 570 rpt strip_tac >> 571 `k + 1 < n - k` by rw[GSYM LESS_HALF_IFF] >> 572 `0 < k + 1 /\ 0 < n - k` by decide_tac >> 573 `(k + 1) * binomial n (k + 1) = (n - k) * binomial n k` by rw[binomial_right_eqn] >> 574 `HALF n <= n` by rw[DIV_LESS_EQ] >> 575 `0 < binomial n k` by rw[binomial_pos] >> 576 `0 < binomial n (k + 1)` by metis_tac[MULT_0, MULT_EQ_0, NOT_ZERO_LT_ZERO] >> 577 metis_tac[MULT_EQ_LESS_TO_MORE]); 578 579(* Theorem: binomial n k <= binomial n (HALF n) *) 580(* Proof: 581 Since (k + 1) * binomial n (k + 1) = (n - k) * binomial n k by binomial_right_eqn 582 binomial n (k + 1) / binomial n k = (n - k) / (k + 1) 583 As k varies from 0, 1, to (n-1), n 584 the ratio varies from n/1, (n-1)/2, (n-2)/3, ...., 1/n, 0/(n+1). 585 The ratio is greater than 1 when (n - k) / (k + 1) > 1 586 or n - k > k + 1 587 or n > 2 * k + 1 588 or HALF n >= k + (HALF 1) 589 or k <= HALF n 590 Thus (binomial n (HALF n)) is greater than all preceding coefficients. 591 For k > HALF n, note that (binomial n k = binomial n (n - k)) by binomial_sym 592 Hence (binomial n (HALF n)) is greater than all succeeding coefficients, too. 593 594 If n = 0, 595 binomial 0 k = 1 or 0 by binomial_0_n 596 binomial 0 (HALF 0) = 1 by binomial_0_n, ZERO_DIV 597 Hence true. 598 If n <> 0, 599 If k = HALF n, trivially true. 600 If k < HALF n, 601 Then binomial n k < binomial n (HALF n) by binomial_monotone, MONOTONE_MAX 602 Hence true. 603 If ~(k < HALF n), HALF n < k. 604 Then n - k <= HALF n by MORE_HALF_IMP 605 If k > n, 606 Then binomial n k = 0, hence true by binomial_less_0 607 If ~(k > n), then k <= n. 608 Then binomial n k = binomial n (n - k) by binomial_sym, k <= n 609 If n - k = HALF n, trivially true. 610 Otherwise, n - k < HALF n, 611 Thus binomial n (n - k) < binomial n (HALF n) by binomial_monotone, MONOTONE_MAX 612 Hence true. 613*) 614val binomial_max = store_thm( 615 "binomial_max", 616 ``!n k. binomial n k <= binomial n (HALF n)``, 617 rpt strip_tac >> 618 Cases_on `n = 0` >- 619 rw[binomial_0_n] >> 620 Cases_on `k = HALF n` >- 621 rw[] >> 622 Cases_on `k < HALF n` >| [ 623 `binomial n k < binomial n (HALF n)` by rw[binomial_monotone, MONOTONE_MAX] >> 624 decide_tac, 625 `HALF n < k` by decide_tac >> 626 `n - k <= HALF n` by rw[MORE_HALF_IMP] >> 627 Cases_on `k > n` >- 628 rw[binomial_less_0] >> 629 `k <= n` by decide_tac >> 630 `binomial n k = binomial n (n - k)` by rw[GSYM binomial_sym] >> 631 Cases_on `n - k = HALF n` >- 632 rw[] >> 633 `n - k < HALF n` by decide_tac >> 634 `binomial n (n - k) < binomial n (HALF n)` by rw[binomial_monotone, MONOTONE_MAX] >> 635 decide_tac 636 ]); 637 638(* ------------------------------------------------------------------------- *) 639(* Primes and Binomial Coefficients *) 640(* ------------------------------------------------------------------------- *) 641 642(* Theorem: n is prime ==> n divides C(n,k) for all 0 < k < n *) 643(* Proof: 644 C(n,k) = n!/k!/(n-k)! 645 or n! = C(n,k) k! (n-k)! 646 n divides n!, so n divides the product C(n,k) k!(n-k)! 647 For a prime n, n cannot divide k!(n-k)!, all factors less than prime n. 648 By Euclid's lemma, a prime divides a product must divide a factor. 649 So p divides C(n,k). 650*) 651val prime_divides_binomials = store_thm( 652 "prime_divides_binomials", 653 ``!n. prime n ==> 1 < n /\ (!k. 0 < k /\ k < n ==> n divides (binomial n k))``, 654 rpt strip_tac >- 655 metis_tac[ONE_LT_PRIME] >> 656 `(n = n-k + k) /\ (n-k) < n` by decide_tac >> 657 `FACT n = (binomial n k) * (FACT (n-k) * FACT k)` by metis_tac[binomial_formula] >> 658 `~(n divides (FACT k)) /\ ~(n divides (FACT (n-k)))` by metis_tac[PRIME_BIG_NOT_DIVIDES_FACT] >> 659 `n divides (FACT n)` by metis_tac[DIVIDES_FACT, LESS_TRANS] >> 660 metis_tac[P_EUCLIDES]); 661 662(* Theorem: n is prime ==> n divides C(n,k) for all 0 < k < n *) 663(* Proof: by prime_divides_binomials *) 664val prime_divides_binomials_alt = store_thm( 665 "prime_divides_binomials_alt", 666 ``!n k. prime n /\ 0 < k /\ k < n ==> n divides (binomial n k)``, 667 rw[prime_divides_binomials]); 668 669(* Theorem: If prime p divides n, p does not divide (n-1)!/(n-p)! *) 670(* Proof: 671 By contradiction. 672 (n-1)...(n-p+1)/p cannot be an integer 673 as p cannot divide any of the numerator. 674 Note: when p divides n, the nearest multiples for p are n+/-p. 675*) 676val prime_divisor_property = store_thm( 677 "prime_divisor_property", 678 ``!n p. 1 < n /\ p < n /\ prime p /\ p divides n ==> 679 ~(p divides ((FACT (n-1)) DIV (FACT (n-p))))``, 680 spose_not_then strip_assume_tac >> 681 `1 < p` by metis_tac[ONE_LT_PRIME] >> 682 `n-p < n-1` by decide_tac >> 683 `(FACT (n-1)) DIV (FACT (n-p)) = PROD_SET (IMAGE SUC ((count (n-1)) DIFF (count (n-p))))` 684 by metis_tac[FACT_REDUCTION, MULT_DIV, FACT_LESS] >> 685 `(count (n-1)) DIFF (count (n-p)) = {x | (n-p) <= x /\ x < (n-1)}` 686 by srw_tac[ARITH_ss][EXTENSION, EQ_IMP_THM] >> 687 `IMAGE SUC {x | (n-p) <= x /\ x < (n-1)} = {x | (n-p) < x /\ x < n}` by 688 (srw_tac[ARITH_ss][EXTENSION, EQ_IMP_THM] >> 689 qexists_tac `x-1` >> 690 decide_tac) >> 691 `FINITE (count (n - 1) DIFF count (n - p))` by rw[] >> 692 `?y. y IN {x| n - p < x /\ x < n} /\ p divides y` by metis_tac[PROD_SET_EUCLID, IMAGE_FINITE] >> 693 `!m n y. y IN {x | m < x /\ x < n} ==> m < y /\ y < n` by rw[] >> 694 `n-p < y /\ y < n` by metis_tac[] >> 695 `y < n + p` by decide_tac >> 696 `y = n` by metis_tac[MULTIPLE_INTERVAL] >> 697 decide_tac); 698 699(* Theorem: n divides C(n,k) for all 0 < k < n ==> n is prime *) 700(* Proof: 701 By contradiction. Let p be a proper factor of n, 1 < p < n. 702 Then C(n,p) = n(n-1)...(n-p+1)/p(p-1)..1 703 is divisible by n/p, but not n, since 704 C(n,p)/n = (n-1)...(n-p+1)/p(p-1)...1 705 cannot be an integer as p cannot divide any of the numerator. 706 Note: when p divides n, the nearest multiples for p are n+/-p. 707*) 708val divides_binomials_imp_prime = store_thm( 709 "divides_binomials_imp_prime", 710 ``!n. 1 < n /\ (!k. 0 < k /\ k < n ==> n divides (binomial n k)) ==> prime n``, 711 (spose_not_then strip_assume_tac) >> 712 `?p. prime p /\ p < n /\ p divides n` by metis_tac[PRIME_FACTOR_PROPER] >> 713 `n divides (binomial n p)` by metis_tac[PRIME_POS] >> 714 `0 < p` by metis_tac[PRIME_POS] >> 715 `(n = n-p + p) /\ (n-p) < n` by decide_tac >> 716 `FACT n = (binomial n p) * (FACT (n-p) * FACT p)` by metis_tac[binomial_formula] >> 717 `(n = SUC (n-1)) /\ (p = SUC (p-1))` by decide_tac >> 718 `(FACT n = n * FACT (n-1)) /\ (FACT p = p * FACT (p-1))` by metis_tac[FACT] >> 719 `n * FACT (n-1) = (binomial n p) * (FACT (n-p) * (p * FACT (p-1)))` by metis_tac[] >> 720 `0 < n` by decide_tac >> 721 `?q. binomial n p = n * q` by metis_tac[divides_def, MULT_COMM] >> 722 `0 <> n` by decide_tac >> 723 `FACT (n-1) = q * (FACT (n-p) * (p * FACT (p-1)))` 724 by metis_tac[EQ_MULT_LCANCEL, MULT_ASSOC] >> 725 `_ = q * ((FACT (p-1) * p)* FACT (n-p))` by metis_tac[MULT_COMM] >> 726 `_ = q * FACT (p-1) * p * FACT (n-p)` by metis_tac[MULT_ASSOC] >> 727 `FACT (n-1) DIV FACT (n-p) = q * FACT (p-1) * p` by metis_tac[MULT_DIV, FACT_LESS] >> 728 metis_tac[divides_def, prime_divisor_property]); 729 730(* Theorem: n is prime iff n divides C(n,k) for all 0 < k < n *) 731(* Proof: 732 By prime_divides_binomials and 733 divides_binomials_imp_prime. 734*) 735val prime_iff_divides_binomials = store_thm( 736 "prime_iff_divides_binomials", 737 ``!n. prime n <=> 1 < n /\ (!k. 0 < k /\ k < n ==> n divides (binomial n k))``, 738 metis_tac[prime_divides_binomials, divides_binomials_imp_prime]); 739 740(* Theorem: prime n <=> 1 < n /\ !k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0) *) 741(* Proof: by prime_iff_divides_binomials *) 742val prime_iff_divides_binomials_alt = store_thm( 743 "prime_iff_divides_binomials_alt", 744 ``!n. prime n <=> 1 < n /\ !k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0)``, 745 rw[prime_iff_divides_binomials, DIVIDES_MOD_0]); 746 747(* ------------------------------------------------------------------------- *) 748(* Binomial Theorem *) 749(* ------------------------------------------------------------------------- *) 750 751(* Theorem: Binomial Index Shifting, for 752 SUM (k=1..n) C(n,k)x^(n+1-k)y^k 753 = SUM (k=0..n-1) C(n,k+1)x^(n-k)y^(k+1) 754 *) 755(* Proof: 756SUM (k=1..n) C(n,k)x^(n+1-k)y^k 757= SUM (MAP (\k. (binomial n k)* x**(n+1-k) * y**k) (GENLIST SUC n)) 758= SUM (GENLIST (\k. (binomial n k)* x**(n+1-k) * y**k) o SUC n) 759 760SUM (k=0..n-1) C(n,k+1)x^(n-k)y^(k+1) 761= SUM (MAP (\k. (binomial n (k+1)) * x**(n-k) * y**(k+1)) (GENLIST I n)) 762= SUM (GENLIST (\k. (binomial n (k+1)) * x**(n-k) * y**(k+1)) o I n) 763= SUM (GENLIST (\k. (binomial n (k+1)) * x**(n-k) * y**(k+1)) n) 764 765i.e. 766 767(\k. (binomial n k)* x**(n-k+1) * y**k) o SUC 768= (\k. (binomial n (k+1)) * x**(n-k) * y**(k+1)) 769*) 770(* Theorem: Binomial index shift for GENLIST *) 771val GENLIST_binomial_index_shift = store_thm( 772 "GENLIST_binomial_index_shift", 773 ``!n x y. GENLIST ((\k. binomial n k * x ** SUC(n - k) * y ** k) o SUC) n = 774 GENLIST (\k. binomial n (SUC k) * x ** (n-k) * y**(SUC k)) n``, 775 rw_tac std_ss[GENLIST_FUN_EQ] >> 776 `SUC (n - SUC k) = n - k` by decide_tac >> 777 rw_tac std_ss[]); 778 779(* This is closely related to above, with (SUC n) replacing (n), 780 but does not require k < n. *) 781(* Proof: by function equality. *) 782val binomial_index_shift = store_thm( 783 "binomial_index_shift", 784 ``!n x y. (\k. binomial (SUC n) k * x ** ((SUC n) - k) * y ** k) o SUC = 785 (\k. binomial (SUC n) (SUC k) * x ** (n-k) * y ** (SUC k))``, 786 rw_tac std_ss[FUN_EQ_THM]); 787 788(* Pattern for binomial expansion: 789 790 (x+y)(x^3 + 3x^2y + 3xy^2 + y^3) 791 = x(x^3) + 3x(x^2y) + 3x(xy^2) + x(y^3) + 792 y(x^3) + 3y(x^2y) + 3y(xy^2) + y(y^3) 793 = x^4 + (3+1)x^3y + (3+3)(x^2y^2) + (1+3)(xy^3) + y^4 794 = x^4 + 4x^3y + 6x^2y^2 + 4xy^3 + y^4 795 796*) 797 798(* Theorem: multiply x into a binomial term *) 799(* Proof: by function equality and EXP. *) 800val binomial_term_merge_x = store_thm( 801 "binomial_term_merge_x", 802 ``!n x y. (\k. x * k) o (\k. binomial n k * x ** (n - k) * y ** k) = 803 (\k. binomial n k * x ** (SUC(n - k)) * y ** k)``, 804 rw_tac std_ss[FUN_EQ_THM] >> 805 `x * (binomial n k * x ** (n - k) * y ** k) = 806 binomial n k * (x * x ** (n - k)) * y ** k` by decide_tac >> 807 metis_tac[EXP]); 808 809(* Theorem: multiply y into a binomial term *) 810(* Proof: by functional equality and EXP. *) 811val binomial_term_merge_y = store_thm( 812 "binomial_term_merge_y", 813 ``!n x y. (\k. y * k) o (\k. binomial n k * x ** (n - k) * y ** k) = 814 (\k. binomial n k * x ** (n - k) * y ** (SUC k))``, 815 rw_tac std_ss[FUN_EQ_THM] >> 816 `y * (binomial n k * x ** (n - k) * y ** k) = 817 binomial n k * x ** (n - k) * (y * y ** k)` by decide_tac >> 818 metis_tac[EXP]); 819 820(* Theorem: [Binomial Theorem] (x + y)^n = SUM (k=0..n) C(n,k)x^(n-k)y^k *) 821(* Proof: 822 By induction on n. 823 Base case: to prove (x + y)^0 = SUM (k=0..0) C(0,k)x^(0-k)y^k 824 (x + y)^0 = 1 by EXP 825 SUM (k=0..0) C(0,k)x^(n-k)y^k = C(0,0)x^(0-0)y^0 = C(0,0) = 1 by EXP, binomial_def 826 Step case: assume (x + y)^n = SUM (k=0..n) C(n,k)x^(n-k)y^k 827 to prove: (x + y)^SUC n = SUM (k=0..(SUC n)) C(SUC n,k)x^((SUC n)-k)y^k 828 (x + y)^SUC n 829 = (x + y)(x + y)^n by EXP 830 = (x + y) SUM (k=0..n) C(n,k)x^(n-k)y^k by induction hypothesis 831 = x (SUM (k=0..n) C(n,k)x^(n-k)y^k) + 832 y (SUM (k=0..n) C(n,k)x^(n-k)y^k) by RIGHT_ADD_DISTRIB 833 = SUM (k=0..n) C(n,k)x^(n+1-k)y^k + 834 SUM (k=0..n) C(n,k)x^(n-k)y^(k+1) by moving factor into SUM 835 = C(n,0)x^(n+1) + SUM (k=1..n) C(n,k)x^(n+1-k)y^k + 836 SUM (k=0..n-1) C(n,k)x^(n-k)y^(k+1) + C(n,n)y^(n+1) 837 by breaking sum 838 839 = C(n,0)x^(n+1) + SUM (k=0..n-1) C(n,k+1)x^(n-k)y^(k+1) + 840 SUM (k=0..n-1) C(n,k)x^(n-k)y^(k+1) + C(n,n)y^(n+1) 841 by index shifting 842 = C(n,0)x^(n+1) + 843 SUM (k=0..n-1) [C(n,k+1) + C(n,k)] x^(n-k)y^(k+1) + 844 C(n,n)y^(n+1) by merging sums 845 = C(n,0)x^(n+1) + 846 SUM (k=0..n-1) C(n+1,k+1) x^(n-k)y^(k+1) + 847 C(n,n)y^(n+1) by binomial recurrence 848 = C(n,0)x^(n+1) + 849 SUM (k=1..n) C(n+1,k) x^(n+1-k)y^k + 850 C(n,n)y^(n+1) by index shifting again 851 = C(n+1,0)x^(n+1) + 852 SUM (k=1..n) C(n+1,k) x^(n+1-k)y^k + 853 C(n+1,n+1)y^(n+1) by binomial identities 854 = SUM (k=0..(SUC n))C(SUC n,k) x^((SUC n)-k)y^k 855 by synthesis of sum 856*) 857val binomial_thm = store_thm( 858 "binomial_thm", 859 ``!n x y. (x + y) ** n = SUM (GENLIST (\k. (binomial n k) * x ** (n-k) * y ** k) (SUC n))``, 860 Induct_on `n` >- 861 rw[EXP, binomial_n_n] >> 862 rw_tac std_ss[EXP] >> 863 `(x + y) * SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** k) (SUC n)) = 864 x * SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** k) (SUC n)) + 865 y * SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** k) (SUC n))` 866 by metis_tac[RIGHT_ADD_DISTRIB] >> 867 `_ = SUM (GENLIST ((\k. x * k) o (\k. binomial n k * x ** (n - k) * y ** k)) (SUC n)) + 868 SUM (GENLIST ((\k. y * k) o (\k. binomial n k * x ** (n - k) * y ** k)) (SUC n))` 869 by metis_tac[SUM_MULT, MAP_GENLIST] >> 870 `_ = SUM (GENLIST (\k. binomial n k * x ** SUC(n - k) * y ** k) (SUC n)) + 871 SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** (SUC k)) (SUC n))` 872 by rw[binomial_term_merge_x, binomial_term_merge_y] >> 873 `_ = (\k. binomial n k * x ** SUC (n - k) * y ** k) 0 + 874 SUM (GENLIST ((\k. binomial n k * x ** SUC (n - k) * y ** k) o SUC) n) + 875 SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** (SUC k)) (SUC n))` 876 by rw[SUM_DECOMPOSE_FIRST] >> 877 `_ = (\k. binomial n k * x ** SUC (n - k) * y ** k) 0 + 878 SUM (GENLIST ((\k. binomial n k * x ** SUC (n - k) * y ** k) o SUC) n) + 879 (SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n) + 880 (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n )` 881 by rw[SUM_DECOMPOSE_LAST] >> 882 `_ = (\k. binomial n k * x ** SUC(n - k) * y ** k) 0 + 883 SUM (GENLIST (\k. binomial n (SUC k) * x ** (n - k) * y ** (SUC k)) n) + 884 (SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n) + 885 (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n )` 886 by metis_tac[GENLIST_binomial_index_shift] >> 887 `_ = (\k. binomial n k * x ** SUC(n - k) * y ** k) 0 + 888 (SUM (GENLIST (\k. binomial n (SUC k) * x ** (n - k) * y ** (SUC k)) n) + 889 SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n)) + 890 (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n` 891 by decide_tac >> 892 `_ = (\k. binomial n k * x ** SUC (n - k) * y ** k) 0 + 893 SUM (GENLIST (\k. (binomial n (SUC k) * x ** (n - k) * y ** (SUC k) + 894 binomial n k * x ** (n - k) * y ** (SUC k))) n) + 895 (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n` 896 by metis_tac[SUM_ADD_GENLIST] >> 897 `_ = (\k. binomial n k * x ** SUC(n - k) * y ** k) 0 + 898 SUM (GENLIST (\k. (binomial n (SUC k) + binomial n k) * x ** (n - k) * y ** (SUC k)) n) + 899 (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n` 900 by rw[RIGHT_ADD_DISTRIB, MULT_ASSOC] >> 901 `_ = (\k. binomial n k * x ** SUC(n - k) * y ** k) 0 + 902 SUM (GENLIST (\k. binomial (SUC n) (SUC k) * x ** (n - k) * y ** (SUC k)) n) + 903 (\k. binomial n k * x ** (n - k) * y ** (SUC k)) n` 904 by rw[binomial_recurrence, ADD_COMM] >> 905 `_ = binomial (SUC n) 0 * x ** (SUC n) * y ** 0 + 906 SUM (GENLIST (\k. binomial (SUC n) (SUC k) * x ** (n - k) * y ** (SUC k)) n) + 907 binomial (SUC n) (SUC n) * x ** 0 * y ** (SUC n)` 908 by rw[binomial_n_0, binomial_n_n] >> 909 `_ = binomial (SUC n) 0 * x ** (SUC n) * y ** 0 + 910 SUM (GENLIST ((\k. binomial (SUC n) k * x ** ((SUC n) - k) * y ** k) o SUC) n) + 911 binomial (SUC n) (SUC n) * x ** 0 * y ** (SUC n)` 912 by rw[binomial_index_shift] >> 913 `_ = SUM (GENLIST (\k. binomial (SUC n) k * x ** (SUC n - k) * y ** k) (SUC n)) + 914 (\k. binomial (SUC n) k * x ** (SUC n - k) * y ** k) (SUC n)` 915 by rw[SUM_DECOMPOSE_FIRST] >> 916 `_ = SUM (GENLIST (\k. binomial (SUC n) k * x ** (SUC n - k) * y ** k) (SUC (SUC n)))` 917 by rw[SUM_DECOMPOSE_LAST] >> 918 decide_tac); 919 920(* This is a milestone theorem. *) 921 922(* Derive an alternative form. *) 923val binomial_thm_alt = save_thm("binomial_thm_alt", 924 binomial_thm |> SIMP_RULE bool_ss [ADD1]); 925(* val binomial_thm_alt = 926 |- !n x y. (x + y) ** n = 927 SUM (GENLIST (\k. binomial n k * x ** (n - k) * y ** k) (n + 1)): thm *) 928 929(* Theorem: SUM (GENLIST (binomial n) (SUC n)) = 2 ** n *) 930(* Proof: by binomial_sum_alt and function equality. *) 931(* Proof: 932 Put x = 1, y = 1 in binomial_thm, 933 (1 + 1) ** n = SUM (GENLIST (\k. binomial n k * 1 ** (n - k) * 1 ** k) (SUC n)) 934 (1 + 1) ** n = SUM (GENLIST (\k. binomial n k) (SUC n)) by EXP_1 935 or 2 ** n = SUM (GENLIST (binomial n) (SUC n)) by FUN_EQ_THM 936*) 937Theorem binomial_sum: 938 !n. SUM (GENLIST (binomial n) (SUC n)) = 2 ** n 939Proof 940 rpt strip_tac >> 941 `!n. (\k. binomial n k * 1 ** (n - k) * 1 ** k) = binomial n` by rw[FUN_EQ_THM] >> 942 `SUM (GENLIST (binomial n) (SUC n)) = 943 SUM (GENLIST (\k. binomial n k * 1 ** (n - k) * 1 ** k) (SUC n))` by fs[] >> 944 `_ = (1 + 1) ** n` by rw[GSYM binomial_thm] >> 945 simp[] 946QED 947 948(* Derive an alternative form. *) 949val binomial_sum_alt = save_thm("binomial_sum_alt", 950 binomial_sum |> SIMP_RULE bool_ss [ADD1]); 951(* val binomial_sum_alt = |- !n. SUM (GENLIST (binomial n) (n + 1)) = 2 ** n: thm *) 952 953(* ------------------------------------------------------------------------- *) 954(* Binomial Horizontal List *) 955(* ------------------------------------------------------------------------- *) 956 957(* Define Horizontal List in Pascal Triangle *) 958(* 959val binomial_horizontal_def = Define ` 960 binomial_horizontal n = GENLIST (binomial n) (SUC n) 961`; 962*) 963 964(* Use overloading for binomial_horizontal n. *) 965val _ = overload_on("binomial_horizontal", ``\n. GENLIST (binomial n) (n + 1)``); 966 967(* Theorem: binomial_horizontal 0 = [1] *) 968(* Proof: 969 binomial_horizontal 0 970 = GENLIST (binomial 0) (0 + 1) by notation 971 = SNOC (binomial 0 0) [] by GENLIST, ONE 972 = [binomial 0 0] by SNOC 973 = [1] by binomial_n_0 974*) 975val binomial_horizontal_0 = store_thm( 976 "binomial_horizontal_0", 977 ``binomial_horizontal 0 = [1]``, 978 rw[binomial_n_0]); 979 980(* Theorem: LENGTH (binomial_horizontal n) = n + 1 *) 981(* Proof: 982 LENGTH (binomial_horizontal n) 983 = LENGTH (GENLIST (binomial n) (n + 1)) by notation 984 = n + 1 by LENGTH_GENLIST 985*) 986val binomial_horizontal_len = store_thm( 987 "binomial_horizontal_len", 988 ``!n. LENGTH (binomial_horizontal n) = n + 1``, 989 rw[]); 990 991(* Theorem: k < n + 1 ==> MEM (binomial n k) (binomial_horizontal n) *) 992(* Proof: by MEM_GENLIST *) 993val binomial_horizontal_mem = store_thm( 994 "binomial_horizontal_mem", 995 ``!n k. k < n + 1 ==> MEM (binomial n k) (binomial_horizontal n)``, 996 metis_tac[MEM_GENLIST]); 997 998(* Theorem: MEM (binomial n k) (binomial_horizontal n) <=> k <= n *) 999(* Proof: 1000 If part: MEM (binomial n k) (binomial_horizontal n) ==> k <= n 1001 By contradiction, suppose n < k. 1002 Then binomial n k = 0 by binomial_less_0, ~(k <= n) 1003 But ?m. m < n + 1 ==> 0 = binomial n m by MEM_GENLIST 1004 or m <= n ==> binomial n m = 0 by m < n + 1 1005 Yet binomial n m <> 0 by binomial_eq_0 1006 This is a contradiction. 1007 Only-if part: k <= n ==> MEM (binomial n k) (binomial_horizontal n) 1008 By MEM_GENLIST, this is to show: 1009 ?m. m < n + 1 /\ (binomial n k = binomial n m) 1010 Note k <= n ==> k < n + 1, 1011 Take m = k, the result follows. 1012*) 1013val binomial_horizontal_mem_iff = store_thm( 1014 "binomial_horizontal_mem_iff", 1015 ``!n k. MEM (binomial n k) (binomial_horizontal n) <=> k <= n``, 1016 rw[EQ_IMP_THM] >| [ 1017 spose_not_then strip_assume_tac >> 1018 `binomial n k = 0` by rw[binomial_less_0] >> 1019 fs[MEM_GENLIST] >> 1020 `m <= n` by decide_tac >> 1021 fs[binomial_eq_0], 1022 rw[MEM_GENLIST] >> 1023 `k < n + 1` by decide_tac >> 1024 metis_tac[] 1025 ]); 1026 1027(* Theorem: MEM x (binomial_horizontal n) <=> ?k. k <= n /\ (x = binomial n k) *) 1028(* Proof: 1029 By MEM_GENLIST, this is to show: 1030 (?m. m < n + 1 /\ (x = binomial n m)) <=> ?k. k <= n /\ (x = binomial n k) 1031 Since m < n + 1 <=> m <= n by LE_LT1 1032 This is trivially true. 1033*) 1034val binomial_horizontal_member = store_thm( 1035 "binomial_horizontal_member", 1036 ``!n x. MEM x (binomial_horizontal n) <=> ?k. k <= n /\ (x = binomial n k)``, 1037 metis_tac[MEM_GENLIST, LE_LT1]); 1038 1039(* Theorem: k <= n ==> (EL k (binomial_horizontal n) = binomial n k) *) 1040(* Proof: by EL_GENLIST *) 1041val binomial_horizontal_element = store_thm( 1042 "binomial_horizontal_element", 1043 ``!n k. k <= n ==> (EL k (binomial_horizontal n) = binomial n k)``, 1044 rw[EL_GENLIST]); 1045 1046(* Theorem: EVERY (\x. 0 < x) (binomial_horizontal n) *) 1047(* Proof: 1048 EVERY (\x. 0 < x) (binomial_horizontal n) 1049 <=> EVERY (\x. 0 < x) (GENLIST (binomial n) (n + 1)) by notation 1050 <=> !k. k < n + 1 ==> 0 < binomial n k by EVERY_GENLIST 1051 <=> !k. k <= n ==> 0 < binomial n k by arithmetic 1052 <=> T by binomial_pos 1053*) 1054val binomial_horizontal_pos = store_thm( 1055 "binomial_horizontal_pos", 1056 ``!n. EVERY (\x. 0 < x) (binomial_horizontal n)``, 1057 rpt strip_tac >> 1058 `!k n. k < n + 1 <=> k <= n` by decide_tac >> 1059 rw_tac std_ss[EVERY_GENLIST, LESS_EQ_IFF_LESS_SUC, binomial_pos]); 1060 1061(* Theorem: MEM x (binomial_horizontal n) ==> 0 < x *) 1062(* Proof: by binomial_horizontal_pos, EVERY_MEM *) 1063val binomial_horizontal_pos_alt = store_thm( 1064 "binomial_horizontal_pos_alt", 1065 ``!n x. MEM x (binomial_horizontal n) ==> 0 < x``, 1066 metis_tac[binomial_horizontal_pos, EVERY_MEM]); 1067 1068(* Theorem: SUM (binomial_horizontal n) = 2 ** n *) 1069(* Proof: 1070 SUM (binomial_horizontal n) 1071 = SUM (GENLIST (binomial n) (n + 1)) by notation 1072 = 2 ** n by binomial_sum, ADD1 1073*) 1074val binomial_horizontal_sum = store_thm( 1075 "binomial_horizontal_sum", 1076 ``!n. SUM (binomial_horizontal n) = 2 ** n``, 1077 rw_tac std_ss[binomial_sum, GSYM ADD1]); 1078 1079(* Theorem: MAX_LIST (binomial_horizontal n) = binomial n (HALF n) *) 1080(* Proof: 1081 Let l = binomial_horizontal n, m = binomial n (HALF n). 1082 Then l <> [] by binomial_horizontal_len, LENGTH_NIL 1083 and HALF n <= n by DIV_LESS_EQ, 0 < 2 1084 or HALF n < n + 1 by arithmetic 1085 Also MEM m l by binomial_horizontal_mem 1086 and !x. MEM x l ==> x <= m by binomial_max, MEM_GENLIST 1087 Thus m = MAX_LIST l by MAX_LIST_TEST 1088*) 1089val binomial_horizontal_max = store_thm( 1090 "binomial_horizontal_max", 1091 ``!n. MAX_LIST (binomial_horizontal n) = binomial n (HALF n)``, 1092 rpt strip_tac >> 1093 qabbrev_tac `l = binomial_horizontal n` >> 1094 qabbrev_tac `m = binomial n (HALF n)` >> 1095 `l <> []` by metis_tac[binomial_horizontal_len, LENGTH_NIL, DECIDE``n + 1 <> 0``] >> 1096 `HALF n <= n` by rw[DIV_LESS_EQ] >> 1097 `HALF n < n + 1` by decide_tac >> 1098 `MEM m l` by rw[binomial_horizontal_mem, Abbr`l`, Abbr`m`] >> 1099 metis_tac[binomial_max, MEM_GENLIST, MAX_LIST_TEST]); 1100 1101(* Theorem: MAX_SET (IMAGE (binomial n) (count (n + 1))) = binomial n (HALF n) *) 1102(* Proof: 1103 Let f = binomial n, s = IMAGE f (count (n + 1)). 1104 Note FINITE (count (n + 1)) by FINITE_COUNT 1105 so FINITE s by IMAGE_FINITE 1106 Also count (n + 1) <> {} by COUNT_EQ_EMPTY, n + 1 <> 0 1107 so s <> {} by IMAGE_EQ_EMPTY 1108 Now !k. k IN (count (n + 1)) ==> f k <= f (HALF n) by binomial_max 1109 ==> !x. x IN s ==> x <= f (HALF n) by IN_IMAGE 1110 Also HALF n <= n by DIV_LESS_EQ, 0 < 2 1111 so HALF n IN (count (n + 1)) by IN_COUNT 1112 ==> f (HALF n) IN s by IN_IMAGE 1113 Thus MAX_SET s = f (HALF n) by MAX_SET_TEST 1114*) 1115val binomial_row_max = store_thm( 1116 "binomial_row_max", 1117 ``!n. MAX_SET (IMAGE (binomial n) (count (n + 1))) = binomial n (HALF n)``, 1118 rpt strip_tac >> 1119 qabbrev_tac `f = binomial n` >> 1120 qabbrev_tac `s = IMAGE f (count (n + 1))` >> 1121 `FINITE s` by rw[Abbr`s`] >> 1122 `s <> {}` by rw[COUNT_EQ_EMPTY, Abbr`s`] >> 1123 `!k. k IN (count (n + 1)) ==> f k <= f (HALF n)` by rw[binomial_max, Abbr`f`] >> 1124 `!x. x IN s ==> x <= f (HALF n)` by metis_tac[IN_IMAGE] >> 1125 `HALF n <= n` by rw[DIV_LESS_EQ] >> 1126 `HALF n IN (count (n + 1))` by rw[] >> 1127 `f (HALF n) IN s` by metis_tac[IN_IMAGE] >> 1128 rw[MAX_SET_TEST]); 1129 1130(* Theorem: k <= m /\ m <= n ==> 1131 ((binomial m k) * (binomial n m) = (binomial n k) * (binomial (n - k) (m - k))) *) 1132(* Proof: 1133 Using binomial_formula2, 1134 1135 (binomial m k) * (binomial n m) 1136 n! m! 1137 = ----------- * ------------------ binomial formula 1138 m! (n - m)! k! (m - k)! 1139 n! m! 1140 = ----------- * ------------------ cancel m! 1141 k! m! (m - k)! (n - m)! 1142 n! (n - k)! 1143 = ----------- * ------------------ replace by (n - k)! 1144 k! (n - k)! (m - k)! (n - m)! 1145 1146 = (binomial n k) * (binomial (n - k) (m - k)) binomial formula 1147*) 1148val binomial_product_identity = store_thm( 1149 "binomial_product_identity", 1150 ``!m n k. k <= m /\ m <= n ==> 1151 ((binomial m k) * (binomial n m) = (binomial n k) * (binomial (n - k) (m - k)))``, 1152 rpt strip_tac >> 1153 `m - k <= n - k` by decide_tac >> 1154 `(n - k) - (m - k) = n - m` by decide_tac >> 1155 `FACT m = binomial m k * (FACT (m - k) * FACT k)` by rw[binomial_formula2] >> 1156 `FACT n = binomial n m * (FACT (n - m) * FACT m)` by rw[binomial_formula2] >> 1157 `FACT n = binomial n k * (FACT (n - k) * FACT k)` by rw[binomial_formula2] >> 1158 `FACT (n - k) = binomial (n - k) (m - k) * (FACT (n - m) * FACT (m - k))` by metis_tac[binomial_formula2] >> 1159 `FACT n = FACT (n - m) * (FACT k * (FACT (m - k) * ((binomial m k) * (binomial n m))))` by metis_tac[MULT_ASSOC, MULT_COMM] >> 1160 `FACT n = FACT (n - m) * (FACT k * (FACT (m - k) * ((binomial n k) * (binomial (n - k) (m - k)))))` by metis_tac[MULT_ASSOC, MULT_COMM] >> 1161 metis_tac[MULT_LEFT_CANCEL, FACT_LESS, NOT_ZERO_LT_ZERO]); 1162 1163(* Theorem: binomial n (HALF n) <= 4 ** (HALF n) *) 1164(* Proof: 1165 Let m = HALF n, l = binomial_horizontal n 1166 Note LENGTH l = n + 1 by binomial_horizontal_len 1167 If EVEN n, 1168 Then n = 2 * m by EVEN_HALF 1169 and m <= n by m <= 2 * m 1170 Note EL m l <= SUM l by SUM_LE_EL, m < n + 1 1171 Now EL m l = binomial n m by binomial_horizontal_element, m <= n 1172 and SUM l 1173 = 2 ** n by binomial_horizontal_sum 1174 = 4 ** m by EXP_EXP_MULT 1175 Hence binomial n m <= 4 ** m. 1176 If ~EVEN n, 1177 Then ODD n by EVEN_ODD 1178 and n = 2 * m + 1 by ODD_HALF 1179 so m + 1 <= n by m + 1 <= 2 * m + 1 1180 with m <= n by m + 1 <= n 1181 Note EL m l = binomial n m by binomial_horizontal_element, m <= n 1182 and EL (m + 1) l = binomial n (m + 1) by binomial_horizontal_element, m + 1 <= n 1183 Note binomial n (m + 1) = binomial n m by binomial_sym 1184 Thus 2 * binomial n m 1185 = binomial n m + binomial n (m + 1) by above 1186 = EL m l + EL (m + 1) l 1187 <= SUM l by SUM_LE_SUM_EL, m < m + 1, m + 1 < n + 1 1188 and SUM l 1189 = 2 ** n by binomial_horizontal_sum 1190 = 2 * 2 ** (2 * m) by EXP, ADD1 1191 = 2 * 4 ** m by EXP_EXP_MULT 1192 Hence binomial n m <= 4 ** m. 1193*) 1194val binomial_middle_upper_bound = store_thm( 1195 "binomial_middle_upper_bound", 1196 ``!n. binomial n (HALF n) <= 4 ** (HALF n)``, 1197 rpt strip_tac >> 1198 qabbrev_tac `m = HALF n` >> 1199 qabbrev_tac `l = binomial_horizontal n` >> 1200 `LENGTH l = n + 1` by rw[binomial_horizontal_len, Abbr`l`] >> 1201 Cases_on `EVEN n` >| [ 1202 `n = 2 * m` by rw[EVEN_HALF, Abbr`m`] >> 1203 `m < n + 1` by decide_tac >> 1204 `EL m l <= SUM l` by rw[SUM_LE_EL] >> 1205 `EL m l = binomial n m` by rw[binomial_horizontal_element, Abbr`l`] >> 1206 `SUM l = 2 ** n` by rw[binomial_horizontal_sum, Abbr`l`] >> 1207 `_ = 4 ** m` by rw[EXP_EXP_MULT] >> 1208 decide_tac, 1209 `ODD n` by metis_tac[EVEN_ODD] >> 1210 `n = 2 * m + 1` by rw[ODD_HALF, Abbr`m`] >> 1211 `EL m l = binomial n m` by rw[binomial_horizontal_element, Abbr`l`] >> 1212 `EL (m + 1) l = binomial n (m + 1)` by rw[binomial_horizontal_element, Abbr`l`] >> 1213 `binomial n (m + 1) = binomial n m` by rw[Once binomial_sym] >> 1214 `EL m l + EL (m + 1) l <= SUM l` by rw[SUM_LE_SUM_EL] >> 1215 `SUM l = 2 ** n` by rw[binomial_horizontal_sum, Abbr`l`] >> 1216 `_ = 2 * 2 ** (2 * m)` by metis_tac[EXP, ADD1] >> 1217 `_ = 2 * 4 ** m` by rw[EXP_EXP_MULT] >> 1218 decide_tac 1219 ]); 1220 1221(* ------------------------------------------------------------------------- *) 1222(* Stirling's Approximation *) 1223(* ------------------------------------------------------------------------- *) 1224 1225(* Stirling's formula: n! ~ sqrt(2 pi n) (n/e)^n. *) 1226val _ = overload_on("Stirling", 1227 ``(!n. FACT n = (SQRT (2 * pi * n)) * (n DIV e) ** n) /\ 1228 (!n. SQRT n = n ** h) /\ (2 * h = 1) /\ (0 < pi) /\ (0 < e) /\ 1229 (!a b x y. (a * b) DIV (x * y) = (a DIV x) * (b DIV y)) /\ 1230 (!a b c. (a DIV c) DIV (b DIV c) = a DIV b)``); 1231 1232(* Theorem: Stirling ==> 1233 !n. 0 < n /\ EVEN n ==> (binomial n (HALF n) = (2 ** (n + 1)) DIV (SQRT (2 * pi * n))) *) 1234(* Proof: 1235 Note HALF n <= n by DIV_LESS_EQ, 0 < 2 1236 Let k = HALF n, then n = 2 * k by EVEN_HALF 1237 Note 0 < k by 0 < n = 2 * k 1238 so (k * 2) DIV k = 2 by MULT_TO_DIV, 0 < k 1239 or n DIV k = 2 by MULT_COMM 1240 Also 0 < pi * n by MULT_EQ_0, 0 < pi, 0 < n 1241 so 0 < 2 * pi * n by arithmetic 1242 1243 Some theorems on the fly: 1244 Claim: !a b j. (a ** j) DIV (b ** j) = (a DIV b) ** j [1] 1245 Proof: By induction on j. 1246 Base: (a ** 0) DIV (b ** 0) = (a DIV b) ** 0 1247 (a ** 0) DIV (b ** 0) 1248 = 1 DIV 1 = 1 by EXP, DIVMOD_ID, 0 < 1 1249 = (a DIV b) ** 0 by EXP 1250 Step: (a ** j) DIV (b ** j) = (a DIV b) ** j ==> 1251 (a ** (SUC j)) DIV (b ** (SUC j)) = (a DIV b) ** (SUC j) 1252 (a ** (SUC j)) DIV (b ** (SUC j)) 1253 = (a * a ** j) DIV (b * b ** j) by EXP 1254 = (a DIV b) * ((a ** j) DIV (b ** j)) by assumption 1255 = (a DIV b) * (a DIV b) ** j by induction hypothesis 1256 = (a DIV b) ** (SUC j) by EXP 1257 1258 Claim: !a b c. (a DIV b) * c = (a * c) DIV b [2] 1259 Proof: (a DIV b) * c 1260 = (a DIV b) * (c DIV 1) by DIV_1 1261 = (a * c) DIV (b * 1) by assumption 1262 = (a * c) DIV b by MULT_RIGHT_1 1263 1264 Claim: !a b. a DIV b = 2 * (a DIV (2 * b)) [3] 1265 Proof: a DIV b 1266 = 1 * (a DIV b) by MULT_LEFT_1 1267 = (n DIV n) * (a DIV b) by DIVMOD_ID, 0 < n 1268 = (n * a) DIV (n * b) by assumption 1269 = (n * a) DIV (k * (2 * b)) by arithmetic, n = 2 * k 1270 = (n DIV k) * (a DIV (2 * b)) by assumption 1271 = 2 * (a DIV (2 * b)) by n DIV k = 2 1272 1273 Claim: !a b. 0 < b ==> (a * (b ** h DIV b) = a DIV (b ** h)) [4] 1274 Proof: Let c = b ** h. 1275 Then b = c * c by EXP_EXP_MULT 1276 so 0 < c by MULT_EQ_0, 0 < b 1277 a * (c DIV b) 1278 = (c DIV b) * a by MULT_COMM 1279 = (a * c) DIV b by [2] 1280 = (a * c) DIV (c * c) by b = c * c 1281 = (a DIV c) * (c DIV c) by assumption 1282 = a DIV c by DIVMOD_ID, c DIV c = 1, 0 < c 1283 1284 Note (FACT k) ** 2 1285 = (SQRT (2 * pi * k)) ** 2 * ((k DIV e) ** k) ** 2 by EXP_BASE_MULT 1286 = (SQRT (2 * pi * k)) ** 2 * (k DIV e) ** n by EXP_EXP_MULT, n = 2 * k 1287 = (SQRT (pi * n)) ** 2 * (k DIV e) ** n by MULT_ASSOC, 2 * k = n 1288 = ((pi * n) ** h) ** 2 * (k DIV e) ** n by assumption 1289 = (pi * n) * (k DIV e) ** n by EXP_EXP_MULT, h * 2 = 1 1290 1291 binomial n (HALF n) 1292 = binomial n k by k = HALF n 1293 = FACT n DIV (FACT k * FACT (n - k)) by binomial_formula3, k <= n 1294 = FACT n DIV (FACT k * FACT k) by arithmetic, n - k = 2 * k - k = k 1295 = FACT n DIV ((FACT k) ** 2) by EXP_2 1296 = FACT n DIV ((pi * n) * (k DIV e) ** n) by above 1297 = ((2 * pi * n) ** h * (n DIV e) ** n) DIV ((pi * n) * (k DIV e) ** n) by assumption 1298 = ((2 * pi * n) ** h DIV (pi * n)) * ((n DIV e) ** n DIV ((k DIV e) ** n)) by (a * b) DIV (x * y) = (a DIV x) * (b DIV y) 1299 = ((2 * pi * n) ** h DIV (pi * n)) * ((n DIV e) DIV (k DIV e)) ** n by (a ** n) DIV (b ** n) = (a DIV b) ** n) 1300 = 2 * ((2 * pi * n) ** h DIV (2 * pi * n)) * ((n DIV e) DIV (k DIV e)) ** n by MULT_ASSOC, a DIV b = 2 * a DIV (2 * b) 1301 = 2 * ((2 * pi * n) ** h DIV (2 * pi * n)) * (n DIV k) ** n by assumption, apply DIV_DIV_DIV_MULT 1302 = 2 DIV (2 * pi * n) ** h * (n DIV k) ** n by 2 * x ** h DIV x = 2 DIV (x ** h) 1303 = 2 DIV (2 * pi * n) ** h * 2 ** n by n DIV k = 2 1304 = 2 * 2 ** n DIV (2 * pi * n) ** h by (a DIV b) * c = a * c DIV b 1305 = 2 ** (SUC n) DIV (2 * pi * n) ** h by EXP 1306 = 2 ** (n + 1)) DIV (SQRT (2 * pi * n)) by ADD1, assumption 1307*) 1308val binomial_middle_by_stirling = store_thm( 1309 "binomial_middle_by_stirling", 1310 ``Stirling ==> !n. 0 < n /\ EVEN n ==> (binomial n (HALF n) = (2 ** (n + 1)) DIV (SQRT (2 * pi * n)))``, 1311 rpt strip_tac >> 1312 `HALF n <= n /\ (n = 2 * HALF n)` by rw[DIV_LESS_EQ, EVEN_HALF] >> 1313 qabbrev_tac `k = HALF n` >> 1314 `0 < k` by decide_tac >> 1315 `n DIV k = 2` by metis_tac[MULT_TO_DIV, MULT_COMM] >> 1316 `0 < pi * n` by metis_tac[MULT_EQ_0, NOT_ZERO_LT_ZERO] >> 1317 `0 < 2 * pi * n` by decide_tac >> 1318 `(FACT k) ** 2 = (SQRT (2 * pi * k)) ** 2 * ((k DIV e) ** k) ** 2` by rw[EXP_BASE_MULT] >> 1319 `_ = (SQRT (2 * pi * k)) ** 2 * (k DIV e) ** n` by rw[GSYM EXP_EXP_MULT] >> 1320 `_ = (pi * n) * (k DIV e) ** n` by rw[GSYM EXP_EXP_MULT] >> 1321 (`!a b j. (a ** j) DIV (b ** j) = (a DIV b) ** j` by (Induct_on `j` >> rw[EXP])) >> 1322 `!a b c. (a DIV b) * c = (a * c) DIV b` by metis_tac[DIV_1, MULT_RIGHT_1] >> 1323 `!a b. a DIV b = 2 * (a DIV (2 * b))` by metis_tac[DIVMOD_ID, MULT_LEFT_1] >> 1324 `!a b. 0 < b ==> (a * (b ** h DIV b) = a DIV (b ** h))` by 1325 (rpt strip_tac >> 1326 qabbrev_tac `c = b ** h` >> 1327 `b = c * c` by rw[GSYM EXP_EXP_MULT, Abbr`c`] >> 1328 `0 < c` by metis_tac[MULT_EQ_0, NOT_ZERO_LT_ZERO] >> 1329 `a * (c DIV b) = (a * c) DIV (c * c)` by metis_tac[MULT_COMM] >> 1330 `_ = (a DIV c) * (c DIV c)` by metis_tac[] >> 1331 metis_tac[DIVMOD_ID, MULT_RIGHT_1]) >> 1332 `binomial n k = (FACT n) DIV (FACT k * FACT (n - k))` by metis_tac[binomial_formula3] >> 1333 `_ = (FACT n) DIV (FACT k) ** 2` by metis_tac[EXP_2, DECIDE``2 * k - k = k``] >> 1334 `_ = ((2 * pi * n) ** h * (n DIV e) ** n) DIV ((pi * n) * (k DIV e) ** n)` by prove_tac[] >> 1335 `_ = ((2 * pi * n) ** h DIV (pi * n)) * ((n DIV e) ** n DIV ((k DIV e) ** n))` by metis_tac[] >> 1336 `_ = ((2 * pi * n) ** h DIV (pi * n)) * ((n DIV e) DIV (k DIV e)) ** n` by metis_tac[] >> 1337 `_ = 2 * ((2 * pi * n) ** h DIV (2 * pi * n)) * ((n DIV e) DIV (k DIV e)) ** n` by metis_tac[MULT_ASSOC] >> 1338 `_ = 2 * ((2 * pi * n) ** h DIV (2 * pi * n)) * (n DIV k) ** n` by metis_tac[] >> 1339 `_ = 2 DIV (2 * pi * n) ** h * (n DIV k) ** n` by metis_tac[] >> 1340 `_ = 2 DIV (2 * pi * n) ** h * 2 ** n` by metis_tac[] >> 1341 `_ = (2 * 2 ** n DIV (2 * pi * n) ** h)` by metis_tac[] >> 1342 metis_tac[EXP, ADD1]); 1343 1344(* ------------------------------------------------------------------------- *) 1345(* Useful theorems for Binomial *) 1346(* ------------------------------------------------------------------------- *) 1347 1348(* Theorem: !k. 0 < k /\ k < n ==> (binomial n k MOD n = 0) <=> 1349 !h. 0 <= h /\ h < PRE n ==> (binomial n (SUC h) MOD n = 0) *) 1350(* Proof: by h = PRE k, or k = SUC h. 1351 If part: put k = SUC h, 1352 then 0 < SUC h ==> 0 <= h, 1353 and SUC h < n ==> PRE (SUC h) = h < PRE n by prim_recTheory.PRE 1354 Only-if part: put h = PRE k, 1355 then 0 <= PRE k ==> 0 < k 1356 and PRE k < PRE n ==> k < n by INV_PRE_LESS 1357*) 1358val binomial_range_shift = store_thm( 1359 "binomial_range_shift", 1360 ``!n . 0 < n ==> ((!k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0)) <=> 1361 (!h. h < PRE n ==> ((binomial n (SUC h)) MOD n = 0)))``, 1362 rw_tac std_ss[EQ_IMP_THM] >| [ 1363 `0 < SUC h /\ SUC h < n` by decide_tac >> 1364 rw_tac std_ss[], 1365 `k <> 0` by decide_tac >> 1366 `?h. k = SUC h` by metis_tac[num_CASES] >> 1367 `h < PRE n` by decide_tac >> 1368 rw_tac std_ss[] 1369 ]); 1370 1371(* Theorem: binomial n k MOD n = 0 <=> (binomial n k * x ** (n-k) * y ** k) MOD n = 0 *) 1372(* Proof: 1373 (binomial n k * x ** (n-k) * y ** k) MOD n = 0 1374 <=> (binomial n k * (x ** (n-k) * y ** k)) MOD n = 0 by MULT_ASSOC 1375 <=> (((binomial n k) MOD n) * ((x ** (n - k) * y ** k) MOD n)) MOD n = 0 by MOD_TIMES2 1376 If part, apply 0 * z = 0 by MULT. 1377 Only-if part, pick x = 1, y = 1, apply EXP_1. 1378*) 1379val binomial_mod_zero = store_thm( 1380 "binomial_mod_zero", 1381 ``!n. 0 < n ==> !k. (binomial n k MOD n = 0) <=> (!x y. (binomial n k * x ** (n-k) * y ** k) MOD n = 0)``, 1382 rw_tac std_ss[EQ_IMP_THM] >- 1383 metis_tac[MOD_TIMES2, ZERO_MOD, MULT] >> 1384 metis_tac[EXP_1, MULT_RIGHT_1]); 1385 1386 1387(* Theorem: (!k. 0 < k /\ k < n ==> (!x y. ((binomial n k * x ** (n - k) * y ** k) MOD n = 0))) <=> 1388 (!h. h < PRE n ==> (!x y. ((binomial n (SUC h) * x ** (n - (SUC h)) * y ** (SUC h)) MOD n = 0))) *) 1389(* Proof: by h = PRE k, or k = SUC h. *) 1390val binomial_range_shift_alt = store_thm( 1391 "binomial_range_shift_alt", 1392 ``!n . 0 < n ==> ((!k. 0 < k /\ k < n ==> (!x y. ((binomial n k * x ** (n - k) * y ** k) MOD n = 0))) <=> 1393 (!h. h < PRE n ==> (!x y. ((binomial n (SUC h) * x ** (n - (SUC h)) * y ** (SUC h)) MOD n = 0))))``, 1394 rw_tac std_ss[EQ_IMP_THM] >| [ 1395 `0 < SUC h /\ SUC h < n` by decide_tac >> 1396 rw_tac std_ss[], 1397 `k <> 0` by decide_tac >> 1398 `?h. k = SUC h` by metis_tac[num_CASES] >> 1399 `h < PRE n` by decide_tac >> 1400 rw_tac std_ss[] 1401 ]); 1402 1403(* Theorem: !k. 0 < k /\ k < n ==> (binomial n k) MOD n = 0 <=> 1404 !x y. SUM (GENLIST ((\k. (binomial n k * x ** (n - k) * y ** k) MOD n) o SUC) (PRE n)) = 0 *) 1405(* Proof: 1406 !k. 0 < k /\ k < n ==> (binomial n k) MOD n = 0 1407 <=> !k. 0 < k /\ k < n ==> !x y. ((binomial n k * x ** (n - k) * y ** k) MOD n = 0) by binomial_mod_zero 1408 <=> !h. h < PRE n ==> !x y. ((binomial n (SUC h) * x ** (n - (SUC h)) * y ** (SUC h)) MOD n = 0) by binomial_range_shift_alt 1409 <=> !x y. EVERY (\z. z = 0) (GENLIST (\k. (binomial n (SUC k) * x ** (n - (SUC k)) * y ** (SUC k)) MOD n) (PRE n)) by EVERY_GENLIST 1410 <=> !x y. EVERY (\x. x = 0) (GENLIST ((\k. binomial n k * x ** (n - k) * y ** k) o SUC) (PRE n) by FUN_EQ_THM 1411 <=> !x y. SUM (GENLIST ((\k. (binomial n k * x ** (n - k) * y ** k) MOD n) o SUC) (PRE n)) = 0 by SUM_EQ_0 1412*) 1413val binomial_mod_zero_alt = store_thm( 1414 "binomial_mod_zero_alt", 1415 ``!n. 0 < n ==> ((!k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0)) <=> 1416 !x y. SUM (GENLIST ((\k. (binomial n k * x ** (n - k) * y ** k) MOD n) o SUC) (PRE n)) = 0)``, 1417 rpt strip_tac >> 1418 `!x y. (\k. (binomial n (SUC k) * x ** (n - SUC k) * y ** (SUC k)) MOD n) = (\k. (binomial n k * x ** (n - k) * y ** k) MOD n) o SUC` by rw_tac std_ss[FUN_EQ_THM] >> 1419 `(!k. 0 < k /\ k < n ==> ((binomial n k) MOD n = 0)) <=> 1420 (!k. 0 < k /\ k < n ==> (!x y. ((binomial n k * x ** (n - k) * y ** k) MOD n = 0)))` by rw_tac std_ss[binomial_mod_zero] >> 1421 `_ = (!h. h < PRE n ==> (!x y. ((binomial n (SUC h) * x ** (n - (SUC h)) * y ** (SUC h)) MOD n = 0)))` by rw_tac std_ss[binomial_range_shift_alt] >> 1422 `_ = !x y h. h < PRE n ==> (((binomial n (SUC h) * x ** (n - (SUC h)) * y ** (SUC h)) MOD n = 0))` by metis_tac[] >> 1423 rw_tac std_ss[EVERY_GENLIST, SUM_EQ_0]); 1424 1425 1426(* ------------------------------------------------------------------------- *) 1427(* Binomial Theorem with prime exponent *) 1428(* ------------------------------------------------------------------------- *) 1429 1430 1431(* Theorem: [Binomial Expansion for prime exponent] (x + y)^p = x^p + y^p (mod p) *) 1432(* Proof: 1433 (x+y)^p (mod p) 1434 = SUM (k=0..p) C(p,k)x^(p-k)y^k (mod p) by binomial theorem 1435 = (C(p,0)x^py^0 + SUM (k=1..(p-1)) C(p,k)x^(p-k)y^k + C(p,p)x^0y^p) (mod p) by breaking sum 1436 = (x^p + SUM (k=1..(p-1)) C(p,k)x^(p-k)y^k + y^k) (mod p) by binomial_n_0, binomial_n_n 1437 = ((x^p mod p) + (SUM (k=1..(p-1)) C(p,k)x^(p-k)y^k) (mod p) + (y^p mod p)) mod p by MOD_PLUS 1438 = ((x^p mod p) + (SUM (k=1..(p-1)) (C(p,k)x^(p-k)y^k) (mod p)) + (y^p mod p)) mod p 1439 = (x^p mod p + 0 + y^p mod p) mod p by prime_iff_divides_binomials 1440 = (x^p + y^p) (mod p) by MOD_PLUS 1441*) 1442val binomial_thm_prime = store_thm( 1443 "binomial_thm_prime", 1444 ``!p. prime p ==> (!x y. (x + y) ** p MOD p = (x ** p + y ** p) MOD p)``, 1445 rpt strip_tac >> 1446 `0 < p` by rw_tac std_ss[PRIME_POS] >> 1447 `!k. 0 < k /\ k < p ==> ((binomial p k) MOD p = 0)` by metis_tac[prime_iff_divides_binomials, DIVIDES_MOD_0] >> 1448 `SUM (GENLIST ((\k. binomial p k * x ** (p - k) * y ** k) o SUC) (PRE p)) MOD p = 0` by metis_tac[SUM_GENLIST_MOD, binomial_mod_zero_alt, ZERO_MOD] >> 1449 `(x + y) ** p MOD p = (x ** p + SUM (GENLIST ((\k. binomial p k * x ** (p - k) * y ** k) o SUC) (PRE p)) + y ** p) MOD p` by rw_tac std_ss[binomial_thm, SUM_DECOMPOSE_FIRST_LAST, binomial_n_0, binomial_n_n, EXP] >> 1450 metis_tac[MOD_PLUS3, ADD_0, MOD_PLUS]); 1451 1452(* ------------------------------------------------------------------------- *) 1453 1454(* export theory at end *) 1455val _ = export_theory(); 1456 1457(*===========================================================================*) 1458