1(* ------------------------------------------------------------------------- *) 2(* Primality Tests *) 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 "primes"; 12 13(* ------------------------------------------------------------------------- *) 14 15 16 17(* open dependent theories *) 18(* val _ = load "logPowerTheory"; *) 19open logPowerTheory; (* for SQRT *) 20open helperNumTheory helperFunctionTheory; 21open arithmeticTheory dividesTheory; 22 23 24(* ------------------------------------------------------------------------- *) 25(* Primality Tests Documentation *) 26(* ------------------------------------------------------------------------- *) 27(* Overloading: 28*) 29(* 30 31 Two Factors Properties: 32 two_factors_property_1 |- !n a b. (n = a * b) /\ a < SQRT n ==> SQRT n <= b 33 two_factors_property_2 |- !n a b. (n = a * b) /\ SQRT n < a ==> b <= SQRT n 34 two_factors_property |- !n a b. (n = a * b) ==> a <= SQRT n \/ b <= SQRT n 35 36 Primality or Compositeness based on SQRT: 37 prime_by_sqrt_factors |- !p. prime p <=> 38 1 < p /\ !q. 1 < q /\ q <= SQRT p ==> ~(q divides p) 39 prime_factor_estimate |- !n. 1 < n ==> 40 (~prime n <=> ?p. prime p /\ p divides n /\ p <= SQRT n) 41 42 Primality Testing Algorithm: 43 factor_seek_def |- !q n c. factor_seek n c q = 44 if c <= q then n 45 else if 1 < q /\ (n MOD q = 0) then q 46 else factor_seek n c (q + 1) 47 prime_test_def |- !n. prime_test n <=> 48 if n <= 1 then F else factor_seek n (1 + SQRT n) 2 = n 49 factor_seek_bound |- !n c q. 0 < n ==> factor_seek n c q <= n 50 factor_seek_thm |- !n c q. 1 < q /\ q <= c /\ c <= n ==> 51 (factor_seek n c q = n <=> !p. q <= p /\ p < c ==> ~(p divides n)) 52 prime_test_thm |- !n. prime n <=> prime_test n 53 54*) 55 56(* ------------------------------------------------------------------------- *) 57(* Helper Theorems *) 58(* ------------------------------------------------------------------------- *) 59 60(* ------------------------------------------------------------------------- *) 61(* Two Factors Properties *) 62(* ------------------------------------------------------------------------- *) 63 64(* Theorem: (n = a * b) /\ a < SQRT n ==> SQRT n <= b *) 65(* Proof: 66 If n = 0, then a = 0 or b = 0 by MULT_EQ_0 67 But SQRT 0 = 0 by SQRT_0 68 so a <> 0, making b = 0, and SQRT n <= b true. 69 If n <> 0, a <> 0 and b <> 0 by MULT_EQ_0 70 By contradiction, suppose b < SQRT n. 71 Then n = a * b < a * SQRT n by LT_MULT_LCANCEL, 0 < a. 72 and a * SQRT n < SQRT n * SQRT n by LT_MULT_RCANCEL, 0 < SQRT n. 73 making n < (SQRT n) ** 2 by LESS_TRANS, EXP_2 74 This contradicts (SQRT n) ** 2 <= n by SQRT_PROPERTY 75*) 76val two_factors_property_1 = store_thm( 77 "two_factors_property_1", 78 ``!n a b. (n = a * b) /\ a < SQRT n ==> SQRT n <= b``, 79 rpt strip_tac >> 80 Cases_on `n = 0` >| [ 81 `a <> 0 /\ (b = 0) /\ (SQRT n = 0)` by metis_tac[MULT_EQ_0, SQRT_0, DECIDE``~(0 < 0)``] >> 82 decide_tac, 83 `a <> 0 /\ b <> 0` by metis_tac[MULT_EQ_0] >> 84 spose_not_then strip_assume_tac >> 85 `b < SQRT n` by decide_tac >> 86 `0 < a /\ 0 < b /\ 0 < SQRT n` by decide_tac >> 87 `n < a * SQRT n` by rw[] >> 88 `a * SQRT n < SQRT n * SQRT n` by rw[] >> 89 `n < (SQRT n) ** 2` by metis_tac[LESS_TRANS, EXP_2] >> 90 `(SQRT n) ** 2 <= n` by rw[SQRT_PROPERTY] >> 91 decide_tac 92 ]); 93 94(* Theorem: (n = a * b) /\ SQRT n < a ==> b <= SQRT n *) 95(* Proof: 96 If n = 0, then a = 0 or b = 0 by MULT_EQ_0 97 But SQRT 0 = 0 by SQRT_0 98 so a <> 0, making b = 0, and b <= SQRT n true. 99 If n <> 0, a <> 0 and b <> 0 by MULT_EQ_0 100 By contradiction, suppose SQRT n < b. 101 Then SUC (SQRT n) ** 2 102 = SUC (SQRT n) * SUC (SQRT n) by EXP_2 103 <= a * SUC (SQRT n) by LT_MULT_RCANCEL, 0 < SUC (SQRT n). 104 <= a * b = n by LT_MULT_LCANCEL, 0 < a. 105 Giving (SUC (SQRT n)) ** 2 <= n by LESS_EQ_TRANS 106 or SQRT ((SUC (SQRT n)) ** 2) <= SQRT n by SQRT_LE 107 or SUC (SQRT n) <= SQRT n by SQRT_OF_SQ 108 which is a contradiction to !m. SUC m > m by LESS_SUC_REFL 109 *) 110val two_factors_property_2 = store_thm( 111 "two_factors_property_2", 112 ``!n a b. (n = a * b) /\ SQRT n < a ==> b <= SQRT n``, 113 rpt strip_tac >> 114 Cases_on `n = 0` >| [ 115 `a <> 0 /\ (b = 0) /\ (SQRT 0 = 0)` by metis_tac[MULT_EQ_0, SQRT_0, DECIDE``~(0 < 0)``] >> 116 decide_tac, 117 `a <> 0 /\ b <> 0` by metis_tac[MULT_EQ_0] >> 118 spose_not_then strip_assume_tac >> 119 `SQRT n < b` by decide_tac >> 120 `SUC (SQRT n) <= a /\ SUC (SQRT n) <= b` by decide_tac >> 121 `SUC (SQRT n) * SUC (SQRT n) <= a * SUC (SQRT n)` by rw[] >> 122 `a * SUC (SQRT n) <= n` by rw[] >> 123 `SUC (SQRT n) ** 2 <= n` by metis_tac[LESS_EQ_TRANS, EXP_2] >> 124 `SUC (SQRT n) <= SQRT n` by metis_tac[SQRT_LE, SQRT_OF_SQ] >> 125 decide_tac 126 ]); 127 128(* Theorem: (n = a * b) ==> a <= SQRT n \/ b <= SQRT n *) 129(* Proof: 130 By contradiction, suppose SQRT n < a /\ SQRT n < b. 131 Then (n = a * b) /\ SQRT n < a ==> b <= SQRT n by two_factors_property_2 132 which contradicts SQRT n < b. 133 *) 134val two_factors_property = store_thm( 135 "two_factors_property", 136 ``!n a b. (n = a * b) ==> a <= SQRT n \/ b <= SQRT n``, 137 rpt strip_tac >> 138 spose_not_then strip_assume_tac >> 139 `SQRT n < a` by decide_tac >> 140 metis_tac[two_factors_property_2]); 141 142(* ------------------------------------------------------------------------- *) 143(* Primality or Compositeness based on SQRT *) 144(* ------------------------------------------------------------------------- *) 145 146(* Theorem: prime p <=> 1 < p /\ !q. 1 < q /\ q <= SQRT p ==> ~(q divides p) *) 147(* Proof: 148 If part: prime p ==> 1 < p /\ !q. 1 < q /\ q <= SQRT p ==> ~(q divides p) 149 First one: prime p ==> 1 < p is true by ONE_LT_PRIME 150 Second one: by contradiction, suppose q divides p. 151 But prime p /\ q divides p ==> (q = p) or (q = 1) by prime_def 152 Given 1 < q, q <> 1, hence q = p. 153 This means p <= SQRT p, giving p = 0 or p = 1 by SQRT_GE_SELF 154 which contradicts p <> 0 and p <> 1 by PRIME_POS, prime_def 155 Only-if part: 1 < p /\ !q. 1 < q /\ q <= SQRT p ==> ~(q divides p) ==> prime p 156 By prime_def, this is to show: 157 (1) p <> 1, true since 1 < p. 158 (2) b divides p ==> (b = p) \/ (b = 1) 159 By contradiction, suppose b <> p /\ b <> 1. 160 b divides p ==> ?a. p = a * b by divides_def 161 which means a <= SQRT p \/ b <= SQRT p by two_factors_property 162 If a <= SQRT p, 163 1 < p ==> p <> 0, so a <> 0 by MULT 164 also b <> p ==> a <> 1 by MULT_LEFT_1 165 so 1 < a, and a divides p by prime_def, MULT_COMM 166 The implication gives ~(a divides p), a contradiction. 167 If b <= SQRT p, 168 1 < p ==> p <> 0, so b <> 0 by MULT_0 169 also b <> 1, so 1 < b, and b divides p. 170 The implication gives ~(b divides p), a contradiction. 171 *) 172val prime_by_sqrt_factors = store_thm( 173 "prime_by_sqrt_factors", 174 ``!p. prime p <=> 1 < p /\ !q. 1 < q /\ q <= SQRT p ==> ~(q divides p)``, 175 rw[EQ_IMP_THM] >- 176 rw[ONE_LT_PRIME] >- 177 (spose_not_then strip_assume_tac >> 178 `0 < p` by rw[PRIME_POS] >> 179 `p <> 0 /\ q <> 1` by decide_tac >> 180 `(q = p) /\ p <> 1` by metis_tac[prime_def] >> 181 metis_tac[SQRT_GE_SELF]) >> 182 rw[prime_def] >> 183 spose_not_then strip_assume_tac >> 184 `?a. p = a * b` by rw[GSYM divides_def] >> 185 `a <= SQRT p \/ b <= SQRT p` by rw[two_factors_property] >| [ 186 `a <> 1` by metis_tac[MULT_LEFT_1] >> 187 `p <> 0` by decide_tac >> 188 `a <> 0` by metis_tac[MULT] >> 189 `1 < a` by decide_tac >> 190 metis_tac[divides_def, MULT_COMM], 191 `p <> 0` by decide_tac >> 192 `b <> 0` by metis_tac[MULT_0] >> 193 `1 < b` by decide_tac >> 194 metis_tac[] 195 ]); 196 197(* Theorem: 1 < n ==> (~prime n <=> ?p. prime p /\ p divides n /\ p <= SQRT n) *) 198(* Proof: 199 If part ~prime n ==> ?p. prime p /\ p divides n /\ p <= SQRT n 200 Given n <> 1, ?p. prime p /\ p divides n by PRIME_FACTOR 201 If p <= SQRT n, take this p. 202 If ~(p <= SQRT n), SQRT n < p. 203 Since p divides n, ?q. n = p * q by divides_def, MULT_COMM 204 Hence q <= SQRT n by two_factors_property_2 205 Since prime p but ~prime n, q <> 1 by MULT_RIGHT_1 206 so ?t. prime t /\ t divides q by PRIME_FACTOR 207 Since 1 < n, n <> 0, so q <> 0 by MULT_0 208 so t divides q ==> t <= q by DIVIDES_LE, 0 < q. 209 Take t, then t divides n by DIVIDES_TRANS 210 and t <= SQRT n by LESS_EQ_TRANS 211 Only-if part: ?p. prime p /\ p divides n /\ p <= SQRT n ==> ~prime n 212 By contradiction, assume prime n. 213 Then p divides n ==> p = 1 or p = n by prime_def 214 But prime p ==> p <> 1, so p = n by ONE_LT_PRIME 215 Giving p <= SQRT p, 216 thus forcing p = 0 or p = 1 by SQRT_GE_SELF 217 Both are impossible for prime p. 218*) 219val prime_factor_estimate = store_thm( 220 "prime_factor_estimate", 221 ``!n. 1 < n ==> (~prime n <=> ?p. prime p /\ p divides n /\ p <= SQRT n)``, 222 rpt strip_tac >> 223 `n <> 1` by decide_tac >> 224 rw[EQ_IMP_THM] >| [ 225 `?p. prime p /\ p divides n` by rw[PRIME_FACTOR] >> 226 Cases_on `p <= SQRT n` >- 227 metis_tac[] >> 228 `SQRT n < p` by decide_tac >> 229 `?q. n = q * p` by rw[GSYM divides_def] >> 230 `_ = p * q` by rw[MULT_COMM] >> 231 `q <= SQRT n` by metis_tac[two_factors_property_2] >> 232 `q <> 1` by metis_tac[MULT_RIGHT_1] >> 233 `?t. prime t /\ t divides q` by rw[PRIME_FACTOR] >> 234 `n <> 0` by decide_tac >> 235 `q <> 0` by metis_tac[MULT_0] >> 236 `0 < q ` by decide_tac >> 237 `t <= q` by rw[DIVIDES_LE] >> 238 `q divides n` by metis_tac[divides_def] >> 239 metis_tac[DIVIDES_TRANS, LESS_EQ_TRANS], 240 spose_not_then strip_assume_tac >> 241 `1 < p` by rw[ONE_LT_PRIME] >> 242 `p <> 1 /\ p <> 0` by decide_tac >> 243 `p = n` by metis_tac[prime_def] >> 244 metis_tac[SQRT_GE_SELF] 245 ]); 246 247(* ------------------------------------------------------------------------- *) 248(* Primality Testing Algorithm *) 249(* ------------------------------------------------------------------------- *) 250 251(* Seek the prime factor of number n, starting with q, up to cutoff c. *) 252val factor_seek_def = tDefine "factor_seek" ` 253 factor_seek n c q = 254 if c <= q then n 255 else if 1 < q /\ (n MOD q = 0) then q 256 else factor_seek n c (q + 1) 257`(WF_REL_TAC `measure (\(n,c,q). c - q)` >> simp[]); 258(* Use 1 < q so that, for prime n, it gives a result n for any initial q, including q = 1. *) 259 260(* Primality test by seeking a factor exceeding (SQRT n). *) 261val prime_test_def = Define` 262 prime_test n = 263 if n <= 1 then F 264 else factor_seek n (1 + SQRT n) 2 = n 265`; 266 267(* 268EVAL ``MAP prime_test [1 .. 15]``; = [F; T; T; F; T; F; T; F; F; F; T; F; T; F; F]: thm 269*) 270 271(* Theorem: 0 < n ==> factor_seek n c q <= n *) 272(* Proof: 273 By induction from factor_seek_def. 274 If c <= q, 275 Then factor_seek n c q = n, hence true by factor_seek_def 276 If q = 0, 277 Then factor_seek n c 0 = 0, hence true by factor_seek_def 278 If n MOD q = 0, 279 Then factor_seek n c q = q by factor_seek_def 280 Thus q divides n by DIVIDES_MOD_0, q <> 0 281 hence q <= n by DIVIDES_LE, 0 < n 282 Otherwise, 283 factor_seek n c q 284 = factor_seek n c (q + 1) by factor_seek_def 285 <= n by induction hypothesis 286*) 287val factor_seek_bound = store_thm( 288 "factor_seek_bound", 289 ``!n c q. 0 < n ==> factor_seek n c q <= n``, 290 ho_match_mp_tac (theorem "factor_seek_ind") >> 291 rw[] >> 292 rw[Once factor_seek_def] >> 293 `q divides n` by rw[DIVIDES_MOD_0] >> 294 rw[DIVIDES_LE]); 295 296(* Theorem: 1 < q /\ q <= c /\ c <= n ==> 297 ((factor_seek n c q = n) <=> (!p. q <= p /\ p < c ==> ~(p divides n))) *) 298(* Proof: 299 By induction from factor_seek_def, this is to show: 300 (1) n MOD q = 0 ==> ?p. (q <= p /\ p < c) /\ p divides n 301 Take p = q, then n MOD q = 0 ==> q divides n by DIVIDES_MOD_0, 0 < q 302 (2) n MOD q <> 0 ==> factor_seek n c (q + 1) = n <=> 303 !p. q <= p /\ p < c ==> ~(p divides n) 304 factor_seek n c (q + 1) = n 305 <=> !p. q + 1 <= p /\ p < c ==> ~(p divides n)) by induction hypothesis 306 or !p. q < p /\ p < c ==> ~(p divides n)) 307 But n MOD q <> 0 gives ~(q divides n) by DIVIDES_MOD_0, 0 < q 308 Thus !p. q <= p /\ p < c ==> ~(p divides n)) 309*) 310val factor_seek_thm = store_thm( 311 "factor_seek_thm", 312 ``!n c q. 1 < q /\ q <= c /\ c <= n ==> 313 ((factor_seek n c q = n) <=> (!p. q <= p /\ p < c ==> ~(p divides n)))``, 314 ho_match_mp_tac (theorem "factor_seek_ind") >> 315 rw[] >> 316 rw[Once factor_seek_def] >| [ 317 qexists_tac `q` >> 318 rw[DIVIDES_MOD_0], 319 rw[EQ_IMP_THM] >> 320 spose_not_then strip_assume_tac >> 321 `0 < q` by decide_tac >> 322 `p <> q` by metis_tac[DIVIDES_MOD_0] >> 323 `q + 1 <= p` by decide_tac >> 324 metis_tac[] 325 ]); 326 327(* Theorem: prime n = prime_test n *) 328(* Proof: 329 prime n 330 <=> 1 < n /\ !q. 1 < q /\ n <= SQRT n ==> ~(n divides p) by prime_by_sqrt_factors 331 <=> 1 < n /\ !q. 2 <= q /\ n < c ==> ~(n divides p) where c = 1 + SQRT n 332 Under 1 < n, 333 Note SQRT n <> 0 by SQRT_EQ_0, n <> 0 334 so 1 < 1 + SQRT n = c, or 2 <= c. 335 Also SQRT n <= n by SQRT_LE_SELF 336 but SQRT n <> n by SQRT_EQ_SELF, 1 < n 337 so SQRT n < n, or c <= n. 338 Thus 1 < n /\ !q. 2 <= q /\ n < c ==> ~(n divides p) 339 <=> factor_seek n c q = n by factor_seek_thm 340 <=> prime_test n by prime_test_def 341*) 342val prime_test_thm = store_thm( 343 "prime_test_thm", 344 ``!n. prime n = prime_test n``, 345 rw[prime_test_def, prime_by_sqrt_factors] >> 346 rw[EQ_IMP_THM] >| [ 347 qabbrev_tac `c = SQRT n + 1` >> 348 `0 < 2` by decide_tac >> 349 `SQRT n <> 0` by rw[SQRT_EQ_0] >> 350 `2 <= c` by rw[Abbr`c`] >> 351 `SQRT n <= n` by rw[SQRT_LE_SELF] >> 352 `SQRT n <> n` by rw[SQRT_EQ_SELF] >> 353 `c <= n` by rw[Abbr`c`] >> 354 `!q. 2 <= q /\ q < c ==> ~(q divides n)` by fs[Abbr`c`] >> 355 rw[factor_seek_thm], 356 qabbrev_tac `c = SQRT n + 1` >> 357 `0 < 2` by decide_tac >> 358 `SQRT n <> 0` by rw[SQRT_EQ_0] >> 359 `2 <= c` by rw[Abbr`c`] >> 360 `SQRT n <= n` by rw[SQRT_LE_SELF] >> 361 `SQRT n <> n` by rw[SQRT_EQ_SELF] >> 362 `c <= n` by rw[Abbr`c`] >> 363 fs[factor_seek_thm] >> 364 `!p. 1 < p /\ p <= SQRT n ==> ~(p divides n)` by fs[Abbr`c`] >> 365 rw[] 366 ]); 367 368 369(* ------------------------------------------------------------------------- *) 370 371(* export theory at end *) 372val _ = export_theory(); 373 374(*===========================================================================*) 375