1(* ------------------------------------------------------------------------- *) 2(* Basic Computations *) 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 "computeBasic"; 12 13(* ------------------------------------------------------------------------- *) 14 15 16 17(* val _ = load "jcLib"; *) 18open jcLib; 19 20(* val _ = load "SatisfySimps"; (* for SatisfySimps.SATISFY_ss *) *) 21 22(* Get dependent theories local *) 23 24(* Get dependent theories in lib *) 25(* (* val _ = load "helperNumTheory"; -- in monoidTheory *) *) 26(* (* val _ = load "helperSetTheory"; -- in monoidTheory *) *) 27(* (* val _ = load "helperFunctionTheory"; -- in ringTheory *) *) 28(* (* val _ = load "helperListTheory"; -- in polyRingTheory *) *) 29(* val _ = load "logPowerTheory"; *) 30open helperNumTheory helperSetTheory helperListTheory helperFunctionTheory; 31open pred_setTheory listTheory arithmeticTheory; 32 33open logPowerTheory; (* for LOG2, SQRT, and Perfect Power, Power Free *) 34open logrootTheory; 35 36(* (* val _ = load "dividesTheory"; -- in helperNumTheory *) *) 37(* (* val _ = load "gcdTheory"; -- in helperNumTheory *) *) 38open dividesTheory gcdTheory; 39 40(* val _ = load "GaussTheory"; *) 41open EulerTheory; 42open GaussTheory; 43 44(* val _ = load "whileTheory"; *) 45open whileTheory; 46 47 48(* ------------------------------------------------------------------------- *) 49(* Basic Computations Documentation *) 50(* ------------------------------------------------------------------------- *) 51(* Datatype and Overloading: 52 sqrt_compute = root_compute 2 53*) 54(* Definitions and Theorems (# are exported, ! in computeLib): 55 56 LOG2 Computation: 57 log_compute_def |- !n. log_compute n = 58 if n = 0 then LOG2 0 59 else if n = 1 then 0 60 else 1 + log_compute (HALF n) 61 log_compute_eqn |- !n. log_compute n = LOG2 n 62 ulog_step_def |- !n m k. ulog_step n m k = 63 if m = 0 then 0 64 else if n <= m then k 65 else ulog_step n (TWICE m) (SUC k) 66 ulog_step_eq_count_up |- !n m k. ulog_step n m k = count_up n m k 67 ulog_compute_def |- !n. ulog_compute n = ulog_step n 1 0 68 ulog_compute_eqn |- !n. ulog_compute n = ulog n 69 70 71 Power Computation: 72 exp_binary_def |- !n m. exp_binary m n = 73 if n = 0 then 1 74 else (let r = exp_binary (m * m) (HALF n) in 75 if EVEN n then r else m * r) 76 exp_binary_0 |- !m. exp_binary m 0 = 1 77 exp_binary_1 |- !m. exp_binary m 1 = m 78 exp_binary_even |- !m n. EVEN n ==> (exp_binary m n = exp_binary (m * m) (HALF n)) 79 exp_binary_odd |- !m n. ODD n ==> (exp_binary m n = m * exp_binary (m * m) (HALF n)) 80 exp_binary_suc |- !m n. exp_binary m (SUC n) = m * exp_binary m n 81 exp_binary_eqn |- !m n. exp_binary m n = m ** n 82 exp_binary_of_0 |- !n. exp_binary 0 n = if n = 0 then 1 else 0 83 84 Fast Exponentiation: 85 exp_step_def |- !r n m. exp_step m n r = if n = 0 then r 86 else exp_step (SQ m) (HALF n) (if EVEN n then r else r * m) 87 exp_compute_def |- !m n. exp_compute m n = exp_step m n 1 88 exp_step_0 |- !m r. exp_step m 0 r = r 89 exp_step_1 |- !m r. exp_step m 1 r = r * m 90 exp_step_2 |- !m r. exp_step m 2 r = r * m * m 91 exp_step_even |- !n. EVEN n ==> !m r. exp_step m n r = exp_step (SQ m) (HALF n) r 92 exp_step_odd |- !n. ODD n ==> !m r. exp_step m n r = exp_step (SQ m) (HALF n) (r * m) 93 exp_step_eqn |- !m n r. exp_step m n r = r * m ** n 94 exp_compute_0 |- !m. exp_compute m 0 = 1 95 exp_compute_1 |- !m. exp_compute m 1 = m 96 exp_compute_2 |- !m. exp_compute m 2 = SQ m 97 exp_compute_even |- !n. EVEN n ==> !m. exp_compute m n = exp_compute (SQ m) (HALF n) 98 exp_compute_eqn |- !m n. exp_compute m n = m ** n 99 100 Modular Exponentiation: 101 exp_mod_binary_def |- !n m a. exp_mod_binary a n m = 102 if m = 0 then a ** n MOD 0 103 else if m = 1 then 0 104 else if n = 0 then 1 105 else (let b = exp_mod_binary ((a * a) MOD m) (HALF n) m 106 in if EVEN n then b else (a * b) MOD m) 107 exp_mod_binary_0 |- !a m. 0 < m ==> (exp_mod_binary a 0 m = if m = 1 then 0 else 1) 108 exp_mod_binary_1 |- !a m. exp_mod_binary a 1 m = a MOD m 109 exp_mod_binary_even |- !m n. 0 < m /\ EVEN n ==> 110 !a. exp_mod_binary a n m = exp_mod_binary ((a * a) MOD m) (HALF n) m 111 exp_mod_binary_odd |- !m n. 0 < m /\ ODD n ==> 112 !a. exp_mod_binary a n m = (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m 113 exp_mod_binary_suc |- !m n. 0 < m ==> !a. exp_mod_binary a (SUC n) m = (a * exp_mod_binary a n m) MOD m 114 exp_mod_binary_eqn |- !m n a. exp_mod_binary a n m = a ** n MOD m 115 exp_mod_binary_0_n |- !m n. exp_mod_binary 0 n m = (if n = 0 then 1 else 0) MOD m 116 117 Fast Modular Exponentiation: 118 exp_mod_step_def |- !r n m k. exp_mod_step m n k r = if k = 0 then (r * m ** n) MOD k 119 else if n = 0 then r MOD k 120 else exp_mod_step (SQ m MOD k) (HALF n) k 121 (if EVEN n then r else (r * m) MOD k) 122 exp_mod_compute_def |- !m n k. exp_mod_compute m n k = exp_mod_step m n k 1 123 exp_mod_step_0 |- !m r k. exp_mod_step m 0 k r = r MOD k 124 exp_mod_step_1 |- !m r k. exp_mod_step m 1 k r = (r * m) MOD k 125 exp_mod_step_2 |- !m r k. exp_mod_step m 2 k r = (r * m * m) MOD k 126 exp_mod_step_even |- !k. 0 < k ==> !n. EVEN n ==> 127 !m r. exp_mod_step m n k r = exp_mod_step (SQ m MOD k) (HALF n) k r 128 exp_mod_step_odd |- !k. 0 < k ==> !n. ODD n ==> 129 !m r. exp_mod_step m n k r = exp_mod_step (SQ m MOD k) (HALF n) k ((r * m) MOD k) 130 exp_mod_step_eqn |- !m n k r. exp_mod_step m n k r = (r * m ** n) MOD k 131 exp_mod_compute_0 |- !m k. exp_mod_compute m 0 k = 1 MOD k 132 exp_mod_compute_1 |- !m k. exp_mod_compute m 1 k = m MOD k 133 exp_mod_compute_2 |- !m k. exp_mod_compute m 2 k = SQ m MOD k 134 exp_mod_compute_even |- !k. 0 < k ==> !n. EVEN n ==> 135 !m. exp_mod_compute m n k = exp_mod_compute (SQ m MOD k) (HALF n) k 136 exp_mod_compute_odd |- !k. 0 < k ==> !n. ODD n ==> 137 !m. exp_mod_compute m n k = exp_mod_step (SQ m MOD k) (HALF n) k (m MOD k) 138 exp_mod_compute_eqn |- !m n k. exp_mod_compute m n k = m ** n MOD k 139 140 ROOT computation: 141 root_compute_def |- !r n. root_compute r n = 142 if 0 < r then 143 if n = 0 then 0 else 144 (let x = 2 * root_compute r (n DIV exp_compute 2 r) 145 in if n < exp_compute (SUC x) r then x else SUC x) 146 else ROOT 0 n 147 root_compute_alt |- !r n. root_compute r n = 148 if 0 < r then 149 if n = 0 then 0 else 150 (let x = 2 * root_compute r (n DIV 2 ** r) 151 in if n < SUC x ** r then x else SUC x) 152 else ROOT 0 n 153 root_compute_0 |- !r. 0 < r ==> (root_compute r 0 = 0) 154 root_compute_1 |- !n. root_compute 1 n = n 155 root_compute_eqn |- !r n. root_compute r n = ROOT r n 156 sqrt_compute_eqn |- !n. sqrt_compute n = SQRT n 157 158 Power Free Test: 159 power_index_compute_def |- !n k. power_index_compute n k = 160 if k <= 1 then 1 161 else if exp_compute (root_compute k n) k = n then k 162 else power_index_compute n (k - 1) 163 power_index_compute_alt |- !n k. power_index_compute n k = 164 if k <= 1 then 1 165 else if ROOT k n ** k = n then k 166 else power_index_compute n (k - 1) 167 power_index_compute_eqn |- !n k. power_index_compute n k = power_index n k 168 power_free_check_def |- !n. power_free_check n <=> 169 1 < n /\ (power_index_compute n (ulog_compute n) = 1) 170 power_free_check_0 |- power_free_check 0 <=> F 171 power_free_check_1 |- power_free_check 1 <=> F 172 power_free_check_alt |- !n. power_free_check n <=> 1 < n /\ (power_index n (ulog n) = 1) 173 power_free_check_eqn |- !n. power_free_check n <=> power_free n 174! power_free_eval |- !n. power_free n <=> power_free_check n 175 power_free_check_by_ulog |- !n. power_free_check n <=> 1 < n /\ (power_index n (ulog n) = 1) 176 power_free_check_by_LOG2 |- !n. power_free_check n <=> 1 < n /\ (power_index n (LOG2 n) = 1) 177 178 GCD computation: 179 ODD_IMP_GCD_CANCEL_EVEN |- !m n. ODD n ==> (gcd n (2 * m) = gcd n m) 180 BINARY_GCD |- !m n. (EVEN m /\ EVEN n ==> (gcd m n = 2 * gcd (HALF m) (HALF n))) /\ 181 (EVEN m /\ ODD n ==> (gcd m n = gcd (HALF m) n)) 182 gcd_compute_def |- !n m. gcd_compute n m = 183 if n = 0 then m 184 else if m = 0 then n 185 else if n = m then n 186 else if EVEN n then 187 if EVEN m then 2 * gcd_compute (HALF n) (HALF m) 188 else gcd_compute (HALF n) m 189 else if EVEN m then gcd_compute n (HALF m) 190 else if n < m then gcd_compute n (m - n) 191 else gcd_compute (n - m) m 192 gcd_compute_0 |- !n. (gcd_compute 0 n = n) /\ (gcd_compute n 0 = n) 193 gcd_compute_id |- !n. gcd_compute n n = n 194 gcd_compute_even_even |- !m n. EVEN n /\ EVEN m ==> (gcd_compute n m = 2 * gcd_compute (HALF n) (HALF m)) 195 gcd_compute_even_odd |- !m n. EVEN n /\ ODD m ==> (gcd_compute n m = gcd_compute (HALF n) m) 196 gcd_compute_odd_even |- !m n. ODD n /\ EVEN m ==> (gcd_compute n m = gcd_compute n (HALF m)) 197 gcd_compute_odd_odd |- !m n. ODD n /\ ODD m ==> 198 (gcd_compute n m = if n < m then gcd_compute n (m - n) else gcd_compute (n - m) m) 199 gcd_compute_eqn |- !m n. gcd_compute n m = gcd n m 200 gcd_compute_1 |- !n. (gcd_compute 1 n = 1) /\ (gcd_compute n 1 = 1) 201 gcd_compute_sym |- !m n. gcd_compute n m = gcd_compute m n 202 203 Phi Computation: 204 count_coprime_def |- !n j. count_coprime n j = 205 if j = 0 then 0 206 else if j = 1 then 1 207 else count_coprime n (j - 1) + if gcd_compute j n = 1 then 1 else 0 208 phi_compute_def |- !n. phi_compute n = count_coprime n n 209 phi_compute_0 |- phi_compute 0 = 0 210 phi_compute_1 |- phi_compute 1 = 1 211 count_coprime_0 |- !n. count_coprime n 0 = 0 212 count_coprime_1 |- !n. count_coprime n 1 = 1 213 count_coprime_alt |- !n j. count_coprime n j = 214 if j = 0 then 0 215 else if j = 1 then 1 216 else count_coprime n (j - 1) + if coprime j n then 1 else 0 217 count_coprime_suc |- !n k. count_coprime n (SUC k) = count_coprime n k + if coprime (SUC k) n then 1 else 0 218 count_coprime_eqn |- !n k. k <= n ==> (count_coprime n k = CARD (coprimes n INTER natural k)) 219 phi_compute_eqn |- !n. phi_compute n = phi n 220 221*) 222 223(* ------------------------------------------------------------------------- *) 224(* Helper Theorems *) 225(* ------------------------------------------------------------------------- *) 226 227(* ------------------------------------------------------------------------- *) 228(* Computations *) 229(* ------------------------------------------------------------------------- *) 230 231(* There are several computations to investigate: 232 233** LOG2 computation 234** Phi computation 235** Fast exponentiation 236** Fast modulo exponentiation 237** Order computation by fast exponentiation 238 239** Root computation 240** Perfect-power check by root computation 241** Special square-root computation? 242 243** Binary GCD computation 244** Euler phi computation by binary GCD computation 245** Coprime check by binary GCD computation? 246 247** Polynomial Addition under modulus (modulus is irrelevant) 248** Polynomial Multiplication under modulus (by rewriting) 249** Polynomial Introspective efficient check under modulus X ** k - |1| 250 251*) 252 253(* ------------------------------------------------------------------------- *) 254(* LOG2 Computation *) 255(* ------------------------------------------------------------------------- *) 256(* Pseudocode: 257 258Input: n with 0 < n 259Output: log2 n 260 261c <- 0 // initial count 262while (1 < n) { 263 c <- c + 1 // increment count 264 n <- n DIV 2 // reduce n by half 265} 266return c // the result 267 268*) 269 270(* 271Question: Can (countdivs n) be expressed in a WHILE loop? 272 Or, can we prove that (countdivs n) invokes itself at most (LOG2 n) times? 273*) 274 275(* Compute LOG2 by counting the number of divisions down to 1 *) 276val log_compute_def = Define ` 277 log_compute n = if n = 0 then LOG2 0 278 else if n = 1 then 0 279 else 1 + log_compute (HALF n) 280`; 281 282(* 283> EVAL ``log_compute 4``; --> 2 284> EVAL ``log_compute 5``; --> 2 285> EVAL ``log_compute 6``; --> 2 286> EVAL ``log_compute 7``; --> 2 287> EVAL ``log_compute 8``; --> 3 288> EVAL ``log_compute 9``; --> 3 289> EVAL ``log_compute 10``; --> 3 290> EVAL ``log_compute 1024``; --> 10 291*) 292 293(* Theorem: log_compute n = LOG2 n *) 294(* Proof: 295 By complete induction on n. 296 Apply log_compute_def, this is to show: 297 (1) 0 = LOG2 1, true by LOG_1 298 (1) If 1 < n ==> 1 + log_compute (HALF n) = LOG2 n 299 Note 2 <= n by 1 < n 300 so (0 + 1) * 2 <= n by arithmetic 301 Thus 0 < HALF n by X_LT_DIV, 0 < 2 302 LOG2 n 303 = 1 + LOG2 (HALF n) by LOG_DIV, 1 < 2 304 = 1 + log_compute (HALF n) by induction hypothesis, 0 < HALF n 305*) 306val log_compute_eqn = store_thm( 307 "log_compute_eqn", 308 ``!n. log_compute n = LOG2 n``, 309 completeInduct_on `n` >> 310 rw_tac std_ss[Once log_compute_def] >- 311 rw[LOG_1] >> 312 rw[LOG_DIV, X_LT_DIV]); 313 314(* For the cost: 315f(n) --> cost_f(n) 316add1(n) --> cost_add1(n) = log n 317half(n) --> cost_half(n) = 1 318log(n) --> cost_log(n) = ??? by log_compute_def 319 = cost_add1(n) + ... + cost_log(half n) 320*) 321 322(* Compute ulog by counting the number of doubles than equal or exceeds n *) 323(* This is just count_up renamed to ulog_step *) 324val ulog_step_def = tDefine "ulog_step" ` 325 ulog_step n m k = 326 if m = 0 then 0 (* just to provide m <> 0 for the next one *) 327 else if n <= m then k else ulog_step n (2 * m) (SUC k) 328` (WF_REL_TAC `measure (\(n, m, k). n - m)`); 329 330(* Theorem: ulog_step n m k = count_up n m k *) 331(* Proof: by ulog_step_def, count_up_def *) 332val ulog_step_eq_count_up = store_thm( 333 "ulog_step_eq_count_up", 334 ``!n m k. ulog_step n m k = count_up n m k``, 335 ntac 2 strip_tac >> 336 completeInduct_on `n - m` >> 337 simp[Once ulog_step_def] >> 338 rpt strip_tac >> 339 metis_tac[count_up_def]); 340 341(* Define upper LOG2 n by ulog_step *) 342val ulog_compute_def = Define` 343 ulog_compute n = ulog_step n 1 0 344`; 345 346(* 347> EVAL ``ulog_compute 4``; --> 2 348> EVAL ``ulog_compute 5``; --> 3 349> EVAL ``ulog_compute 6``; --> 3 350> EVAL ``ulog_compute 7``; --> 3 351> EVAL ``ulog_compute 8``; --> 3 352> EVAL ``ulog_compute 9``; --> 4 353> EVAL ``ulog_compute 10``; --> 4 354> EVAL ``ulog_compute 1024``; --> 10 355*) 356 357(* Theorem: ulog_compute n = ulog n *) 358(* Proof: by ulog_compute_def, ulog_def, ulog_step_eq_count_up *) 359val ulog_compute_eqn = store_thm( 360 "ulog_compute_eqn", 361 ``!n. ulog_compute n = ulog n``, 362 rw[ulog_compute_def, ulog_def, ulog_step_eq_count_up]); 363 364(* ------------------------------------------------------------------------- *) 365(* Power Computation (Fast Exponential Computation) *) 366(* ------------------------------------------------------------------------- *) 367(* Pseudocode: 368 369Input: m n 370Output: m ** n 371 372if (m == 1) return 1 // quick check for convenience only 373if (n == 1) return m // quick check for convenience only 374 375r <- 1 // initial result 376while (n <> 0) { 377 if (EVEN n) { // even exponent 378 n <- HALF n // reduce exponent by half 379 m <- m * m // update base by its square 380 } 381 else { // odd exponent 382 n <- n - 1 // make exponent even 383 r <- r * m // update result by multiplying with base 384 } 385} 386return r // the result 387 388Since the exponent n is decreasing (either by half or by one) in each iteration, the while-loop terminates. 389 390Recursive version: 391if (n = 0) then 1 392else (* precompute *) r <- (m * m) ** (HALF n) 393 (* now the result *) if EVEN n then r else m * r 394*) 395 396(* Define exponentiation by repeated squaring. *) 397val exp_binary_def = Define` 398 exp_binary (m:num) n = 399 if (n = 0) then 1 (* m ** 0 = 1 *) 400 (* either n or (n - 1) is EVEN, precompute the repeated square *) 401 else let r = exp_binary (m * m) (HALF n) in 402 if EVEN n then r else m * r (* if EVEN, return reuslt, otherwise multiply by base for ODD *) 403`; 404 405(* 406> EVAL ``exp_binary 3 2``; --> 9 407> EVAL ``exp_binary 3 3``; --> 27 408> EVAL ``exp_binary 3 4``; --> 81 409> EVAL ``exp_binary 3 5``; --> 243 410> EVAL ``exp_binary 2 10``; --> 1024 411> EVAL ``exp_binary 3 10``; --> 59049 412*) 413 414(* Theorem: exp_binary m 0 = 1 *) 415(* Proof: by exp_binary_def *) 416val exp_binary_0 = store_thm( 417 "exp_binary_0", 418 ``!m. exp_binary m 0 = 1``, 419 rw[Once exp_binary_def]); 420 421(* Theorem: exp_binary m 1 = m *) 422(* Proof: by exp_binary_def *) 423val exp_binary_1 = store_thm( 424 "exp_binary_1", 425 ``!m. exp_binary m 1 = m``, 426 rw[Once exp_binary_def] >> 427 rw[Once exp_binary_def]); 428 429(* Theorem: EVEN n ==> (exp_binary m n = exp_binary (m * m) (HALF n)) *) 430(* Proof: 431 If n = 0, true by exp_binary_0, HALF 0 = 0. 432 If n <> 0, true by exp_binary_def 433*) 434val exp_binary_even = store_thm( 435 "exp_binary_even", 436 ``!m n. EVEN n ==> (exp_binary m n = exp_binary (m * m) (HALF n))``, 437 rpt strip_tac >> 438 Cases_on `n = 0` >- 439 rw[exp_binary_0] >> 440 rw[Once exp_binary_def]); 441 442(* Theorem: ODD n ==> (exp_binary m n = m * exp_binary (m * m) (HALF n)) *) 443(* Proof: 444 Note ODD 0 = F by ODD 445 and ODD n = ~(EVEN n) by ODD_EVEN 446 Thus true by exp_binary_def 447*) 448val exp_binary_odd = store_thm( 449 "exp_binary_odd", 450 ``!m n. ODD n ==> (exp_binary m n = m * exp_binary (m * m) (HALF n))``, 451 rw_tac std_ss[Once exp_binary_def] >> 452 metis_tac[ODD_EVEN]); 453 454(* Theorem: exp_binary m (SUC n) = m * exp_binary m n *) 455(* Proof: 456 By complete induction on n. 457 If EVEN n, 458 Then ODD (SUC n) by EVEN_ODD_SUC 459 exp_binary m (SUC n) 460 = m * exp_binary (m * m) (HALF (SUC n)) by exp_binary_odd 461 = m * exp_binary (m * m) (HALF n) by EVEN_SUC_HALF 462 = m * exp_binary m n by exp_binary_even 463 464 If ~(EVEN n), then ODD n by EVEN_ODD 465 Then EVEN (SUC n) by EVEN_ODD_SUC 466 Note n <> 0 by ODD, ODD 0 = F. 467 Thus HALF n < n by DIV_LESS, 0 < n, 1 < 2. 468 m * exp_binary m n 469 = m * m * exp_binary (m * m) (HALF n) by exp_binary_odd 470 = exp_binary (m * m) (SUC (HALF n)) by induction hypothesis, HALF n < n. 471 = exp_binary (m * m) (HALF (SUC n)) by ODD_SUC_HALF 472 = exp_binary m (SUC n) by exp_binary_even 473*) 474val exp_binary_suc = store_thm( 475 "exp_binary_suc", 476 ``!m n. exp_binary m (SUC n) = m * exp_binary m n``, 477 completeInduct_on `n` >> 478 rpt strip_tac >> 479 Cases_on `EVEN n` >| [ 480 `ODD (SUC n)` by rw[GSYM EVEN_ODD_SUC] >> 481 `exp_binary m (SUC n) = m * exp_binary (m * m) (HALF (SUC n))` by rw[exp_binary_odd] >> 482 `_ = m * exp_binary (m * m) (HALF n)` by rw[EVEN_SUC_HALF] >> 483 `_ = m * exp_binary m n` by rw[exp_binary_even] >> 484 rw[], 485 `ODD n` by metis_tac[EVEN_ODD] >> 486 `EVEN (SUC n)` by rw[GSYM EVEN_ODD_SUC] >> 487 `n <> 0` by metis_tac[ODD] >> 488 `HALF n < n` by rw[] >> 489 `m * exp_binary m n = m * m * exp_binary (m * m) (HALF n)` by rw[exp_binary_odd] >> 490 `_ = exp_binary (m * m) (SUC (HALF n))` by rw[] >> 491 `_ = exp_binary (m * m) (HALF (SUC n))` by rw[ODD_SUC_HALF] >> 492 `_ = exp_binary m (SUC n)` by rw[exp_binary_even] >> 493 rw[] 494 ]); 495 496(* Theorem: exp_binary m n = m ** n *) 497(* Proof: 498 By induction on n. 499 Base: exp_binary m 0 = m ** 0 500 exp_binary m 0 501 = 1 by exp_binary_0 502 = m ** 0 by EXP 503 Step: exp_binary m n = m ** n ==> exp_binary m (SUC n) = m ** SUC n 504 exp_binary m (SUC n) 505 = m * exp_binary m n by exp_binary_suc 506 = m * m ** n by induction hypothesis 507 = m ** SUC n by EXP 508*) 509val exp_binary_eqn = store_thm( 510 "exp_binary_eqn", 511 ``!m n. exp_binary m n = m ** n``, 512 rpt strip_tac >> 513 Induct_on `n` >- 514 rw[exp_binary_0] >> 515 rw[exp_binary_suc, EXP]); 516 517(* Theorem: exp_binary 0 n = if n = 0 then 1 else 0 *) 518(* Proof: 519 exp_binary 0 n 520 = 0 ** n by exp_binary_eqn 521 = if n = 0 then 1 else 0 by ZERO_EXP 522*) 523val exp_binary_of_0 = store_thm( 524 "exp_binary_of_0", 525 ``!n. exp_binary 0 n = if n = 0 then 1 else 0``, 526 rw[exp_binary_eqn]); 527 528(* ------------------------------------------------------------------------- *) 529(* Fast Exponentiation (Repeated Squares Method) *) 530(* ------------------------------------------------------------------------- *) 531 532(* Fast Exponentiation -- iterative form *) 533(* Pseudo-Code: 534 m ** n = r <- 1 535 loop: 536 if (n == 0) return r 537 if (ODD n) r <- r * m 538 n <- HALF n 539 m <- SQ m 540 goto loop 541*) 542 543(* Define fast exponentiation *) 544val exp_step_def = Define` 545 exp_step m n r = (* r = m ** n *) 546 if n = 0 then r else 547 exp_step (SQ m) (HALF n) (if EVEN n then r else r * m) 548`; 549val exp_compute_def = Define` 550 exp_compute m n = exp_step m n 1 551`; 552 553(* 554EVAL ``exp_compute 2 10``; --> 1024 555EVAL ``exp_compute 3 10``; --> 59049 556EVAL ``exp_compute 3 8 = 3 ** 8``; --> T 557*) 558 559(* Theorem: exp_step m 0 r = r *) 560(* Proof: by exp_step_def *) 561val exp_step_0 = store_thm( 562 "exp_step_0", 563 ``!m r. exp_step m 0 r = r``, 564 rw[Once exp_step_def]); 565 566(* Theorem: exp_step m 1 r = r * m *) 567(* Proof: by exp_step_def *) 568val exp_step_1 = store_thm( 569 "exp_step_1", 570 ``!m r. exp_step m 1 r = r * m``, 571 rw[Once exp_step_def, Once exp_step_def]); 572 573(* Theorem: exp_step m 2 r = r * m * m *) 574(* Proof: 575 exp_step m 2 r 576 = exp_step (SQ m) (HALF 2) r by exp_step_def, EVEN 2 577 = exp_step (SQ m) 1 r by arithmetic 578 = r * SQ m by exp_step_1 579 = r * (m * m) by EXP_2 580 = r * m * m by MULT_ASSOC 581*) 582val exp_step_2 = store_thm( 583 "exp_step_2", 584 ``!m r. exp_step m 2 r = r * m * m``, 585 rw[Once exp_step_def, Once exp_step_def, Once exp_step_def]); 586 587(* Theorem: EVEN n ==> !m r. exp_step m n r = exp_step (SQ m) (HALF n) r *) 588(* Proof: 589 If n = 0, HALF 0 = 0 by HALF_EQ_0 590 Thus LHS = r = RHS by exp_step_def 591 If n <> 0, true by exp_step_def 592*) 593val exp_step_even = store_thm( 594 "exp_step_even", 595 ``!n. EVEN n ==> !m r. exp_step m n r = exp_step (SQ m) (HALF n) r``, 596 rpt strip_tac >> 597 Cases_on `n = 0` >- 598 metis_tac[exp_step_def, HALF_EQ_0] >> 599 rw[Once exp_step_def]); 600 601(* Theorem: ODD n ==> !m r. exp_step m n r = exp_step (SQ m) (HALF n) (r * m) *) 602(* Proof: 603 Note ODD n ==> ~EVEN n by ODD_EVEN 604 and n <> 0 by EVEN 605 The result follows by exp_step_def 606*) 607val exp_step_odd = store_thm( 608 "exp_step_odd", 609 ``!n. ODD n ==> !m r. exp_step m n r = exp_step (SQ m) (HALF n) (r * m)``, 610 rpt strip_tac >> 611 `~EVEN n /\ n <> 0` by metis_tac[ODD_EVEN, EVEN] >> 612 rw[Once exp_step_def]); 613 614(* Theorem: exp_step m n r = r * m ** n *) 615(* Proof: 616 By complete induction on n. 617 If n = 0, 618 exp_step m 0 r 619 = r by exp_step_0 620 = r * 1 by MULT_RIGHT_1 621 = r * m ** 0 by EXP_0 622 If n <> 0, 623 then HALF n < n by HALF_LESS, 0 < n 624 If EVEN n, 625 exp_step m n r 626 = exp_step (SQ m) (HALF n) r by exp_step_even 627 = r * (SQ m) ** (HALF n) by induction hypothesis, HALF n < n 628 = r * m ** n by EXP_EVEN 629 If ~EVEN n, 630 then ODD n by EVEN_ODD 631 exp_step m n r 632 = exp_step (SQ m) (HALF n) (r * m) by exp_step_odd 633 = (r * m) * (SQ m) ** (HALF n) by induction hypothesis, HALF n < n 634 = r * (m * (SQ m) ** (HALF n)) by MULT_ASSOC 635 = r * m ** n by EXP_ODD 636*) 637val exp_step_eqn = store_thm( 638 "exp_step_eqn", 639 ``!m n r. exp_step m n r = r * m ** n``, 640 completeInduct_on `n` >> 641 rpt strip_tac >> 642 Cases_on `n = 0` >- 643 rw[exp_step_0] >> 644 `HALF n < n` by rw[] >> 645 Cases_on `EVEN n` >- 646 rw[exp_step_even, GSYM EXP_EVEN] >> 647 metis_tac[ODD_EVEN, exp_step_odd, MULT_ASSOC, EXP_ODD]); 648 649(* Theorem: exp_compute m 0 = 1 *) 650(* Proof: by exp_compute_def, exp_step_0 *) 651val exp_compute_0 = store_thm( 652 "exp_compute_0", 653 ``!m. exp_compute m 0 = 1``, 654 rw[exp_compute_def, exp_step_0]); 655 656(* Theorem: exp_compute m 1 = m *) 657(* Proof: by exp_compute_def, exp_step_1 *) 658val exp_compute_1 = store_thm( 659 "exp_compute_1", 660 ``!m. exp_compute m 1 = m``, 661 rw[exp_compute_def, exp_step_1]); 662 663(* Theorem: exp_compute m 2 = SQ m *) 664(* Proof: by exp_compute_def, exp_step_2 *) 665val exp_compute_2 = store_thm( 666 "exp_compute_2", 667 ``!m. exp_compute m 2 = SQ m``, 668 rw[exp_compute_def, exp_step_2]); 669 670(* Theorem: EVEN n ==> !m. exp_compute m n = exp_compute (SQ m) (HALF n) *) 671(* Proof: 672 exp_compute m n 673 = exp_step m n 1 by exp_compute_def 674 = exp_step (SQ m) (HALF n) 1 by exp_step_even, EVEN n 675 = exp_compute (SQ m) (HALF n) by exp_compute_def 676*) 677val exp_compute_even = store_thm( 678 "exp_compute_even", 679 ``!n. EVEN n ==> !m. exp_compute m n = exp_compute (SQ m) (HALF n)``, 680 rw[exp_compute_def, exp_step_even]); 681 682(* Theorem: exp_compute m n = m ** n *) 683(* Proof: 684 exp_compute m n 685 = exp_step m n 1 by exp_compute_def 686 = 1 * m ** n by exp_step_eqn 687 = m ** n by MULT_LEFT_1 688*) 689val exp_compute_eqn = store_thm( 690 "exp_compute_eqn", 691 ``!m n. exp_compute m n = m ** n``, 692 rw[exp_compute_def, exp_step_eqn]); 693 694(* ------------------------------------------------------------------------- *) 695(* Modular Exponentiation. *) 696(* ------------------------------------------------------------------------- *) 697(* Pseudocode: 698 699Input: b n m 700Output: (b ** n) MOD m 701 702Recursive version: 703if (n = 0) then return 1 704if EVEN n then return (b * b MOD m) ** (HALF n) 705else return b * ((b * b MOD m) ** (HALF n)) MOD m 706 707val EXP_MOD_EQN = 708 |- !b n m. 1 < m ==> 709 b ** n MOD m = 710 if n = 0 then 1 711 else (let result = SQ b ** HALF n MOD m 712 in if EVEN n then result else (b * result) MOD m): thm 713*) 714 715(* Pseudocode: 716 717Input: a n m 718Output: (a ** n) MOD m 719 720Iterative version: 721 722r <- 1 // initial result 723while (n <> 0) { 724 if (EVEN n) { // even exponent 725 n <- HALF n // reduce exponent by half 726 a <- a * a MOD m // update base by its square 727 r <- a // update result by new base 728 } 729 else { // odd exponent 730 b <- a // save current base 731 n <- n - 1 // make exponent even 732 n <- HALF n // reduce exponent by half 733 a <- a * a MOD m // update base by its square 734 r <- b * a MOD m // update result by multiplying with base 735 } 736} 737return r // the result 738 739Since the exponent n is decreasing by half in each iteration, the while-loop terminates. 740 741*) 742 743(* Smart definition of: a ** n (mod m) by repeated squaring *) 744(* 745val exp_mod_binary_def = Define` 746 exp_mod_binary a n m = 747 if m = 0 then (a ** n) MOD 0 (* whatever that is! *) 748 else if m = 1 then 0 (* all MOD 1 = 0 *) 749 else if n = 0 then 1 (* a ** 0 (mod m) = 1 *) 750 (* if EVEN n then a ** n (mod m) = (a * a MOD m) ** (HALF n) *) 751 else if EVEN n then exp_mod_binary ((a * a) MOD m) (HALF n) m 752 (* if ODD n then a ** n (mod m) = a * ((a * a MOD m) ** (HALF n)) MOD m *) 753 else (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m 754`; 755*) 756val exp_mod_binary_def = Define` 757 exp_mod_binary a n m = 758 if m = 0 then (a ** n) MOD 0 (* whatever that is! *) 759 else if m = 1 then 0 (* all MOD 1 = 0 *) 760 else if n = 0 then 1 (* a ** 0 (mod m) = 1 *) 761 (* if EVEN n then a ** n (mod m) = (a * a MOD m) ** (HALF n) *) 762 (* if ODD n then a ** n (mod m) = a * ((a * a MOD m) ** (HALF n)) MOD m *) 763 else (let b = exp_mod_binary ((a * a) MOD m) (HALF n) m 764 in if EVEN n then b else (a * b) MOD m) 765`; 766(* Note the for-all order: n m a, which is good. 767val it = 768 |- !n m a. exp_mod_binary a n m = 769 if m = 0 then a ** n MOD 0 770 else if m = 1 then 0 771 else if n = 0 then 1 772 else (let b = exp_mod_binary ((a * a) MOD m) (HALF n) m in if EVEN n then b else (a * b) MOD m): thm 773*) 774(* 775- type_of ``exp_mod_binary a n m``; 776> val it = ``:num`` : hol_type 777- type_of ``exp_mod_binary``; 778> val it = ``:num -> num -> num -> num`` : hol_type 779*) 780 781(* 782> EVAL ``exp_mod_binary 3 0 10``; --> 1 783> EVAL ``exp_mod_binary 3 1 10``; --> 3 784> EVAL ``exp_mod_binary 3 2 10``; --> 9 785> EVAL ``exp_mod_binary 3 3 10``; --> 7 786> EVAL ``exp_mod_binary 3 4 10``; --> 1 787*) 788 789 790(* Theorem: 0 < m ==> (exp_mod_binary a 0 m = if m = 1 then 0 else 1) *) 791(* Proof: by exp_mod_binary_def. *) 792val exp_mod_binary_0 = store_thm( 793 "exp_mod_binary_0", 794 ``!a m. 0 < m ==> (exp_mod_binary a 0 m = if m = 1 then 0 else 1)``, 795 rw[Once exp_mod_binary_def]); 796 797(* Theorem: exp_mod_binary a 1 m = a MOD m *) 798(* Proof: by exp_mod_binary_def, exp_mod_binary_0. *) 799val exp_mod_binary_1 = store_thm( 800 "exp_mod_binary_1", 801 ``!a m. exp_mod_binary a 1 m = a MOD m``, 802 rw[Once exp_mod_binary_def] >> 803 rw[exp_mod_binary_0]); 804 805(* Theorem: 0 < m /\ EVEN n ==> !a. exp_mod_binary a n m = exp_mod_binary ((a * a) MOD m) (HALF n) m *) 806(* Proof: 807 If n = 0, 808 exp_mod_binary a 0 m 809 = 1 by exp_mod_binary_0 810 = exp_mod_binary (((a * a) MOD m)) 0 m by exp_mod_binary_0 811 = exp_mod_binary (((a * a) MOD m)) (HALF 0) m by ZERO_MOD, 0 < 2 812 If n <> 0, 813 If m = 1, 814 exp_mod_binary a n 1 815 = 0 by exp_mod_binary_def 816 = exp_mod_binary ((a * a) MOD 1) (HALF n) 1 by exp_mod_binary_def 817 If m <> 1, true for EVEN n by exp_mod_binary_def 818*) 819val exp_mod_binary_even = store_thm( 820 "exp_mod_binary_even", 821 ``!m n. 0 < m /\ EVEN n ==> !a. exp_mod_binary a n m = exp_mod_binary ((a * a) MOD m) (HALF n) m``, 822 rpt strip_tac >> 823 Cases_on `n = 0` >- 824 simp[exp_mod_binary_0] >> 825 rw[Once exp_mod_binary_def] >> 826 rw[Once exp_mod_binary_def]); 827 828(* Theorem: 0 < m /\ ODD n ==> !a. exp_mod_binary a n m = (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m *) 829(* Proof: 830 If m = 1, 831 exp_mod_binary a n 1 832 = 0 by exp_mod_binary_def 833 = (a * exp_mod_binary ((a * a) MOD 1) (HALF n) 1) MOD 1 by MOD_1 834 If m <> 1, 835 and ODD n <=> ~(EVEN n) by ODD_EVEN 836 so n <> 0 by EVEN 837 Thus true for ~(EVEN n) by exp_mod_binary_def 838*) 839val exp_mod_binary_odd = store_thm( 840 "exp_mod_binary_odd", 841 ``!m n. 0 < m /\ ODD n ==> !a. exp_mod_binary a n m = (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m``, 842 rpt strip_tac >> 843 Cases_on `m = 1` >- 844 rw[Once exp_mod_binary_def] >> 845 rpt strip_tac >> 846 `~EVEN n` by rw[GSYM ODD_EVEN] >> 847 `n <> 0` by metis_tac[EVEN] >> 848 rw[Once exp_mod_binary_def]); 849 850(* Theorem: 0 < m ==> !a. exp_mod_binary a (SUC n) m = (a * exp_mod_binary a n m) MOD m *) 851(* Proof: 852 By complete induction on n. 853 If EVEN n, 854 Then ODD (SUC n) by EVEN_ODD_SUC 855 exp_mod_binary a (SUC n) m 856 = (a * exp_mod_binary ((a * a) MOD m) (HALF (SUC n)) m) MOD m by exp_mod_binary_odd 857 = (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m by EVEN_SUC_HALF 858 = (a * exp_mod_binary a n m) MOD m by exp_mod_binary_even 859 If ~(EVEN n), then ODD n. by ODD_EVEN 860 Note ODD n ==> n <> 0 by ODD 861 Thus HALF n < n by DIV_LESS, 0 < n, 1 < 2 862 and EVEN (SUC n) by EVEN_ODD_SUC 863 (a * exp_mod_binary a n m) MOD m 864 = (a * (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m) MOD m by exp_mod_binary_odd, 0 < m 865 = (a * (a * exp_mod_binary ((a * a) MOD m) (HALF n) m)) MOD m by MOD_MOD, MOD_TIMES2, 0 < m 866 = (a * a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m by MULT_ASSOC 867 = (((a * a) MOD m) * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m by MOD_MOD, MOD_TIMES2, 0 < m 868 = exp_mod_binary ((a * a) MOD m) (SUC (HALF n)) m by induction hypothesis 869 = exp_mod_binary ((a * a) MOD m) (HALF (SUC n)) m by ODD_SUC_HALF 870 = exp_mod_binary a (SUC n) m by exp_mod_binary_even, 0 < m 871*) 872val exp_mod_binary_suc = store_thm( 873 "exp_mod_binary_suc", 874 ``!m n. 0 < m ==> !a. exp_mod_binary a (SUC n) m = (a * exp_mod_binary a n m) MOD m``, 875 ntac 3 strip_tac >> 876 completeInduct_on `n` >> 877 rpt strip_tac >> 878 Cases_on `EVEN n` >| [ 879 `ODD (SUC n)` by metis_tac[EVEN_ODD_SUC] >> 880 `exp_mod_binary a (SUC n) m = (a * exp_mod_binary ((a * a) MOD m) (HALF (SUC n)) m) MOD m` by rw[exp_mod_binary_odd] >> 881 `_ = (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m` by rw[EVEN_SUC_HALF] >> 882 `_ = (a * exp_mod_binary a n m) MOD m` by rw[exp_mod_binary_even] >> 883 rw[], 884 `ODD n` by rw[ODD_EVEN] >> 885 `n <> 0` by metis_tac[ODD] >> 886 `HALF n < n` by rw[DIV_LESS] >> 887 `EVEN (SUC n)` by rw[GSYM EVEN_ODD_SUC] >> 888 `(a * exp_mod_binary a n m) MOD m = (a * (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m) MOD m` by rw[exp_mod_binary_odd] >> 889 `_ = ((a * a) MOD m * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m` by metis_tac[MOD_MOD, MOD_TIMES2, MULT_ASSOC] >> 890 `_ = exp_mod_binary ((a * a) MOD m) (SUC (HALF n)) m` by rw[] >> 891 `_ = exp_mod_binary ((a * a) MOD m) (HALF (SUC n)) m` by rw[ODD_SUC_HALF] >> 892 `_ = exp_mod_binary a (SUC n) m` by rw[exp_mod_binary_even] >> 893 rw[] 894 ]); 895 896(* Theorem: exp_mod_binary a n m = a ** n MOD m *) 897(* Proof: 898 If m = 0, true trivially by exp_mod_binary_def 899 If m = 1, 900 exp_mod_binary a n m 901 = 1 by exp_mod_binary_def 902 = (a ** n) MOD 1 by MOD_1 903 If m <> 0 and m <> 1, 904 Then 1 < m. 905 By induction on n. 906 Base: exp_mod_binary a 0 m = a ** 0 MOD m 907 exp_mod_binary a 0 m 908 = 1 by exp_mod_binary_0 909 = 1 MOD m by ONE_MOD, 1 < m 910 = a ** 0 MOD m by EXP 911 Step: exp_mod_binary a n m = a ** n MOD m ==> exp_mod_binary a (SUC n) m = a ** SUC n MOD m 912 exp_mod_binary a (SUC n) m 913 = (a * exp_mod_binary a n m) MOD m by exp_mod_binary_suc, 0 < m 914 = (a * a ** n MOD m) MOD m by induction hypothesis 915 = (a * a ** n) MOD m by MOD_MOD, MOD_TIMES2, 0 < m 916 = (a ** SUC n) MOD m by EXP 917*) 918val exp_mod_binary_eqn = store_thm( 919 "exp_mod_binary_eqn", 920 ``!m n a. exp_mod_binary a n m = a ** n MOD m``, 921 rpt strip_tac >> 922 Cases_on `m = 0` >- 923 rw[Once exp_mod_binary_def] >> 924 Cases_on `m = 1` >- 925 rw[Once exp_mod_binary_def] >> 926 `1 < m /\ 0 < m` by decide_tac >> 927 Induct_on `n` >- 928 simp[exp_mod_binary_0, ONE_MOD, EXP] >> 929 metis_tac[exp_mod_binary_suc, MOD_MOD, MOD_TIMES2, EXP]); 930 931(* Another proof of the same result *) 932 933(* Theorem: exp_mod_binary a n m = (a ** n) MOD m *) 934(* Proof: 935 If m = 0, true trivially by exp_mod_binary_def 936 If m = 1, 937 exp_mod_binary a n m 938 = 1 by exp_mod_binary_def 939 = (a ** n) MOD 1 by MOD_1 940 If m <> 0 and m <> 1, 941 Then 1 < m. 942 By complete induction on n. 943 Assume: !j. j < n ==> !a. exp_mod_binary a j m = a ** j MOD m 944 To show: exp_mod_binary a n m = a ** n MOD m 945 If n = 0, 946 exp_mod_binary a 0 m 947 = 1 by exp_mod_binary_0 948 = 1 MOD m by ONE_MOD, 1 < m 949 = (a ** 0) MOD m by EXP 950 If n <> 0, 951 Then HALF n < n by HALF_LESS 952 If EVEN n, 953 Then n MOD 2 = 0 by EVEN_MOD2 954 exp_mod_binary a n m 955 = exp_mod_binary ((a * a) MOD m) (HALF n) m by exp_mod_binary_def, n MOD 2 = 0 956 = exp_mod_binary ((a ** 2) MOD m) (HALF n) m by EXP_2 957 = (((a ** 2) MOD m) ** (HALF n)) MOD m by induction hypothesis, HALF n < n 958 = ((a ** 2) ** (HALF n)) MOD m by EXP_MOD, 0 < m 959 = (a ** (2 * HALF n)) MOD m by EXP_EXP_MULT 960 = (a ** n) MOD m by EVEN_HALF 961 If ~EVEN n, 962 Then ODD n by ODD_EVEN 963 exp_mod_binary a n m 964 = (a * exp_mod_binary ((a * a) MOD m) (HALF n) m) MOD m by exp_mod_binary_def, n MOD 2 <> 0 965 = (a * exp_mod_binary ((a ** 2) MOD m) (HALF n) m) MOD m by EXP_2 966 = (a * (((a ** 2) MOD m) ** (HALF n)) MOD m) MOD m by induction hypothesis, HALF n < n 967 = (a * ((a ** 2) ** (HALF n)) MOD m) MOD m by EXP_MOD, 0 < m 968 = (a * (a ** (2 * HALF n)) MOD m) MOD m by EXP_EXP_MULT 969 = (a * a ** (2 * HALF n)) MOD m by MOD_TIMES2, 0 < m 970 = (a ** (1 + 2 * HALF n)) MOD m by EXP_ADD 971 = (a ** (2 * HALF n + 1)) MOD m by arithmetic 972 = (a ** n) MOD m by ODD_HALF 973*) 974val exp_mod_binary_eqn = store_thm( 975 "exp_mod_binary_eqn", 976 ``!m n a. exp_mod_binary a n m = (a ** n) MOD m``, 977 ntac 2 strip_tac >> 978 Cases_on `m = 0` >- 979 rw[Once exp_mod_binary_def] >> 980 Cases_on `m = 1` >- 981 rw[Once exp_mod_binary_def] >> 982 `1 < m /\ 0 < m` by decide_tac >> 983 completeInduct_on `n` >> 984 rpt strip_tac >> 985 Cases_on `n = 0` >- 986 rw[exp_mod_binary_0, EXP] >> 987 `0 < m` by decide_tac >> 988 `HALF n < n` by rw[HALF_LESS] >> 989 rw[Once exp_mod_binary_def] >| [ 990 `((a ** 2) ** HALF n) MOD m = (a ** (2 * HALF n)) MOD m` by rw[EXP_EXP_MULT] >> 991 `_ = (a ** n) MOD m` by rw[GSYM EVEN_HALF, EVEN_MOD2] >> 992 rw[], 993 `ODD n` by rw[ODD_EVEN] >> 994 `(a * (a ** 2) ** HALF n) MOD m = (a * (a ** (2 * HALF n) MOD m)) MOD m` by rw[EXP_EXP_MULT] >> 995 `_ = (a * a ** (2 * HALF n)) MOD m` by metis_tac[MOD_TIMES2, MOD_MOD] >> 996 `_ = (a ** (2 * HALF n + 1)) MOD m` by rw[EXP_ADD] >> 997 `_ = a ** n MOD m` by metis_tac[ODD_HALF] >> 998 rw[] 999 ]); 1000 1001(* Theorem: exp_mod_binary 0 n m = (if n = 0 then 1 else 0) MOD m *) 1002(* Proof: 1003 exp_mod_binary 0 n m 1004 = (0 ** n) MOD m by exp_mod_binary_eqn 1005 = (if n = 0 then 1 else 0) MOD m by ZERO_EXP 1006*) 1007val exp_mod_binary_0_n = store_thm( 1008 "exp_mod_binary_0_n", 1009 ``!m n. exp_mod_binary 0 n m = (if n = 0 then 1 else 0) MOD m``, 1010 rw[exp_mod_binary_eqn, ZERO_EXP]); 1011 1012(* ------------------------------------------------------------------------- *) 1013(* Fast Modular Exponentiation *) 1014(* ------------------------------------------------------------------------- *) 1015 1016(* Fast Modular Exponentiation -- iterative form *) 1017(* Pseudo-Code: 1018 (m ** n) MOD k = r <- 1 1019 loop: 1020 if (n == 0) return (r MOD k) 1021 if (ODD n) r <- (r * m) MOD k 1022 n <- HALF n 1023 m <- (SQ m) MOD k 1024 goto loop 1025*) 1026 1027(* Define fast modulo exponentiation *) 1028val exp_mod_step_def = Define` 1029 exp_mod_step m n k r = 1030 if k = 0 then (r * m ** n) MOD k 1031 else if n = 0 then (r MOD k) else 1032 exp_mod_step ((SQ m) MOD k) (HALF n) k (if EVEN n then r else (r * m) MOD k) 1033`; 1034val exp_mod_compute_def = Define` 1035 exp_mod_compute m n k = exp_mod_step m n k 1 1036`; 1037 1038(* 1039EVAL ``exp_mod_compute 2 10 3``; --> 1024 MOD 3 = 1 1040EVAL ``exp_mod_compute 3 10 17``; --> 59049 MOD 17 = 8 1041EVAL ``exp_mod_compute 3 8 19 = (3 ** 8) MOD 19``; --> T 1042*) 1043 1044(* Theorem: exp_mod_step m 0 k r = r MOD k *) 1045(* Proof: by exp_mod_step_def *) 1046val exp_mod_step_0 = store_thm( 1047 "exp_mod_step_0", 1048 ``!m r k. exp_mod_step m 0 k r = r MOD k``, 1049 rw[Once exp_mod_step_def]); 1050 1051(* Theorem: exp_mod_step m 1 k r = (r * m) MOD k *) 1052(* Proof: by exp_mod_step_def *) 1053val exp_mod_step_1 = store_thm( 1054 "exp_mod_step_1", 1055 ``!m r k. exp_mod_step m 1 k r = (r * m) MOD k``, 1056 rw[Once exp_mod_step_def, Once exp_mod_step_def]); 1057 1058(* Theorem: exp_mod_step m 2 k r = (r * m * m) MOD k *) 1059(* Proof: 1060 exp_mod_step m 2 k r 1061 = exp_mod_step ((SQ m) MOD k) (HALF 2) k r by exp_mod_step_def, EVEN 2 1062 = exp_mod_step ((SQ m) MOD k) 1 r by arithmetic 1063 = (r * (SQ m) MOD k) MOD k by exp_mod_step_1 1064 = (r * (SQ m)) MOD k by MOD_TIMES2, MOD_MOD 1065 = (r * (m * m)) MOD k by EXP_2 1066 = (r * m * m) MOD k by MULT_ASSOC 1067*) 1068val exp_mod_step_2 = store_thm( 1069 "exp_mod_step_2", 1070 ``!m r k. exp_mod_step m 2 k r = (r * m * m) MOD k``, 1071 rw[Once exp_mod_step_def, Once exp_mod_step_def, Once exp_mod_step_def]); 1072 1073(* Theorem: 0 < k ==> !n. EVEN n ==> 1074 !m r. exp_mod_step m n k r = exp_mod_step ((SQ m) MOD k) (HALF n) k r *) 1075(* Proof: 1076 if n = 0, HALF 0 = 0 by HALF_EQ_0 1077 Thus LHS = r MOD k = RHS by exp_mod_step_def 1078 If n <> 0, true by exp_mod_step_def 1079*) 1080val exp_mod_step_even = store_thm( 1081 "exp_mod_step_even", 1082 ``!k. 0 < k ==> !n. EVEN n ==> 1083 !m r. exp_mod_step m n k r = exp_mod_step ((SQ m) MOD k) (HALF n) k r``, 1084 rpt strip_tac >> 1085 Cases_on `n = 0` >- 1086 metis_tac[exp_mod_step_def, HALF_EQ_0, NOT_ZERO_LT_ZERO] >> 1087 rw[Once exp_mod_step_def]); 1088 1089(* Theorem: 0 < k ==> !n. ODD n ==> 1090 !m r. exp_mod_step m n k r = exp_mod_step ((SQ m) MOD k) (HALF n) k ((r * m) MOD k) *) 1091(* Proof: 1092 Note ODD n ==> ~EVEN n by ODD_EVEN 1093 and m <> 0 by EVEN 1094 The result follows by exp_mod_step_def 1095*) 1096val exp_mod_step_odd = store_thm( 1097 "exp_mod_step_odd", 1098 ``!k. 0 < k ==> !n. ODD n ==> 1099 !m r. exp_mod_step m n k r = exp_mod_step ((SQ m) MOD k) (HALF n) k ((r * m) MOD k)``, 1100 rpt strip_tac >> 1101 `~EVEN n /\ n <> 0` by metis_tac[ODD_EVEN, EVEN] >> 1102 rw[Once exp_mod_step_def]); 1103 1104(* Theorem: exp_mod_step m n k r = (r * m ** n) MOD k *) 1105(* Proof: 1106 If k = 0, true by exp_mod_step_def 1107 If k <> 0, 0 < k by NOT_ZERO_LT_ZERO 1108 By complete induction on n. 1109 if n = 0, 1110 exp_mod_step m 0 k r 1111 = r MOD k by exp_mod_step_0 1112 = (r * 1) MOD k by MULT_RIGHT_1 1113 = (r * m ** 0) MOD k by EXP_0 1114 If n <> 0, 1115 then HALF n < n by HALF_LESS, 0 < n 1116 If EVEN n, 1117 exp_mod_step m n k r 1118 = exp_mod_step ((SQ m) MOD k) (HALF n) k r by exp_mod_step_even 1119 = (r * ((SQ m) MOD k) ** (HALF n)) MOD k by induction hypothesis, HALF n < n 1120 = (r * (SQ m) ** (HALF n)) MOD k by EXP_MOD, MOD_TIMES2, MOD_MOD, 0 < k 1121 = (r * m ** n) MOD k by EXP_EVEN 1122 If ~EVEN n, 1123 then ODD n by EVEN_ODD 1124 exp_mod_step m n k r 1125 = exp_mod_step ((SQ m) MOD k) (HALF n) k ((r * m) MOD k) by exp_mod_step_odd 1126 = (((r * m) MOD k) * (((SQ m) MOD k) ** (HALF n))) MOD k by induction hypothesis, HALF n < n 1127 = ((r * m) * (SQ m) ** (HALF n)) MOD k by EXP_MOD, MOD_TIMES2, MOD_MOD, 0 < k 1128 = (r * (m * (SQ m) ** (HALF n))) MOD k by MULT_ASSOC 1129 = (r * m ** n) MOD k by EXP_ODD 1130*) 1131val exp_mod_step_eqn = store_thm( 1132 "exp_mod_step_eqn", 1133 ``!m n k r. exp_mod_step m n k r = (r * m ** n) MOD k``, 1134 rpt strip_tac >> 1135 Cases_on `k = 0` >- 1136 rw[Once exp_mod_step_def] >> 1137 `0 < k` by decide_tac >> 1138 qid_spec_tac `r` >> 1139 qid_spec_tac `n` >> 1140 qid_spec_tac `m` >> 1141 completeInduct_on `n` >> 1142 rpt strip_tac >> 1143 Cases_on `n = 0` >- 1144 rw[exp_mod_step_0] >> 1145 `HALF n < n` by rw[] >> 1146 Cases_on `EVEN n` >| [ 1147 `exp_mod_step m n k r = exp_mod_step ((SQ m) MOD k) (HALF n) k r` by rw[exp_mod_step_even] >> 1148 `_ = (r * ((SQ m) MOD k) ** (HALF n)) MOD k` by rw[] >> 1149 `_ = (r * (SQ m) ** (HALF n)) MOD k` by rw[Once EXP_MOD] >> 1150 `_ = (r * m ** n) MOD k` by rw[GSYM EXP_EVEN] >> 1151 rw[], 1152 `ODD n` by rw[ODD_EVEN] >> 1153 `exp_mod_step m n k r = exp_mod_step ((SQ m) MOD k) (HALF n) k (r * m) MOD k` by rw[exp_mod_step_odd] >> 1154 `_ = (((r * m) MOD k) * (((SQ m) MOD k) ** (HALF n))) MOD k` by rw[] >> 1155 `_ = ((r * m) * (SQ m) ** (HALF n)) MOD k` by rw[Once EXP_MOD] >> 1156 `_ = (r * (m * (SQ m) ** (HALF n))) MOD k` by rw[] >> 1157 `_ = (r * m ** n) MOD k` by rw[EXP_ODD] >> 1158 rw[] 1159 ]); 1160 1161(* Theorem: exp_mod_compute m 0 k = 1 MOD k *) 1162(* Proof: 1163 exp_mod_compute m 0 k 1164 = exp_mod_step m 0 k 1 by exp_mod_compute_def 1165 = 1 MOD k by exp_mod_step_0 1166*) 1167val exp_mod_compute_0 = store_thm( 1168 "exp_mod_compute_0", 1169 ``!m k. exp_mod_compute m 0 k = 1 MOD k``, 1170 rw[exp_mod_compute_def, exp_mod_step_0]); 1171 1172(* Theorem: exp_mod_compute m 1 k = m MOD k *) 1173(* Proof: 1174 exp_mod_compute m 1 k 1175 = exp_mod_step m 1 k 1 by exp_mod_compute_def 1176 = (1 * m) MOD k by exp_mod_step_1 1177 = m MOD k by MULT_LEFT_1 1178*) 1179val exp_mod_compute_1 = store_thm( 1180 "exp_mod_compute_1", 1181 ``!m k. exp_mod_compute m 1 k = m MOD k``, 1182 rw[exp_mod_compute_def, exp_mod_step_1]); 1183 1184(* Theorem: exp_mod_compute m 2 k = (SQ m) MOD k *) 1185(* Proof: 1186 exp_mod_compute m 2 k 1187 = exp_mod_step m 2 k 1 by exp_mod_compute_def 1188 = (1 * n * n) MOD k by exp_mod_step_2 1189 = (SQ m) MOD k by EXP_2 1190*) 1191val exp_mod_compute_2 = store_thm( 1192 "exp_mod_compute_2", 1193 ``!m k. exp_mod_compute m 2 k = (SQ m) MOD k``, 1194 rw[exp_mod_compute_def, exp_mod_step_2]); 1195 1196(* Theorem: 0 < k ==> !n. EVEN n ==> 1197 !m. exp_mod_compute m n k = exp_mod_compute ((SQ m) MOD k) (HALF n) k *) 1198(* Proof: 1199 exp_mod_compute m n k 1200 = exp_mod_step m n k 1 by exp_mod_compute_def 1201 = exp_mod_step ((SQ m) MOD k) (HALF n) k 1 by exp_mod_step_even, EVEN n 1202 = exp_mod_compute ((SQ m) MOD k) (HALF n) k by exp_mod_compute_def 1203*) 1204val exp_mod_compute_even = store_thm( 1205 "exp_mod_compute_even", 1206 ``!k. 0 < k ==> !n. EVEN n ==> 1207 !m. exp_mod_compute m n k = exp_mod_compute ((SQ m) MOD k) (HALF n) k``, 1208 rw[exp_mod_compute_def, exp_mod_step_even]); 1209 1210(* Theorem: 0 < k ==> !n. ODD n ==> 1211 !m. exp_mod_compute m n k = exp_mod_step ((SQ m) MOD k) (HALF n) k (m MOD k) *) 1212(* Proof: 1213 exp_mod_compute m n k 1214 = exp_mod_step m n k 1 by exp_mod_compute_def 1215 = exp_mod_step ((SQ m) MOD k) (HALF n) k (1 * m MOD k) by exp_mod_step_odd, ODD n 1216 = exp_mod_step ((SQ m) MOD k) (HALF n) k (m MOD k) by MULT_LEFT_1 1217*) 1218val exp_mod_compute_odd = store_thm( 1219 "exp_mod_compute_odd", 1220 ``!k. 0 < k ==> !n. ODD n ==> 1221 !m. exp_mod_compute m n k = exp_mod_step ((SQ m) MOD k) (HALF n) k (m MOD k)``, 1222 rw[exp_mod_compute_def, exp_mod_step_odd]); 1223 1224(* Theorem: exp_mod_compute m n k = (m ** n) MOD k *) 1225(* Proof: 1226 exp_mod_compute m n k 1227 = exp_mod_step m n k 1 by exp_mod_compute_def 1228 = (1 * m ** n) MOD k by exp_mod_step_eqn 1229 = (m ** n) MOD k by MULT_LEFT_1 1230*) 1231val exp_mod_compute_eqn = store_thm( 1232 "exp_mod_compute_eqn", 1233 ``!m n k. exp_mod_compute m n k = (m ** n) MOD k``, 1234 rw[exp_mod_compute_def, exp_mod_step_eqn]); 1235 1236(* ------------------------------------------------------------------------- *) 1237(* ROOT computation *) 1238(* ------------------------------------------------------------------------- *) 1239(* Pseudocode: 1240 1241Input: n r 1242Output: ROOT r n = r-th root of n. 1243 1244Make use of indentity: 1245n ^ (1/r) = 2 (n/ 2^r) ^(1/r) 1246 1247if n = 0 then 0 1248else (* precompute *) let x = 2 * r-th root of (n DIV (2 ** r)) 1249 (* apply *) in if n < (SUC x) ** r then x else (SUC x) 1250*) 1251 1252 1253(* Define root computation by root identity *) 1254val root_compute_def = tDefine "root_compute" ` 1255 root_compute r n = 1256 if 0 < r then 1257 if n = 0 then 0 1258 else (let x = 2 * root_compute r (n DIV (exp_compute 2 r)) in 1259 if n < exp_compute (SUC x) r then x else (SUC x)) 1260 else ROOT 0 n 1261` (WF_REL_TAC `measure (\(r, n). n)` >> 1262 rw[] >> 1263 `1 < 2 ** r` by rw[ONE_LT_EXP] >> 1264 rw[exp_compute_eqn, DIV_LESS]); 1265 1266(* 1267> EVAL ``root_compute 3 0``; --> 0 1268> EVAL ``root_compute 3 2``; --> 1 1269> EVAL ``root_compute 3 3``; --> 1 1270> EVAL ``root_compute 3 7``; --> 1 1271> EVAL ``root_compute 3 8``; --> 2 1272> EVAL ``root_compute 3 10``; --> 2 1273> EVAL ``root_compute 3 20``; --> 2 1274> EVAL ``root_compute 3 27``; --> 3 1275> EVAL ``root_compute 3 26``; --> 2 1276> EVAL ``root_compute 3 28``; --> 3 1277> EVAL ``root_compute 3 63``; --> 3 1278> EVAL ``root_compute 3 64``; --> 4 1279> EVAL ``root_compute 3 65``; --> 4 1280*) 1281 1282(* Theorem: root_compute r n = 1283 if 0 < r then if n = 0 then 0 else (let x = 2 * root_compute r (n DIV (2 ** r)) 1284 in if n < (SUC x) ** r then x else (SUC x)) 1285 else ROOT 0 n *) 1286(* Proof: by root_compute_def, exp_compute_eqn *) 1287val root_compute_alt = store_thm( 1288 "root_compute_alt", 1289 ``!r n. root_compute r n = if 0 < r then 1290 if n = 0 then 0 else 1291 (let x = 2 * root_compute r (n DIV (2 ** r)) 1292 in if n < (SUC x) ** r then x else (SUC x)) 1293 else ROOT 0 n``, 1294 rw[Once root_compute_def, exp_compute_eqn]); 1295 1296(* Theorem: root_compute r 0 = 0 *) 1297(* Proof: by root_compute_def *) 1298val root_compute_0 = store_thm( 1299 "root_compute_0", 1300 ``!r. 0 < r ==> (root_compute r 0 = 0)``, 1301 rw[Once root_compute_def]); 1302 1303(* Theorem: root_compute 1 n = n *) 1304(* Proof: 1305 By complete induction on n. 1306 Assume !m. m < n ==> (root_compute 1 m = m) 1307 To show: root_compute 1 n = n 1308 1309 Note HALF n < n by HALF_LESS, 0 < n 1310 1311 root_compute 1 n 1312 = let x = 2 * root_compute 1 (n DIV (2 ** 1)) 1313 in if n < (SUC x) ** 1 then x else (SUC x) by root_compute_alt 1314 = let x = 2 * root_compute 1 (HALF n) 1315 in if n < SUC x then x else SUC x by simplification, EXP_1 1316 = let x = 2 * HALF n 1317 in if n < SUC x then x else SUC x by induction hypothesis 1318 = if n < SUC (2 * HALF n) then (2 * HALF n) else SUC (2 * HALF n) 1319 1320 If EVEN n, 1321 Then n = 2 * HALF n by EVEN_HALF 1322 Since n < SUC n by LESS_SUC 1323 root_compute 1 n 1324 = 2 * HALF n by above 1325 = n by EVEN_HALF 1326 If ~(EVEN n), then ODD n by EVEN_ODD 1327 Then n = SUC (2 * HALF n) by ODD_HALF, ADD1 1328 or n < SUC (2 * HALF n) = F 1329 root_compute 1 n 1330 = SUC (2 * HALF n) by above 1331 = n by ODD_HALF 1332*) 1333val root_compute_1 = store_thm( 1334 "root_compute_1", 1335 ``!n. root_compute 1 n = n``, 1336 completeInduct_on `n` >> 1337 rw[Once root_compute_alt] >| [ 1338 spose_not_then strip_assume_tac >> 1339 `ODD n` by metis_tac[EVEN_HALF, ODD_EVEN] >> 1340 `n = SUC (2 * HALF n)` by metis_tac[ODD_HALF, ADD1] >> 1341 decide_tac, 1342 spose_not_then strip_assume_tac >> 1343 `EVEN n` by metis_tac[ODD_HALF, ADD1, EVEN_ODD] >> 1344 `n = 2 * HALF n` by rw[EVEN_HALF] >> 1345 decide_tac 1346 ]); 1347 1348(* Theorem: root_compute r n = ROOT r n *) 1349(* Proof: 1350 If r = 0, true by by root_compute_alt 1351 If r <> 0, then 0 < r. 1352 By complete induction on n. 1353 Assume: !m. m < n ==> (root_compute r m = ROOT r m) 1354 To show: root_compute r n = ROOT r n 1355 1356 If n = 0, 1357 root_compute r 0 1358 = 0 by root_compute_0 1359 = ROOT r 0 by ROOT_COMPUTE 1360 If n <> 0, 1361 Note 1 < 2 ** r by ONE_LT_EXP, 0 < r 1362 Thus n DIV 2 ** r < n by DIV_LESS, 0 < n 1363 This makes the induction hypothesis applicable for m = n DIV 2 ** r. 1364 Let x = 2 * ROOT r (n DIV 2 ** r). 1365 1366 root_compute r n 1367 = if n < SUC x ** r then x else SUC x by root_compute_alt, induction hypothesis 1368 = ROOT r n by ROOT_COMPUTE 1369*) 1370val root_compute_eqn = store_thm( 1371 "root_compute_eqn", 1372 ``!r n. root_compute r n = ROOT r n``, 1373 rpt strip_tac >> 1374 Cases_on `r = 0` >- 1375 rw[Once root_compute_alt] >> 1376 `0 < r` by decide_tac >> 1377 completeInduct_on `n` >> 1378 Cases_on `n = 0` >- 1379 rw[root_compute_0, ROOT_COMPUTE] >> 1380 `1 < 2 ** r` by rw[ONE_LT_EXP] >> 1381 `n DIV 2 ** r < n` by rw[DIV_LESS] >> 1382 qabbrev_tac `x = 2 * ROOT r (n DIV 2 ** r)` >> 1383 rw[Once root_compute_alt] >| [ 1384 `ROOT r n = x` by rw[Once ROOT_COMPUTE, Abbr`x`] >> 1385 decide_tac, 1386 `ROOT r n = SUC x` by rw[Once ROOT_COMPUTE, Abbr`x`] >> 1387 decide_tac 1388 ]); 1389 1390(* ------------------------------------------------------------------------- *) 1391(* Square Root Computation *) 1392(* ------------------------------------------------------------------------- *) 1393 1394(* Overload square-root *) 1395val _ = overload_on("sqrt_compute", ``root_compute 2``); 1396 1397(* Theorem: sqrt_compute n = SQRT n *) 1398(* Proof: by root_compute_eqn *) 1399val sqrt_compute_eqn = store_thm( 1400 "sqrt_compute_eqn", 1401 ``!n. sqrt_compute n = SQRT n``, 1402 rw[root_compute_eqn]); 1403 1404(* 1405> EVAL ``sqrt_compute 1``; --> 1 1406> EVAL ``sqrt_compute 2``; --> 1 1407> EVAL ``sqrt_compute 3``; --> 1 1408> EVAL ``sqrt_compute 4``; --> 2 1409> EVAL ``sqrt_compute 16``; --> 4 1410> EVAL ``sqrt_compute 121``; --> 11 1411> EVAL ``sqrt_compute 1024``; --> 32 1412*) 1413 1414(* ------------------------------------------------------------------------- *) 1415(* Power Free Test *) 1416(* ------------------------------------------------------------------------- *) 1417(* Pseudocode for power-free check. 1418 1419Input: n 1420Output: true if the only way to write n = b ** c is b = n , c = 1. 1421 1422c <- 2 1423while (c <= LOG2 n) { 1424 b <- root n c 1425 if (exp b c = n) return (b, c) 1426 c <- c + 1 1427} 1428or: 1429c <- 2 1430while (c <= ulog n) { 1431 b <- root n c 1432 if (exp b c = n) return (b, c) 1433 c <- c + 1 1434} 1435 1436The following theorems guarantee that exit must occur within the while-loop. 1437 1438perfect_power_bound_LOG2; 1439|- !n. 1 < n ==> !m. perfect_power n m <=> ?k. k <= LOG2 n /\ (n = m ** k) 1440perfect_power_bound_ulog; 1441|- !n. 0 < n ==> !m. perfect_power n m <=> ?k. k <= ulog n /\ (n = m ** k) 1442*) 1443 1444(* Define power-index computation, search downwards from k = LOG2 n *) 1445val power_index_compute_def = Define` 1446 power_index_compute n k = 1447 if k <= 1 then 1 1448 else if exp_compute (root_compute k n) k = n then k 1449 else power_index_compute n (k - 1) 1450`; 1451 1452(* 1453> EVAL ``power_index_compute 8 (log_compute 8)``; --> 3 1454> EVAL ``power_index_compute 7 (log_compute 7)``; --> 1 1455> EVAL ``power_index_compute 9 (log_compute 9)``; --> 2 1456> EVAL ``power_index_compute 10 (log_compute 10)``; --> 1 1457> EVAL ``power_index_compute 1024 (log_compute 1024)``; --> 10 1458*) 1459 1460(* Theorem: power_index_compute n k = 1461 if k <= 1 then 1 else if (ROOT k n) ** k = n then k else power_index_compute n (k - 1) *) 1462(* Proof: by power_index_compute_def, root_compute_eqn, exp_compute_eqn *) 1463val power_index_compute_alt = store_thm( 1464 "power_index_compute_alt", 1465 ``!n k. power_index_compute n k = 1466 if k <= 1 then 1 else if (ROOT k n) ** k = n then k else power_index_compute n (k - 1)``, 1467 rw[Once power_index_compute_def, root_compute_eqn, exp_compute_eqn]); 1468 1469(* Theorem: power_index_compute n k = power_index n k *) 1470(* Proof: 1471 By complete induction on k. 1472 Expand by power_index_compute_alt, this is to show: 1473 (1) k <= 1 ==> 1 = power_index n k 1474 Note k = 0 or k = 1 by k <= 1 1475 and power_index n 0 = 1 by power_index_0 1476 and power_index n 1 = 1 by power_index_1 1477 (2) ~(k <= 1) /\ ROOT k n ** k = n ==> k = power_index n k 1478 That is 1 < k by ~(k <= 1) 1479 Thus power_index n k = k by power_index_exact_root, 0 < k 1480 (3) ROOT k n ** k <> n ==> power_index_compute n (k - 1) = power_index n k 1481 power_index_compute n (k - 1) 1482 = power_index n (k - 1) by induction hypothesis, k - 1 < k 1483 = power_index n k by power_index_not_exact_root 1484*) 1485val power_index_compute_eqn = store_thm( 1486 "power_index_compute_eqn", 1487 ``!n k. power_index_compute n k = power_index n k``, 1488 completeInduct_on `k` >> 1489 rw[Once power_index_compute_alt] >- 1490 metis_tac[power_index_0, power_index_1, DECIDE``k <= 1 <=> (k = 0) \/ (k = 1)``] >- 1491 rw[power_index_exact_root] >> 1492 rw[power_index_not_exact_root]); 1493 1494(* Define power-free check *) 1495val power_free_check_def = Define` 1496 power_free_check n <=> (1 < n) /\ (power_index_compute n (ulog_compute n) = 1) 1497`; 1498 1499(* 1500> EVAL ``power_free_check 6``; --> T 1501> EVAL ``power_free_check 7``; --> T 1502> EVAL ``power_free_check 8``; --> F 1503> EVAL ``power_free_check 9``; --> F 1504> EVAL ``power_free_check 10``; --> T 1505> EVAL ``power_free_check 26``; --> T 1506> EVAL ``power_free_check 27``; --> F 1507> EVAL ``power_free_check 127``; --> T 1508> EVAL ``power_free_check 128``; --> F 1509> EVAL ``power_free_check 0``; --> F 1510> EVAL ``power_free_check 1``; --> F 1511> EVAL ``power_free_check 2``; --> T 1512*) 1513 1514(* Theorem: power_free_check 0 <=> F *) 1515(* Proof: by power_free_check_def *) 1516val power_free_check_0 = store_thm( 1517 "power_free_check_0", 1518 ``power_free_check 0 <=> F``, 1519 rw[power_free_check_def]); 1520 1521(* Theorem: power_free_check 1 <=> F *) 1522(* Proof: by power_free_check_def *) 1523val power_free_check_1 = store_thm( 1524 "power_free_check_1", 1525 ``power_free_check 1 <=> F``, 1526 rw[power_free_check_def]); 1527 1528(* Theorem: power_free_check n <=> (1 < n) /\ (power_index n (ulog n) = 1) *) 1529(* Proof: by power_free_check_def, power_index_compute_eqn, ulog_compute_eqn *) 1530val power_free_check_alt = store_thm( 1531 "power_free_check_alt", 1532 ``!n. power_free_check n <=> (1 < n) /\ (power_index n (ulog n) = 1)``, 1533 rw[power_free_check_def, power_index_compute_eqn, ulog_compute_eqn]); 1534 1535(* Theorem: power_free_check n <=> power_free n *) 1536(* Proof: by power_free_check_alt, power_free_by_power_index_ulog *) 1537val power_free_check_eqn = store_thm( 1538 "power_free_check_eqn", 1539 ``!n. power_free_check n <=> power_free n``, 1540 rw[power_free_check_alt, power_free_by_power_index_ulog]); 1541 1542(* Theorem: power_free n <=> power_free_check n *) 1543(* Proof: by power_free_check_eqn *) 1544val power_free_eval = store_thm( 1545 "power_free_eval[compute]", 1546 ``!n. power_free n <=> power_free_check n``, 1547 rw[power_free_check_eqn]); 1548 1549(* 1550> EVAL ``power_free 0``; = F 1551> EVAL ``power_free 1``; = F 1552> EVAL ``power_free 2``; = T 1553> EVAL ``power_free 3``; = T 1554> EVAL ``power_free 4``; = F 1555> EVAL ``power_free 5``; = T 1556> EVAL ``power_free 6``; = T 1557> EVAL ``power_free 7``; = T 1558> EVAL ``power_free 8``; = F 1559> EVAL ``power_free 9``; = F 1560*) 1561 1562(* Theorem alias *) 1563val power_free_check_by_ulog = save_thm("power_free_check_by_ulog", power_free_check_alt); 1564(* val power_free_check_by_ulog = 1565 |- !n. power_free_check n <=> 1 < n /\ (power_index n (ulog n) = 1): thm *) 1566 1567(* Theorem: power_free_check n <=> (1 < n) /\ (power_index n (LOG2 n) = 1) *) 1568(* Proof: by power_free_check_eqn, power_free_by_power_index_LOG2 *) 1569val power_free_check_by_LOG2 = store_thm( 1570 "power_free_check_by_LOG2", 1571 ``!n. power_free_check n <=> (1 < n) /\ (power_index n (LOG2 n) = 1)``, 1572 rw[power_free_check_eqn, power_free_by_power_index_LOG2]); 1573 1574(* ------------------------------------------------------------------------- *) 1575(* GCD computation *) 1576(* ------------------------------------------------------------------------- *) 1577(* Pseudocode: 1578 1579Input: n m 1580Output: gcd n m 1581 1582Recursive version: 1583 1584if (n = m) return n // gcd n n = n 1585if (n = 0) return m // gcd 0 m = m 1586if (m = 0) return n // gcd n 0 = n 1587 1588if ~(n < m) return gcd m n // gcd n m = gcd m n, ensure first < second. 1589 1590if (EVEN n) then 1591 if (EVEN m) then return 2 * (gcd (HALF n) (HALF m)) // because 2 is a common divisor. 1592 else return gcd (HALF n) m // because 2 is not a common divisor. 1593else 1594 if (EVEN m) then return gcd n (HALF m) // because 2 is not a common divisor. 1595 else return gcd n (HALF (m - n)) // because (m - n) is EVEN, as both are ODD. 1596 1597Iterative version: 1598 1599if (n = 0) return m // gcd 0 m = m 1600if (m = 0) return n // gcd n 0 = n 1601 1602// Let shift k := log K, where K is the greatest power of 2 dividing both n and m. 1603k <- 0 // initialize shift 1604while (EVEN n) and (EVEN m) { 1605 k <- k + 1 // increment shift 1606 n <- HALF n 1607 m <- HALF m 1608} 1609 1610// remove all factors of 2 in n -- they are not common. 1611 while (EVEN n) n <- HALF n 1612 1613// From here on, n is always odd. 1614do { 1615 // remove all factors of 2 in m -- they are not common. 1616 // note: m is not zero, so while will terminate. 1617 while (EVEN m) m <- HALF m 1618 1619 // Now n and m are both odd. Swap if necessary so n <= m, 1620 // then set m = m - n (which is even). 1621 if (n > m) swap (n, m) // ensure n <= m 1622 m <- m - n // set EVEN m 1623} while (m != 0) 1624 1625return n * (2 ** k) // restore common factors of 2 1626 1627*) 1628 1629(* Proof in gcd theory: 1630 1631val BINARY_GCD = store_thm("BINARY_GCD", 1632 ``!m n. 1633 (EVEN m /\ EVEN n ==> (gcd m n = 2 * gcd (m DIV 2) (n DIV 2))) /\ 1634 (EVEN m /\ ODD n ==> (gcd m n = gcd (m DIV 2) n))``, 1635 SIMP_TAC bool_ss [EVEN_EXISTS] THEN REVERSE (REPEAT STRIP_TAC) 1636 THEN `0 < 2` by (MATCH_MP_TAC PRIME_POS THEN REWRITE_TAC [PRIME_2]) 1637 THEN FULL_SIMP_TAC bool_ss [GCD_COMMON_FACTOR, 1638 ONCE_REWRITE_RULE [MULT_COMM] MULT_DIV, 1639 ONCE_REWRITE_RULE [GCD_SYM] ODD_IMP_GCD_CANCEL_EVEN]); 1640*) 1641 1642(* Rework the proof of BINARY_GCD in gcd theory. *) 1643 1644(* Theorem: ODD n ==> (gcd n (2 * m) = gcd n m) *) 1645(* Proof: 1646 Since ODD n, 2 cannot be the common divisor of n and (2 * m). 1647 Therefore, any common divisor d of n and (2 * m) divides n and m, 1648 or gcd n m = gcd n (2 * m) 1649 1650 Let a = gcd n (2 * m), b = gcd n m. 1651 1652 Step 1: show !c. c divides n /\ c divides (2 * m) ==> c divides m 1653 Claim coprime 2 c. 1654 Proof: Let d = gcd 2 c. 1655 Then d divides 2 /\ d divides c by GCD_IS_GREATEST_COMMON_DIVISOR 1656 so d <= 2 by DIVIDES_LE, 0 < 2 1657 and d <> 0` by ZERO_DIVIDES, 2 <> 0 1658 Note ODD n ==> ODD c by DIVIDES_ODD, c divides n 1659 Thus d <> 2 by DIVIDES_MOD_0, ODD_MOD2 1660 or d = 1 by arithmetic 1661 Then c divides m by L_EUCLIDES, coprime 2 c /\ c divides (2 * m). 1662 1663 Step 2: show a divides b 1664 Note a divides n /\ a divides (2 * m) by GCD_IS_GREATEST_COMMON_DIVISOR 1665 or a divides n /\ a divides m by Step 1 1666 so a divides gcd n m = b by GCD_IS_GREATEST_COMMON_DIVISOR 1667 1668 Step 3: show b divides a 1669 Note b divides n /\ b divides m by GCD_IS_GREATEST_COMMON_DIVISOR 1670 or b divides n /\ b divides (2 * m) by DIVIDES_MULTIPLE 1671 or b divides gcd n (2 * m) = a by GCD_IS_GREATEST_COMMON_DIVISOR 1672 1673 Since a divides b /\ b divides a, a = b by DIVIDES_ANTISYM 1674*) 1675val ODD_IMP_GCD_CANCEL_EVEN = store_thm( 1676 "ODD_IMP_GCD_CANCEL_EVEN", 1677 ``!m n. ODD n ==> (gcd n (2 * m) = gcd n m)``, 1678 rpt strip_tac >> 1679 qabbrev_tac `a = gcd n (2 * m)` >> 1680 qabbrev_tac `b = gcd n m` >> 1681 `!c. c divides n /\ c divides (2 * m) ==> c divides m` by 1682 (rpt strip_tac >> 1683 `coprime 2 c` by 1684 (qabbrev_tac `d = gcd 2 c` >> 1685 `d divides 2 /\ d divides c` by rw[GCD_IS_GREATEST_COMMON_DIVISOR, Abbr`d`] >> 1686 `d <= 2` by rw[DIVIDES_LE] >> 1687 `d <> 0` by metis_tac[ZERO_DIVIDES, DECIDE``2 <> 0``] >> 1688 `d <> 2` by metis_tac[DIVIDES_ODD, DIVIDES_MOD_0, ODD_MOD2, DECIDE``1 <> 0 /\ 0 < 2``] >> 1689 decide_tac) >> 1690 metis_tac[L_EUCLIDES]) >> 1691 metis_tac[GCD_IS_GREATEST_COMMON_DIVISOR, DIVIDES_MULTIPLE, DIVIDES_ANTISYM]); 1692 1693(* Proof in gcd theory: 1694 1695val ODD_IMP_GCD_CANCEL_EVEN = prove( 1696 ``!n. ODD n ==> (gcd n (2 * m) = gcd n m)``, 1697 REPEAT STRIP_TAC 1698 THEN MATCH_MP_TAC GCD_CANCEL_MULT 1699 THEN ONCE_REWRITE_TAC [GCD_SYM] 1700 THEN REVERSE (`~divides 2 n` by ALL_TAC) 1701 THEN1 (MP_TAC (Q.SPEC `n` (MATCH_MP PRIME_GCD PRIME_2)) 1702 THEN ASM_REWRITE_TAC []) 1703 THEN REWRITE_TAC [divides_def] 1704 THEN ONCE_REWRITE_TAC [MULT_COMM] 1705 THEN REWRITE_TAC [GSYM EVEN_EXISTS] 1706 THEN FULL_SIMP_TAC bool_ss [ODD_EVEN]); 1707*) 1708 1709(* Theorem: (EVEN m /\ EVEN n ==> (gcd m n = 2 * gcd (HALF m) (HALF n))) /\ 1710 (EVEN m /\ ODD n ==> (gcd m n = gcd (HALF m) n)) *) 1711(* Proof: 1712 This is to show: 1713 (1) EVEN m /\ EVEN n ==> (gcd m n = 2 * gcd (HALF m) (HALF n)) 1714 Let hm = HALF m, then m = 2 * hm by EVEN_HALF 1715 Let hn = HALF n, then n = 2 * hn by EVEN_HALF 1716 gcd m n 1717 = gcd (2 * hm) (2 * hn) by above 1718 = 2 * gcd hm hn by GCD_COMMON_FACTOR 1719 (2) EVEN m /\ ODD n ==> (gcd m n = gcd (HALF m) n) 1720 Let hm = HALF m, then m = 2 * hm by EVEN_HALF 1721 gcd m n 1722 = gcd (2 * hm) n by above 1723 = gcd n (2 * hm) by GCD_SYM 1724 = gcd n hm by ODD_IMP_GCD_CANCEL_EVEN 1725 = gcd hm n by GCD_SYM 1726*) 1727val BINARY_GCD = store_thm( 1728 "BINARY_GCD", 1729 ``!m n. (EVEN m /\ EVEN n ==> (gcd m n = 2 * gcd (HALF m) (HALF n))) /\ 1730 (EVEN m /\ ODD n ==> (gcd m n = gcd (HALF m) n))``, 1731 rpt strip_tac >- 1732 metis_tac[EVEN_HALF, GCD_COMMON_FACTOR] >> 1733 metis_tac[EVEN_HALF, ODD_IMP_GCD_CANCEL_EVEN, GCD_SYM]); 1734 1735(* Note: For ODD m /\ ODD n, 1736 Let n < m, then apply GCD_SUB_L: gcd m n = gcd (m - n) n, with EVEN (m - n) 1737*) 1738 1739(* Convert BINARY_GCD into an algorithm *) 1740val gcd_compute_def = Define` 1741 gcd_compute n m = 1742 if n = 0 then m 1743 else if m = 0 then n 1744 else if n = m then n 1745 else (* n <> m *) 1746 if EVEN n then if EVEN m then 2 * gcd_compute (HALF n) (HALF m) 1747 else (* ODD m *) gcd_compute (HALF n) m 1748 else (* ODD n *) if EVEN m then gcd_compute n (HALF m) 1749 else (* ODD m *) 1750 if n < m then gcd_compute n (m - n) else gcd_compute (n - m) m 1751`; 1752(* Techniques for recursive definition: 1753(1) Keep the parameter order for recursive calls (no swapping of parameter). 1754(2) Provide exit conditions for small enough parameters. 1755(3) If possible, keep the parameters decreasing: 1756 with luck, termination proof is automatic! 1757*) 1758 1759(* 1760> EVAL ``gcd_compute 6 9``; --> 3 1761> EVAL ``gcd_compute 6 8``; --> 2 1762> EVAL ``gcd_compute 6 7``; --> 1 1763> EVAL ``gcd_compute 6 6``; --> 6 1764> EVAL ``gcd_compute 6 5``; --> 1 1765> EVAL ``gcd_compute 6 4``; --> 2 1766> EVAL ``gcd_compute 6 3``; --> 3 1767> EVAL ``gcd_compute 6 2``; --> 2 1768> EVAL ``gcd_compute 6 1``; --> 1 1769> EVAL ``gcd_compute 6 0``; --> 6 1770*) 1771 1772(* Theorem: (gcd_compute 0 n = n) /\ (gcd_compute n 0 = n) *) 1773(* Proof: by gcd_compute_def *) 1774val gcd_compute_0 = store_thm( 1775 "gcd_compute_0", 1776 ``!n. (gcd_compute 0 n = n) /\ (gcd_compute n 0 = n)``, 1777 rw[Once gcd_compute_def] >> 1778 rw[Once gcd_compute_def]); 1779 1780(* Theorem: gcd_compute n n = n *) 1781(* Proof: by gcd_compute_def *) 1782val gcd_compute_id = store_thm( 1783 "gcd_compute_id", 1784 ``!n. gcd_compute n n = n``, 1785 rw[Once gcd_compute_def]); 1786 1787(* Theorem: EVEN n /\ EVEN m ==> (gcd_compute n m = 2 * gcd_compute (HALF n) (HALF m)) *) 1788(* Proof: 1789 If n = m, 1790 Then HALF n = HALF m, 1791 so gcd_compute n m = n by gcd_compute_id 1792 and gcd_compute (HALF n) (HALF m) = HALF n by gcd_compute_id 1793 Thus result is true by EVEN_HALF 1794 If n <> m, 1795 Note EVEN 0 by EVEN 1796 and HALF 0 = 0 by ZERO_DIV 1797 If n = 0, 1798 gcd_compute 0 m = m by gcd_compute_0 1799 gcd_compute (HALF 0) (HALF m) = HALF m by gcd_compute_0 1800 Thus result is true by EVEN_HALF, EVEN m 1801 If m = 0, 1802 gcd_compute n 0 = n by gcd_compute_0 1803 gcd_compute (HALF n) (HALF 0) = HALF n by gcd_compute_0 1804 Thus result is true by EVEN_HALF, EVEN n 1805 If n <> 0 and m <> 0, 1806 The result is true by gcd_compute_def 1807*) 1808val gcd_compute_even_even = store_thm( 1809 "gcd_compute_even_even", 1810 ``!m n. EVEN n /\ EVEN m ==> (gcd_compute n m = 2 * gcd_compute (HALF n) (HALF m))``, 1811 rpt strip_tac >> 1812 Cases_on `n = m` >- 1813 metis_tac[gcd_compute_id, EVEN_HALF] >> 1814 `EVEN 0 /\ (HALF 0 = 0)` by rw[] >> 1815 Cases_on `n = 0` >- 1816 metis_tac[gcd_compute_0, EVEN_HALF] >> 1817 Cases_on `m = 0` >- 1818 metis_tac[gcd_compute_0, EVEN_HALF] >> 1819 rw[Once gcd_compute_def]); 1820 1821(* Theorem: EVEN n /\ ODD m ==> (gcd_compute n m = gcd_compute (HALF n) m) *) 1822(* Proof: 1823 If n = 0, 1824 Then HALF 0 = 0 by ZERO_DIV, 0 < 2 1825 so gcd_compute 0 m = m by gcd_compute_0 1826 and gcd_compute (HALF 0) m = m by gcd_compute_0 1827 Thus result is true. 1828 If n <> 0, 1829 Note m <> 0 by ODD 1830 and ~EVEN m by EVEN_ODD 1831 so n <> m by EVEN n, ~EVEN m 1832 Thus result is true by gcd_compute_def 1833*) 1834val gcd_compute_even_odd = store_thm( 1835 "gcd_compute_even_odd", 1836 ``!m n. EVEN n /\ ODD m ==> (gcd_compute n m = gcd_compute (HALF n) m)``, 1837 rpt strip_tac >> 1838 Cases_on `n = 0` >- 1839 rw[gcd_compute_0] >> 1840 `~EVEN m` by rw[GSYM ODD_EVEN] >> 1841 `n <> m` by metis_tac[] >> 1842 `m <> 0` by metis_tac[ODD] >> 1843 rw[Once gcd_compute_def]); 1844 1845(* Theorem: ODD n /\ EVEN m ==> (gcd_compute n m = gcd_compute n (HALF m)) *) 1846(* Proof: 1847 If m = 0, 1848 Then HALF 0 = 0 by ZERO_DIV, 0 < 2 1849 so gcd_compute n 0 = n by gcd_compute_0 1850 and gcd_compute n (HALF 0) = n by gcd_compute_0 1851 Thus result is true. 1852 If m <> 0, 1853 Note n <> 0 by ODD 1854 and ~EVEN n by EVEN_ODD 1855 so n <> m by EVEN m, ~EVEN n 1856 Thus result is true by gcd_compute_def 1857*) 1858val gcd_compute_odd_even = store_thm( 1859 "gcd_compute_odd_even", 1860 ``!m n. ODD n /\ EVEN m ==> (gcd_compute n m = gcd_compute n (HALF m))``, 1861 rpt strip_tac >> 1862 Cases_on `m = 0` >- 1863 rw[gcd_compute_0] >> 1864 `~EVEN n` by rw[GSYM ODD_EVEN] >> 1865 `n <> m` by metis_tac[] >> 1866 `n <> 0` by metis_tac[ODD] >> 1867 rw[Once gcd_compute_def]); 1868 1869(* Theorem: ODD n /\ ODD m ==> (gcd_compute n m = 1870 if n < m then gcd_compute n (m - n) else gcd_compute (n - m) m) *) 1871(* Proof: 1872 If n = m, 1873 Then gcd_compute n n by gcd_compute_id 1874 and gcd_compute (n - n) n 1875 = gcd_compute 0 n by n - n = 0 1876 = n by gcd_compute_0 1877 Hence true. 1878 If n <> m, 1879 Note ~EVEN n /\ ~EVEN m by ODD_EVEN 1880 also n <> 0 /\ m <> 0 by ODD 1881 The result is true by gcd_compute_def 1882*) 1883val gcd_compute_odd_odd = store_thm( 1884 "gcd_compute_odd_odd", 1885 ``!m n. ODD n /\ ODD m ==> (gcd_compute n m = 1886 if n < m then gcd_compute n (m - n) else gcd_compute (n - m) m)``, 1887 rpt strip_tac >> 1888 Cases_on `n = m` >- 1889 metis_tac[gcd_compute_id, gcd_compute_0, DECIDE``n - n = 0``] >> 1890 `~EVEN n /\ ~EVEN m` by rw[GSYM ODD_EVEN] >> 1891 `n <> 0 /\ m <> 0` by metis_tac[ODD] >> 1892 rw[Once gcd_compute_def]); 1893 1894(* Theorem: gcd_compute n m = gcd n m *) 1895(* Proof: 1896 By complete induction on (n + m). 1897 Apply gcd_compute_def, this is to show: 1898 (1) m = gcd 0 m, true by GCD_0 1899 (2) n <> 0 ==> n = gcd n 0, true by GCD_0 1900 (3) m <> 0 ==> m = gcd m m by GCD_REF 1901 (4) n <> 0 /\ m <> 0 /\ n <> m /\ EVEN n /\ EVEN m ==> 2 * gcd_compute (HALF n) (HALF m) = gcd n m 1902 Note HALF n < n /\ HALF m < m by HALF_LESS, 0 < n, 0 < m 1903 so HALF n + HALF m < n + m by arithmetic 1904 2 * gcd_compute (HALF n) (HALF m) 1905 = 2 * gcd (HALF n) (HALF m) by induction hypothesis 1906 = gcd n m by BINARY_GCD 1907 (5) n <> 0 /\ m <> 0 /\ n <> m /\ EVEN n /\ ~EVEN m ==> gcd_compute (HALF n) m = gcd n m 1908 Note HALF n < n by HALF_LESS, 0 < n 1909 so HALF n + m < n + m by arithmetic 1910 Now ODD m by EVEN_ODD 1911 gcd_compute (HALF n) m 1912 = gcd (HALF n) m by induction hypothesis 1913 = gcd n m by BINARY_GCD 1914 (6) n <> 0 /\ m <> 0 /\ n <> m /\ ~EVEN n /\ EVEN m ==> gcd_compute n (HALF m) = gcd n m 1915 Note ODD n by EVEN_ODD 1916 and HALF m < m by HALF_LESS, 0 < m 1917 so n + HALF m < n + m by arithmetic 1918 gcd_compute n (HALF m) 1919 = gcd n (HALF m) by induction hypothesis 1920 = gcd (HALF m) n by GCD_SYM 1921 = gcd m n by BINARY_GCD 1922 = gcd n m by GCD_SYM 1923 (7) n <> 0 /\ m <> 0 /\ n <> m /\ ~EVEN n /\ ~EVEN m /\ n < m ==> gcd_compute n (m - n) = gcd n m 1924 Note n + (m - n) = m by n < m 1925 so n + (m - n) < n + m by 0 < n 1926 gcd_compute n (m - n) 1927 = gcd n (m - n) by induction hypothesis 1928 = gcd n m by GCD_SUB_R, n < m ==> n <= m 1929 (8) n <> 0 /\ m <> 0 /\ n <> m /\ ~EVEN n /\ ~EVEN m /\ ~(n < m) ==> gcd_compute (n - m) m = gcd n m 1930 Note m <= n by NOT_LESS, ~(n < m) 1931 and (n - m) + m = n by m <= n 1932 so (n - m) + m < n + m by 0 < m 1933 gcd_compute (n - m) m 1934 = gcd (n - m) m by induction hypothesis 1935 = gcd n m by GCD_SUB_L, m <= n 1936*) 1937val gcd_compute_eqn = store_thm( 1938 "gcd_compute_eqn", 1939 ``!m n. gcd_compute n m = gcd n m``, 1940 rpt strip_tac >> 1941 completeInduct_on `n + m` >> 1942 rpt strip_tac >> 1943 rw[Once gcd_compute_def] >| [ 1944 `HALF n < n /\ HALF m < m` by rw[HALF_LESS] >> 1945 `HALF n + HALF m < n + m` by rw[] >> 1946 rw[BINARY_GCD], 1947 metis_tac[BINARY_GCD, EVEN_ODD], 1948 metis_tac[BINARY_GCD, EVEN_ODD, GCD_SYM], 1949 metis_tac[GCD_SUB_R, LESS_IMP_LESS_OR_EQ], 1950 metis_tac[GCD_SUB_L, NOT_LESS] 1951 ]); 1952 1953(* Theorem: (gcd_compute 1 n = 1) /\ (gcd_compute n 1 = 1) *) 1954(* Proof: by gcd_compute_eqn, GCD_1 *) 1955val gcd_compute_1 = store_thm( 1956 "gcd_compute_1", 1957 ``!n. (gcd_compute 1 n = 1) /\ (gcd_compute n 1 = 1)``, 1958 rw_tac std_ss[gcd_compute_eqn, GCD_1]); 1959 1960(* Theorem: gcd_compute n m = gcd_compute m n *) 1961(* Proof: gcd_compute_eqn, GCD_SYM *) 1962val gcd_compute_sym = store_thm( 1963 "gcd_compute_sym", 1964 ``!m n. gcd_compute n m = gcd_compute m n``, 1965 rw[gcd_compute_eqn, GCD_SYM]); 1966 1967(* ------------------------------------------------------------------------- *) 1968(* Phi Computation *) 1969(* ------------------------------------------------------------------------- *) 1970(* Pseudocode: 1971 1972Input: n 1973Output: phi n 1974 1975if (n = 0) return 0 // phi 0 = 0 1976j <- n // initial index 1977c <- 1 // initial count, when j = 1, always coprime j n 1978while (1 < j) { 1979 if (coprime j n) c <- c + 1 // increment count if coprime j n 1980 j <- j - 1 // decrement j 1981} 1982return c // the result 1983 1984*) 1985 1986(* Count the number of coprimes down to 1, accumulator version: 1987val count_coprime_def = Define ` 1988 count_coprime n j k = (* j = n downto 1, k = count from 0 *) 1989 if 1 < j then (* decrement j, and update count k depending on coprime j n *) 1990 count_coprime n (j - 1) (if coprime j n then k + 1 else k) 1991 else k + 1 (* add the last n = 1, always coprime *) 1992`; 1993*) 1994 1995(* Count the number of coprimes down to 1, stack version 1996val count_coprime_def = Define ` 1997 count_coprime n j = (* j = n downto 1, return count from 1 *) 1998 if 1 < j then (* decrement j, and update count k depending on coprime j n *) 1999 count_coprime n (j - 1) + (if (gcd_compute j n = 1) then 1 else 0) 2000 else if (j = 1) then 1 else 0 (* add the last n = 1, always coprime *) 2001`; 2002*) 2003 2004(* Count the number of coprimes, linear stack version *) 2005val count_coprime_def = Define ` 2006 count_coprime n j = 2007 if j = 0 then 0 2008 else if j = 1 then 1 2009 else count_coprime n (j - 1) + (if (gcd_compute j n = 1) then 1 else 0) 2010`; 2011 2012(* Compute phi function *) 2013val phi_compute_def = Define` 2014 phi_compute n = count_coprime n n 2015`; 2016 2017(* 2018> EVAL ``phi_compute 0``; --> 0 2019> EVAL ``phi_compute 1``; --> 1 2020> EVAL ``phi_compute 2``; --> 1 2021> EVAL ``phi_compute 3``; --> 2 2022> EVAL ``phi_compute 4``; --> 2 2023> EVAL ``phi_compute 5``; --> 4 2024> EVAL ``phi_compute 6``; --> 2 2025> EVAL ``phi_compute 7``; --> 6 2026> EVAL ``phi_compute 8``; --> 4 2027> EVAL ``phi_compute 9``; --> 6 2028> EVAL ``phi_compute 10``; --> 4 2029*) 2030 2031(* Theorem: phi_compute 0 = 0 *) 2032(* Proof: 2033 phi_compute 0 2034 = count_coprime 0 0 by phi_compute_def 2035 = 0 by count_coprime_def 2036*) 2037val phi_compute_0 = store_thm( 2038 "phi_compute_0", 2039 ``phi_compute 0 = 0``, 2040 rw[phi_compute_def, Once count_coprime_def]); 2041 2042(* Theorem: phi_compute 1 = 1 *) 2043(* Proof: 2044 phi_compute 1 2045 = count_coprime 1 1 by phi_compute_def 2046 = 1 by count_coprime_def 2047*) 2048val phi_compute_1 = store_thm( 2049 "phi_compute_1", 2050 ``phi_compute 1 = 1``, 2051 rw[phi_compute_def, Once count_coprime_def]); 2052 2053(* Theorem: count_coprime n 0 = 0 *) 2054(* Proof: by count_coprime_def *) 2055val count_coprime_0 = store_thm( 2056 "count_coprime_0", 2057 ``!n. count_coprime n 0 = 0``, 2058 rw[Once count_coprime_def]); 2059 2060(* Theorem: count_coprime n 1 = 1 *) 2061(* Proof: by count_coprime_def *) 2062val count_coprime_1 = store_thm( 2063 "count_coprime_1", 2064 ``!n. count_coprime n 1 = 1``, 2065 rw[Once count_coprime_def]); 2066 2067(* Theorem: count_coprime n j = if j = 0 then 0 else if j = 1 then 1 else 2068 count_coprime n (j - 1) + if coprime j n then 1 else 0 *) 2069(* Proof: by count_coprime_def, gcd_compute_eqn *) 2070val count_coprime_alt = store_thm( 2071 "count_coprime_alt", 2072 ``!n j. count_coprime n j = if j = 0 then 0 else if j = 1 then 1 else 2073 count_coprime n (j - 1) + if coprime j n then 1 else 0``, 2074 rw[Once count_coprime_def, gcd_compute_eqn]); 2075 2076(* Theorem: count_coprime n (SUC k) = count_coprime n k + (if coprime (SUC k) n then 1 else 0) *) 2077(* Proof: 2078 Apply count_coprime_def, this is to show: 2079 (1) coprime (SUC 0) n ==> 1 = count_coprime n 0 + 1, true by count_coprime_0 2080 (2) gcd (SUC 0) n <> 1 ==> 1 = count_coprime n 0, true by GCD_1 2081*) 2082val count_coprime_suc = store_thm( 2083 "count_coprime_suc", 2084 ``!n k. count_coprime n (SUC k) = count_coprime n k + (if coprime (SUC k) n then 1 else 0)``, 2085 rw[Once count_coprime_alt] >- 2086 rw[count_coprime_0] >> 2087 fs[GCD_1]); 2088 2089(* Theorem: k <= n ==> (count_coprime n k = CARD ((coprimes n) INTER (natural k))) *) 2090(* Proof: 2091 By induction on k. 2092 Base: count_coprime n 0 = CARD (coprimes n INTER natural 0) 2093 LHS = count_coprime n 0 = 0 by count_coprime_0 2094 RHS = CARD (coprimes n INTER natural 0) 2095 = CARD (coprimes n INTER {}) by natural_0 2096 = CARD {} by INTER_EMPTY 2097 = 0 = LHS by CARD_EMPTY 2098 Step: count_coprime n k = CARD (coprimes n INTER natural k) ==> 2099 count_coprime n (SUC k) = CARD (coprimes n INTER natural (SUC k)) 2100 2101 If coprime (SUC k) n, then (SUC k) IN (coprimes n) by coprimes_element 2102 But (SUC k) NOTIN (natural k) by natural_element 2103 Thus (SUC k) NOTIN (natural k) INTER (coprimes n) by IN_INTER 2104 Note FINITE (natural k) by natural_finite 2105 and FINITE (coprimes n) by coprimes_finite 2106 ==> FINITE (natural k) INTER (coprimes n) by FINITE_INTER 2107 2108 CARD (coprimes n INTER natural (SUC k)) 2109 = CARD (natural (SUC k) INTER coprimes n) by INTER_COMM 2110 = CARD (((SUC k) INSERT (natural k)) INTER coprimes n) by natural_suc 2111 = CARD ((SUC k) INSERT ((natural k) INTER (coprimes n))) by INSERT_INTER, (SUC k) IN (coprimes n) 2112 = SUC (CARD ((natural k) INTER (coprimes n))) by CARD_INSERT, (SUC k) NOTIN intersection 2113 = CARD ((natural k) INTER (coprimes n)) + 1 by ADD1 2114 = CARD ((coprimes n) INTER (natural k)) + 1 by INTER_COMM 2115 = count_coprime n k + 1 by induction hypothesis 2116 = count_coprime n k + (if coprime (SUC k) n then 1 else 0) 2117 = count_coprime n (SUC k) by count_coprime_suc 2118 2119 If ~(coprime (SUC k) n), then (SUC k) NOTIN (coprimes n) by coprimes_element 2120 CARD (coprimes n INTER natural (SUC k)) 2121 = CARD (natural (SUC k) INTER coprimes n) by INTER_COMM 2122 = CARD (((SUC k) INSERT (natural k)) INTER coprimes n) by natural_suc 2123 = CARD ((natural k) INTER (coprimes n)) by INSERT_INTER, (SUC k) NOTIN (coprimes n) 2124 = CARD (natural k) INTER (coprimes n) + 0 by ADD_0 2125 = CARD ((coprimes n) INTER (natural k)) + 0 by INTER_COMM 2126 = count_coprime n k + 0 by induction hypothesis 2127 = count_coprime n k + (if coprime (SUC k) n then 1 else 0) 2128 = count_coprime n (SUC k) by count_coprime_suc 2129 2130*) 2131val count_coprime_eqn = store_thm( 2132 "count_coprime_eqn", 2133 ``!n k. k <= n ==> (count_coprime n k = CARD ((coprimes n) INTER (natural k)))``, 2134 rpt strip_tac >> 2135 Induct_on `k` >- 2136 rw[count_coprime_0] >> 2137 rw[count_coprime_suc] >| [ 2138 `(SUC k) IN (coprimes n)` by rw[coprimes_element] >> 2139 `(SUC k) NOTIN (natural k)` by rw[natural_element] >> 2140 `(SUC k) NOTIN (natural k) INTER (coprimes n)` by rw[] >> 2141 `CARD (coprimes n INTER natural (SUC k)) = CARD (natural (SUC k) INTER coprimes n)` by rw[INTER_COMM] >> 2142 `_ = CARD (((SUC k) INSERT (natural k)) INTER coprimes n)` by rw[natural_suc] >> 2143 `_ = CARD ((SUC k) INSERT ((natural k) INTER (coprimes n)))` by rw[INSERT_INTER] >> 2144 `_ = SUC (CARD ((natural k) INTER (coprimes n)))` by rw[CARD_INSERT, natural_finite, coprimes_finite] >> 2145 `_ = CARD ((natural k) INTER (coprimes n)) + 1` by rw[] >> 2146 rw[INTER_COMM], 2147 `(SUC k) NOTIN (coprimes n)` by rw[coprimes_element] >> 2148 `CARD (coprimes n INTER natural (SUC k)) = CARD (natural (SUC k) INTER coprimes n)` by rw[INTER_COMM] >> 2149 `_ = CARD (((SUC k) INSERT (natural k)) INTER coprimes n)` by rw[natural_suc] >> 2150 `_ = CARD ((natural k) INTER (coprimes n))` by rw[INSERT_INTER] >> 2151 rw[INTER_COMM] 2152 ]); 2153 2154(* Theorem: phi_compute n = phi n *) 2155(* Proof: 2156 Note (coprimes n) INTER (natural n) by coprimes_subset 2157 and n <= n by LESS_EQ_REFL 2158 phi_compute n 2159 = count_coprime n n by phi_compute_def 2160 = CARD ((coprimes n) INTER (natural n)) by count_coprime_eqn, n <= n. 2161 = CARD (coprimes n) by SUBSET_INTER_ABSORPTION 2162 = phi n by phi_def 2163*) 2164val phi_compute_eqn = store_thm( 2165 "phi_compute_eqn", 2166 ``!n. phi_compute n = phi n``, 2167 metis_tac[phi_compute_def, phi_def, 2168 count_coprime_eqn, coprimes_subset, SUBSET_INTER_ABSORPTION, LESS_EQ_REFL]); 2169 2170(* ------------------------------------------------------------------------- *) 2171 2172(* export theory at end *) 2173val _ = export_theory(); 2174 2175(*===========================================================================*) 2176