1(* ------------------------------------------------------------------------- *) 2(* Windmills of the minds: Fermat's Two Squares Theorem *) 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 "twoSquares"; 12 13(* ------------------------------------------------------------------------- *) 14 15 16(* open dependent theories *) 17(* arithmeticTheory -- load by default *) 18 19(* val _ = load "windmillTheory"; *) 20open helperTheory; 21open helperNumTheory; 22open helperSetTheory; 23open helperFunctionTheory; 24open arithmeticTheory pred_setTheory; 25open dividesTheory; (* for divides_def, prime_def *) 26 27open windmillTheory; 28 29(* val _ = load "involuteFixTheory"; *) 30open involuteTheory; (* for involute_bij *) 31open involuteFixTheory; 32 33(* val _ = load "iterateComposeTheory"; *) 34open iterateTheory; 35open iterateComposeTheory; 36 37(* val _ = load "iterateComputeTheory"; *) 38open iterateComputeTheory; 39 40(* for later computation *) 41open listTheory; 42open rich_listTheory; (* for MAP_REVERSE *) 43open helperListTheory; (* for listRangeINC_LEN *) 44open listRangeTheory; (* for listRangeINC_CONS *) 45 46(* for group action *) 47(* val _ = load "involuteActionTheory"; *) 48open involuteActionTheory; 49open groupActionTheory; 50open groupInstancesTheory; 51 52 53(* ------------------------------------------------------------------------- *) 54(* Windmills of the minds Documentation *) 55(* ------------------------------------------------------------------------- *) 56(* Overloading: 57*) 58(* 59 60 Helper Theorems: 61 62 Two Squares Fixes: 63 zagier_fixes_prime 64 |- !p. prime p /\ (p MOD 4 = 1) ==> 65 (fixes zagier (mills p) = {(1,1,p DIV 4)}) 66 67 Fermat Two-Squares Uniqueness: 68 fermat_two_squares_unique_thm 69 |- !p a b c d. prime p /\ (p = a ** 2 + b ** 2) /\ 70 (p = c ** 2 + d ** 2) ==> ({a; b} = {c; d}) 71 fermat_two_squares_unique_odd_even 72 |- !p a b c d. prime p /\ 73 ODD a /\ EVEN b /\ (p = a ** 2 + b ** 2) /\ 74 ODD c /\ EVEN d /\ (p = c ** 2 + d ** 2) ==> 75 (a = c) /\ (b = d) 76 flip_fixes_prime_card_upper 77 |- !p. prime p /\ (p MOD 4 = 1) ==> CARD (fixes flip (mills p)) <= 1 78 79 80 Fermat Two-Squares Existence: 81 fermat_two_squares_exists_windmill 82 |- !p. prime p /\ (p MOD 4 = 1) ==> ?x y. p = windmill x y y 83 fermat_two_squares_exists_odd_even 84 |- !p. prime p /\ (p MOD 4 = 1) ==> 85 ?u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) 86 87 88 Fermat Two-Squares Theorem: 89 flip_fixes_prime_sing 90 |- !p. prime p /\ (p MOD 4 = 1) ==> SING (fixes flip (mills p)) 91 flip_fixes_prime|- !p. prime p /\ (p MOD 4 = 1) ==> 92 (let u = (1,1,p DIV 4) ; 93 n = iterate_period (zagier o flip) u 94 in fixes flip (mills p) = {FUNPOW (zagier o flip) (HALF n) u}) 95 flip_fixes_prime_alt 96 |- !p u n. prime p /\ (p MOD 4 = 1) /\ (u = (1,1,p DIV 4)) /\ 97 (n = iterate_period (zagier o flip) u) ==> 98 (fixes flip (mills p) = {FUNPOW (zagier o flip) (HALF n) u} 99 fermat_two_squares_thm 100 |- !p. prime p /\ (p MOD 4 = 1) ==> 101 ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) 102 fermat_two_squares_iff 103 |- !p. prime p ==> ((p MOD 4 = 1) <=> 104 ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2)) 105 106 Fermat Two-Squares Algorithm: 107 zagier_flip_1_1_z_period 108 |- !z. (iterate_period (zagier o flip) (1,1,z) = 1) <=> (z = 1) 109 flip_fixes_iterates_prime 110 |- !p u n g. prime p /\ (p MOD 4 = 1) /\ (u = (1,1,p DIV 4)) /\ 111 (n = iterate_period (zagier o flip) u) /\ 112 (g = (\(x,y,z). y <> z)) ==> 113 ~g (FUNPOW (zagier o flip) (HALF n) u) /\ 114 !j. j < HALF n ==> g (FUNPOW (zagier o flip) j u) 115 Computation by WHILE loop: 116 found_def |- !x y z. found (x,y,z) <=> (y = z) 117 found_not |- $~ o found = (\(x,y,z). y <> z) 118 two_sq_def |- !n. two_sq n = WHILE ($~ o found) (zagier o flip) (1,1,n DIV 4) 119 two_sq_alt |- !n. two_sq n = WHILE (\(x,y,z). y <> z) (zagier o flip) (1,1,n DIV 4) 120 two_sq_thm |- !p. prime p /\ (p MOD 4 = 1) ==> two_sq p IN fixes flip (mills p) 121 two_sq_while_hoare 122 |- !p. prime p /\ (p MOD 4 = 1) ==> 123 HOARE_SPEC (fixes zagier (mills p)) 124 (WHILE ($~ o found) (zagier o flip)) 125 (fixes flip (mills p)) 126 two_sq_while_thm|- !p. prime p /\ (p MOD 4 = 1) ==> 127 (let (x,y,z) = WHILE ($~ o found) (zagier o flip) (1,1,p DIV 4) 128 in p = x ** 2 + (y + z) ** 2) 129 two_squares_def |- !n. two_squares n = (let (x,y,z) = two_sq n in (x,y + z)) 130 two_squares_thm |- !p. prime p /\ (p MOD 4 = 1) ==> 131 (let (u,v) = two_squares p in p = u ** 2 + v ** 2) 132 133 Fermat's Two Squares Theorem by Group Action: 134 135 Relation to Fixes and Pairs: 136 involute_fixed_points_eq_fixes 137 |- !f X. f involute X ==> 138 fixed_points (FUNPOW f) Z2 X = fixes f X 139 involute_multi_orbits_union_eq_pairs 140 |- !f X. f involute X ==> 141 BIGUNION (multi_orbits (FUNPOW f) Z2 X) = pairs f X 142 143 Fermat's Two-Squares Theorem (Existence): 144 zagier_fixed_points_sing 145 |- !p. prime p /\ (p MOD 4 = 1) ==> 146 (fixed_points (FUNPOW zagier) Z2 (mills p) = {(1,1,p DIV 4)}) 147 fermat_two_squares_exists_alt 148 |- !p. prime p /\ (p MOD 4 = 1) ==> ?u v. p = u ** 2 + v ** 2 149 150*) 151 152(* ------------------------------------------------------------------------- *) 153(* Helper Theorems *) 154(* ------------------------------------------------------------------------- *) 155 156(* ------------------------------------------------------------------------- *) 157(* Two Squares Fixes. *) 158(* ------------------------------------------------------------------------- *) 159 160 161(* Theorem: prime p /\ (p MOD 4 = 1) ==> (fixes zagier (mills p) = {(1,1,p DIV 4)}) *) 162(* Proof: 163 Let s = mills p, 164 k = p DIV 4, 165 a = fixes zagier s. 166 The goal becomes: a = {(1,1,k)} 167 168 Use EXTENSION to show two sets are equal: 169 (1) (x,y,z) IN a ==> x = 1, y = 1, z = k. 170 Note (x,y,z) IN s /\ 171 zagier (x,y,z) = (x,y,z) by fixes_element 172 but x <> 0 by mills_prime_triple_nonzero 173 so x = y by zagier_fix 174 Now p = windmill x x z by mills_element 175 ==> x = 1, y = 1, z = k by windmill_trivial_prime 176 (2) (1,1,k) IN a 177 Note p = k * 4 + p MOD 4 by DIVISION 178 = 4 * k + 1 by p MOD 4 = 1 179 = windmill 1 1 k by windmill_trivial 180 so (1,1,k) IN s by mills_element 181 and zagier (1,1,k) = (1,1,k) by zagier_1_1_z 182 ==> (1,1,k) IN a by fixes_element 183*) 184Theorem zagier_fixes_prime: 185 !p. prime p /\ (p MOD 4 = 1) ==> (fixes zagier (mills p) = {(1,1,p DIV 4)}) 186Proof 187 rpt strip_tac >> 188 qabbrev_tac `s = mills p` >> 189 qabbrev_tac `k = p DIV 4` >> 190 qabbrev_tac `a = fixes zagier s` >> 191 simp[EXTENSION, pairTheory.FORALL_PROD, EQ_IMP_THM] >> 192 ntac 5 strip_tac >| [ 193 `(p_1,p_1',p_2) IN s /\ (zagier (p_1,p_1',p_2) = (p_1,p_1',p_2))` by metis_tac[fixes_element] >> 194 `p_1 = p_1'` by metis_tac[zagier_fix, mills_prime_triple_nonzero] >> 195 metis_tac[windmill_trivial_prime, mills_element], 196 `p = k * 4 + p MOD 4` by rw[DIVISION, Abbr`k`] >> 197 `_ = windmill 1 1 k` by rw[windmill_trivial] >> 198 `(1,1,k) IN s` by rw[mills_element, Abbr`s`] >> 199 fs[fixes_element, zagier_1_1_z, Abbr`a`] 200 ] 201QED 202 203(* ------------------------------------------------------------------------- *) 204(* Fermat Two-Squares Uniqueness. *) 205(* ------------------------------------------------------------------------- *) 206 207 208(* Theorem: prime p /\ (p = a ** 2 + b ** 2) /\ (p = c ** 2 + d ** 2) ==> 209 ({a; b} = {c; d}) *) 210(* Proof: 211 This is to show: 212 (a = c) /\ (b = d) 213 \/ (a = d) /\ (b = c) by doublet_eq 214 215 Step 1: basic properties 216 Note 0 < p by NOT_PRIME_0 217 and 0 < a /\ 0 < b /\ 0 < c /\ 0 < d by SQ_0, prime_non_square, square_def 218 and a * d < p by squares_sum_inequality, ADD_COMM 219 and b * c < p by squares_sum_inequality, ADD_COMM 220 221 Step 2: use identities and prime divisibility 222 Note (a * d - b * c) * (a * d + b * c) 223 = (d ** 2 - b ** 2) * p by squares_sum_identity_1 224 and (b * c - a * d) * (a * d + b * c) 225 = (b ** 2 - d ** 2) * p by squares_sum_identity_2 226 Thus p divides (a * d - b * c) * (a * d + b * c) by divides_def 227 and p divides (b * c - a * d) * (a * d + b * c) by divides_def 228 ==> (p divides (a * d - b * c) /\ p divides (b * c - a * d)) 229 \/ p divides (a * d + b * c) by euclid_prime 230 231 Case: p divides a * d - b * c /\ p divides b * c - a * d 232 Note a * d - b * c < p by a * d < p 233 and b * c - a * d < p by b * c < p 234 Thus (a * d - b * c = 0) by DIVIDES_LEQ_OR_ZERO 235 and (b * c - a * d = 0) by DIVIDES_LEQ_OR_ZERO 236 ==> a * d = b * c by arithmetic 237 Hence (a = c) and (b = d) by squares_sum_prime_thm 238 239 Case: p divides a * d + b * c 240 Note 0 < (a * d + b * c) < 2 * p by a * d < p, b * c < p, and a,b,c,d <> 0 241 Thus a * d + b * c = p by divides_eq_thm 242 For four_squares_identity, break into cases: 243 If b * d <= a * c, 244 Note p ** 2 245 = p * p by EXP_2 246 = (a ** 2 + b ** 2) * (c ** 2 + d ** 2) 247 = p ** 2 + (a * c - b * d) ** 2 248 by four_squares_identity 249 Thus (a * c - b * d) ** 2 = 0 by EQ_ADD_LCANCEL 250 or a * c - b * d = 0 by EXP_2_EQ_0 251 ==> a * c = b * d by b * d <= a * c 252 ==> (a = c) /\ (b = d) by squares_sum_prime_thm, ADD_COMM 253 If a * c <= b * d, 254 Note p ** 2 255 = p * p by EXP_2 256 = (a ** 2 + b ** 2) * (c ** 2 + d ** 2) 257 = p ** 2 + (b * d - a * c) ** 2 258 by four_squares_identity_alt 259 Thus (b * d - a * c) ** 2 = 0 by EQ_ADD_LCANCEL 260 or b * d - a * c = 0 by EXP_2_EQ_0 261 ==> a * c = b * d by a * c <= b * d 262 ==> (a = c) /\ (b = d) by squares_sum_prime_thm, ADD_COMM 263*) 264Theorem fermat_two_squares_unique_thm: 265 !p a b c d. prime p /\ (p = a ** 2 + b ** 2) /\ (p = c ** 2 + d ** 2) ==> 266 ({a; b} = {c; d}) 267Proof 268 rpt strip_tac >> 269 simp[doublet_eq] >> 270 spose_not_then strip_assume_tac >> 271 `0 < p` by metis_tac[NOT_PRIME_0, NOT_ZERO] >> 272 `0 < a /\ 0 < b /\ 0 < c /\ 0 < d` by metis_tac[SQ_0, prime_non_square, square_def, EXP_2, ADD, ADD_0, NOT_ZERO] >> 273 `a * d < p` by metis_tac[squares_sum_inequality, ADD_COMM] >> 274 `b * c < p` by metis_tac[squares_sum_inequality, ADD_COMM] >> 275 `(a * d - b * c) * (a * d + b * c) = (d ** 2 - b ** 2) * p` by rw[squares_sum_identity_1] >> 276 `(b * c - a * d) * (a * d + b * c) = (b ** 2 - d ** 2) * p` by rw[squares_sum_identity_2] >> 277 `p divides (a * d - b * c) * (a * d + b * c)` by metis_tac[divides_def] >> 278 `p divides (b * c - a * d) * (a * d + b * c)` by metis_tac[divides_def] >> 279 `(p divides (a * d - b * c) /\ p divides (b * c - a * d)) \/ p divides (a * d + b * c)` by metis_tac[euclid_prime] >| [ 280 `a * d - b * c < p /\ b * c - a * d < p` by fs[] >> 281 `(a * d - b * c = 0) /\ (b * c - a * d = 0)` by metis_tac[DIVIDES_LEQ_OR_ZERO, NOT_LESS] >> 282 `a * d = b * c` by fs[] >> 283 metis_tac[squares_sum_prime_thm], 284 `a * d + b * c = p` by 285 (`a * d + b * c < 2 * p` by fs[] >> 286 `0 < a * d + b * c` by metis_tac[MULT, ADD_EQ_0, MULT_EQ_0, NOT_ZERO] >> 287 rw[divides_eq_thm]) >> 288 Cases_on `b * d <= a * c` >| [ 289 `p ** 2 = (a ** 2 + b ** 2) * (c ** 2 + d ** 2)` by fs[] >> 290 `_ = p ** 2 + (a * c - b * d) ** 2` by rfs[four_squares_identity] >> 291 `(a * c - b * d) ** 2 = 0` by metis_tac[EQ_ADD_LCANCEL, ADD_0] >> 292 `a * c - b * d = 0` by metis_tac[EXP_2_EQ_0] >> 293 `a * c = b * d` by fs[] >> 294 metis_tac[squares_sum_prime_thm, ADD_COMM], 295 `p ** 2 = (a ** 2 + b ** 2) * (c ** 2 + d ** 2)` by fs[] >> 296 `_ = p ** 2 + (b * d - a * c) ** 2` by rfs[four_squares_identity_alt] >> 297 `(b * d - a * c) ** 2 = 0` by metis_tac[EQ_ADD_LCANCEL, ADD_0] >> 298 `b * d - a * c = 0` by metis_tac[EXP_2_EQ_0] >> 299 `a * c = b * d` by fs[] >> 300 metis_tac[squares_sum_prime_thm, ADD_COMM] 301 ] 302 ] 303QED 304 305 306(* Theorem: prime p /\ ODD a /\ EVEN b /\ (p = a ** 2 + b ** 2) /\ 307 ODD c /\ EVEN d /\ (p = c ** 2 + d ** 2) ==> ((a = c) /\ (b = d)) *) 308(* Proof: 309 Note {a; b} = {c; d} by fermat_two_squares_unique_thm 310 Thus (a = c) /\ (b = d) \/ (a = d) /\ (b = c) by EXTENSION 311 But (a <> d) \/ (b <> c) by EVEN_ODD 312 so (a = c) /\ (b = d). 313*) 314Theorem fermat_two_squares_unique_odd_even: 315 !p a b c d. prime p /\ ODD a /\ EVEN b /\ (p = a ** 2 + b ** 2) /\ 316 ODD c /\ EVEN d /\ (p = c ** 2 + d ** 2) ==> ((a = c) /\ (b = d)) 317Proof 318 ntac 7 strip_tac >> 319 `{a; b} = {c; d}` by fs[fermat_two_squares_unique_thm] >> 320 fs[EXTENSION] >> 321 metis_tac[EVEN_ODD] 322QED 323 324 325(* Theorem: prime p /\ (p MOD 4 = 1) ==> CARD (fixes flip (mills p)) <= 1 *) 326(* Proof: 327 Let s = mills p, 328 b = fixes flip s. 329 The goal is: CARD b <= 1. 330 By contradiction, suppose CARD b > 1. 331 Then CARD b <> 0, 332 so ?u. u IN b by CARD_EMPTY, MEMBER_NOT_EMPTY 333 Note u IN s by fixes_element 334 and ?x y. u = (x,y,y) by triple_parts, flip_fix 335 Thus p = windmill x y y by mills_element 336 337 Also CARD b <> 1, 338 Then ~SING b by SING_CARD_1 339 so b <> EMPTY by MEMBER_NOT_EMPTY 340 and ?v. v IN b /\ v <> u by SING_ONE_ELEMENT, b <> EMPTY 341 Note v IN s by fixes_element 342 and ?h k. v = (h,k,k) by triple_parts, flip_fix 343 Thus p = windmill h k k by mills_element 344 345 Let c = 2 * y, 346 d = 2 * k. 347 Now ODD p by odd_by_mod_4, p MOD 4 = 1 348 and (p = x ** 2 + c ** 2) /\ ODD x by windmill_x_y_y 349 and (p = h ** 2 + d ** 2) /\ ODD h by windmill_x_y_y 350 but EVEN c /\ EVEN d by EVEN_DOUBLE 351 ==> (x = h) /\ (c = d) by fermat_two_squares_unique_odd_even 352 so (x = h) /\ (y = k) by EQ_MULT_LCANCEL 353 or v = u by PAIR_EQ 354 This contradicts v <> u. 355*) 356Theorem flip_fixes_prime_card_upper: 357 !p. prime p /\ (p MOD 4 = 1) ==> CARD (fixes flip (mills p)) <= 1 358Proof 359 rpt strip_tac >> 360 qabbrev_tac `s = mills p` >> 361 qabbrev_tac `b = fixes flip s` >> 362 spose_not_then strip_assume_tac >> 363 `CARD b <> 0 /\ CARD b <> 1` by decide_tac >> 364 `?u. u IN b` by metis_tac[CARD_EMPTY, MEMBER_NOT_EMPTY] >> 365 `u IN s /\ ?x y. u = (x,y,y)` by metis_tac[fixes_element, triple_parts, flip_fix] >> 366 `p = windmill x y y` by fs[mills_element, Abbr`s`] >> 367 `~SING b` by metis_tac[SING_CARD_1] >> 368 `b <> EMPTY` by metis_tac[MEMBER_NOT_EMPTY] >> 369 `?v. v IN b /\ v <> u` by metis_tac[SING_ONE_ELEMENT] >> 370 `v IN s /\ ?h k. v = (h,k,k)` by metis_tac[fixes_element, triple_parts, flip_fix] >> 371 `p = windmill h k k` by fs[mills_element, Abbr`s`] >> 372 qabbrev_tac `c = 2 * y` >> 373 qabbrev_tac `d = 2 * k` >> 374 `ODD p` by rfs[odd_by_mod_4] >> 375 `(p = x ** 2 + c ** 2) /\ ODD x` by metis_tac[windmill_x_y_y] >> 376 `(p = h ** 2 + d ** 2) /\ ODD h` by metis_tac[windmill_x_y_y] >> 377 `EVEN c /\ EVEN d` by rw[EVEN_DOUBLE, Abbr`c`, Abbr`d`] >> 378 `(x = h) /\ (c = d)` by metis_tac[fermat_two_squares_unique_odd_even] >> 379 `y = k` by fs[Abbr`c`, Abbr`d`] >> 380 metis_tac[pairTheory.PAIR_EQ] 381QED 382 383(* ------------------------------------------------------------------------- *) 384(* Fermat Two-Squares Existence. *) 385(* ------------------------------------------------------------------------- *) 386 387 388(* Theorem: prime p /\ (p MOD 4 = 1) ==> ?x y. p = windmill x y y *) 389(* Proof: 390 Let m = mills p, the solutions (x,y,z) of p = x ** 2 + 4 * y * z. 391 Note ~square p by prime_non_square 392 so FINITE m by mills_non_square_finite 393 and zagier involute m by zagier_closure, zagier_involute, triple_parts, 394 mills_prime_triple_nonzero 395 and flip involute m by flip_closure, flip_involute, triple_parts 396 397 Now work out fixed points: 398 Let a = fixes zagier m, b = fixes flip m. 399 Let k = p DIV 4. 400 Then a = {(1,1,k)} by zagier_fixes_prime 401 402 The punch line: 403 Note ODD (CARD a) <=> ODD (CARD b) by involute_two_fixes_both_odd 404 now CARD a = 1 by CARD_SING 405 so ODD (CARD b) by ODD_1 406 ==> CARD b <> 0 by ODD 407 or b <> EMPTY by CARD_EMPTY 408 Thus ? (x, y, z) IN b by MEMBER_NOT_EMPTY 409 and (x, y, z) IN m /\ 410 (flip (x, y, z) = (x, y, z)) by fixes_element 411 so y = z by flip_fix 412 or (x,y,y) in m by above 413 Thus p = windmill x y y by mills_element 414*) 415Theorem fermat_two_squares_exists_windmill: 416 !p. prime p /\ (p MOD 4 = 1) ==> ?x y. p = windmill x y y 417Proof 418 rpt strip_tac >> 419 qabbrev_tac `m = mills p` >> 420 `~square p` by metis_tac[prime_non_square] >> 421 `FINITE m` by fs[mills_non_square_finite, Abbr`m`] >> 422 `zagier involute m` by metis_tac[zagier_closure, zagier_involute, triple_parts, mills_prime_triple_nonzero] >> 423 `flip involute m` by metis_tac[flip_closure, flip_involute, triple_parts] >> 424 qabbrev_tac `a = fixes zagier m` >> 425 qabbrev_tac `b = fixes flip m` >> 426 qabbrev_tac `k = p DIV 4` >> 427 `a = {(1,1,k)}` by rw[zagier_fixes_prime, Abbr`a`, Abbr`m`] >> 428 `CARD a = 1` by rw[] >> 429 `ODD (CARD a) <=> ODD (CARD b)` by rw[involute_two_fixes_both_odd, Abbr`a`, Abbr`b`] >> 430 `ODD (CARD b)` by metis_tac[ODD_1] >> 431 `CARD b <> 0` by metis_tac[ODD] >> 432 `b <> EMPTY` by metis_tac[CARD_EMPTY] >> 433 `?t. t IN b` by rw[MEMBER_NOT_EMPTY] >> 434 `?x y z. t = (x, y, z)` by rw[triple_parts] >> 435 `t IN m /\ (flip (x, y, z) = (x, y, z))` by metis_tac[fixes_element] >> 436 `y = z` by fs[flip_fix] >> 437 metis_tac[mills_element] 438QED 439 440 441(* Theorem: prime p /\ (p MOD 4 = 1) ==> ?u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) *) 442(* Proof: 443 Note ?x y. p = windmill x y y by fermat_two_squares_exists_windmill 444 = x ** 2 + (2 * y) ** 2 by windmill_by_squares 445 Put u = x, v = 2 * y. 446 Then p = u ** 2 + v ** 2 by above 447 = u * u + v * v by EXP_2 448 or p - u * u = v * v by arithmetic 449 and EVEN v by EVEN_DOUBLE 450 and EVEN (v * v) by EVEN_MULT 451 Now ODD p by odd_by_mod_4, p MOD 4 = 1 452 so ODD (u * u) by EVEN_SUB, EVEN_ODD 453 or ODD u by ODD_MULT 454*) 455Theorem fermat_two_squares_exists_odd_even: 456 !p. prime p /\ (p MOD 4 = 1) ==> 457 ?u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) 458Proof 459 rpt strip_tac >> 460 `?x y. p = windmill x y y` by rw[fermat_two_squares_exists_windmill] >> 461 `_ = x ** 2 + (2 * y) ** 2` by rw[windmill_by_squares] >> 462 qabbrev_tac `u = x` >> 463 qabbrev_tac `v = 2 * y` >> 464 `p = u * u + v * v` by simp[] >> 465 `v * v = p - u * u` by decide_tac >> 466 `EVEN v` by rw[EVEN_DOUBLE, Abbr`v`] >> 467 `EVEN (v * v)` by rw[EVEN_MULT] >> 468 `ODD p` by rw[odd_by_mod_4] >> 469 `u * u <= p` by decide_tac >> 470 `ODD (u * u)` by metis_tac[EVEN_SUB, EVEN_ODD] >> 471 metis_tac[ODD_MULT] 472QED 473 474(* ------------------------------------------------------------------------- *) 475(* Fermat Two-Squares Theorem. *) 476(* ------------------------------------------------------------------------- *) 477 478 479(* Theorem: prime p /\ (p MOD 4 = 1) ==> SING (fixes flip (mills p)) *) 480(* Proof: 481 Let s = mills p, 482 b = fixes flip s. 483 The goal is: SING b. 484 Note ?x y. p = windmill x y y 485 by fermat_two_squares_exists_windmill 486 Let u = (x,y,y). 487 Then u IN s by mills_element 488 and u IN b by fixes_element, flip_fix 489 By SING_DEF, EXTENSION, this is to show: 490 !t. t IN b ==> t = u. 491 Note t IN s by fixes_element 492 and ?h k. t = (h,k,k) by triple_parts, flip_fix 493 Thus p = windmill h k k by mills_element 494 495 Let c = 2 * y, 496 d = 2 * k. 497 Now ODD p by odd_by_mod_4, p MOD 4 = 1 498 and (p = x ** 2 + c ** 2) /\ ODD x by windmill_x_y_y 499 and (p = h ** 2 + d ** 2) /\ ODD h by windmill_x_y_y 500 but EVEN c /\ EVEN d by EVEN_DOUBLE 501 ==> (x = h) /\ (c = d) by fermat_two_squares_unique_odd_even 502 so (x = h) /\ (y = k) by EQ_MULT_LCANCEL 503 or t = u. 504*) 505Theorem flip_fixes_prime_sing: 506 !p. prime p /\ (p MOD 4 = 1) ==> SING (fixes flip (mills p)) 507Proof 508 rpt strip_tac >> 509 `?x y. p = windmill x y y` by rw[fermat_two_squares_exists_windmill] >> 510 qabbrev_tac `s = mills p` >> 511 qabbrev_tac `b = fixes flip s` >> 512 qabbrev_tac `u = (x,y,y)` >> 513 `u IN s` by rw[mills_element, Abbr`u`, Abbr`s`] >> 514 `u IN b` by metis_tac[fixes_element, triple_parts, flip_fix] >> 515 simp[SING_DEF, EXTENSION] >> 516 qexists_tac `u` >> 517 (rewrite_tac[EQ_IMP_THM] >> simp[]) >> 518 rpt strip_tac >> 519 `x' IN s /\ ?h k. x' = (h,k,k)` by metis_tac[fixes_element, triple_parts, flip_fix] >> 520 `p = windmill h k k` by fs[mills_element, Abbr`s`] >> 521 qabbrev_tac `c = 2 * y` >> 522 qabbrev_tac `d = 2 * k` >> 523 `ODD p` by rfs[odd_by_mod_4] >> 524 `(p = x ** 2 + c ** 2) /\ ODD x` by metis_tac[windmill_x_y_y] >> 525 `(p = h ** 2 + d ** 2) /\ ODD h` by metis_tac[windmill_x_y_y] >> 526 `EVEN c /\ EVEN d` by rw[EVEN_DOUBLE, Abbr`c`, Abbr`d`] >> 527 `(x = h) /\ (c = d)` by metis_tac[fermat_two_squares_unique_odd_even] >> 528 `y = k` by fs[Abbr`c`, Abbr`d`] >> 529 simp[Abbr`u`] 530QED 531 532 533(* Theorem: prime p /\ (p MOD 4 = 1) ==> 534 (fixes flip (mills p) = 535 {(let u = (1,1,p DIV 4) ; 536 n = iterate_period (zagier o flip) u 537 in FUNPOW (zagier o flip) (HALF n) u)}) *) 538(* Proof: 539 Let s = mills p, 540 v = FUNPOW (zagier o flip) (HALF n) u, 541 a = fixes zagier s, 542 b = fixes flip s. 543 The goal is: b = {v}. 544 545 Note a = {u} by zagier_fixes_prime 546 so u IN a by IN_SING 547 Also ~square p by prime_non_square 548 so FINITE s by mills_non_square_finite 549 ==> zagier involute s by zagier_involute_mills, p MOD 4 = 1 550 and flip involute s by flip_involute_mills 551 ==> ODD n by involute_involute_fix_sing_period_odd 552 so v IN b by involute_involute_fix_orbit_fix_odd 553 Thus b = {v} by flip_fixes_prime_sing, SING_DEF, IN_SING 554*) 555Theorem flip_fixes_prime: 556 !p. prime p /\ (p MOD 4 = 1) ==> 557 let u = (1, 1, p DIV 4); 558 n = iterate_period (zagier o flip) u 559 in fixes flip (mills p) = {FUNPOW (zagier o flip) (n DIV 2) u} 560Proof 561 rw_tac bool_ss[] >> 562 qabbrev_tac `s = mills p` >> 563 qabbrev_tac `v = FUNPOW (zagier o flip) (HALF n) u` >> 564 qabbrev_tac `a = fixes zagier s` >> 565 qabbrev_tac `b = fixes flip s` >> 566 `a = {u}` by metis_tac[zagier_fixes_prime] >> 567 `u IN a` by fs[] >> 568 `~square p` by metis_tac[prime_non_square] >> 569 `FINITE s` by fs[mills_non_square_finite, Abbr`s`] >> 570 `zagier involute s` by metis_tac[zagier_involute_mills, ONE_NOT_ZERO] >> 571 `flip involute s` by metis_tac[flip_involute_mills] >> 572 drule_then strip_assume_tac involute_involute_fix_sing_period_odd >> 573 last_x_assum (qspecl_then [`zagier`, `flip`, `n`, `u`] strip_assume_tac) >> 574 `ODD n` by rfs[] >> 575 drule_then strip_assume_tac involute_involute_fix_orbit_fix_odd >> 576 last_x_assum (qspecl_then [`zagier`, `flip`, `n`, `u`] strip_assume_tac) >> 577 `v IN b` by rfs[] >> 578 metis_tac[flip_fixes_prime_sing, SING_DEF, IN_SING] 579QED 580 581(* Theorem: prime p /\ (p MOD 4 = 1) /\ (u = (1, 1, p DIV 4)) /\ 582 (n = iterate_period (zagier o flip) u) ==> 583 (fixes flip (mills p) = {FUNPOW (zagier o flip) (n DIV 2) u}) *) 584(* Proof: by flip_fixes_prime *) 585Theorem flip_fixes_prime_alt: 586 !p u n. prime p /\ (p MOD 4 = 1) /\ (u = (1, 1, p DIV 4)) /\ 587 (n = iterate_period (zagier o flip) u) ==> 588 (fixes flip (mills p) = {FUNPOW (zagier o flip) (n DIV 2) u}) 589Proof 590 metis_tac[flip_fixes_prime] 591QED 592 593 594(* Theorem: !p. prime p /\ (p MOD 4 = 1) ==> 595 ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) *) 596(* Proof: 597 Existence part by fermat_two_squares_exists_odd_even 598 Uniqueness part by fermat_two_squares_unique_odd_even 599*) 600Theorem fermat_two_squares_thm: 601 !p. prime p /\ (p MOD 4 = 1) ==> 602 ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) 603Proof 604 rw[Once EXISTS_UNIQUE_THM] >| [ 605 `?u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2)` by rw[fermat_two_squares_exists_odd_even] >> 606 qexists_tac `u` >> 607 rw[Once EXISTS_UNIQUE_THM], 608 metis_tac[fermat_two_squares_unique_odd_even] 609 ] 610QED 611 612 613(* Theorem: prime p ==> 614 ((p MOD 4 = 1) <=> ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2)) *) 615(* Proof: 616 If part: p MOD 4 = 1 ==> ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) 617 This is true by fermat_two_squares_thm 618 Only-if part: ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2) ==> p MOD 4 = 1 619 Note ?u v. ODD u /\ EVEN v /\ 620 (p = u ** 2 + v ** 2) by EXISTS_UNIQUE_THM 621 but ODD (u ** 2) by ODD_EXP 622 and EVEN (v ** 2) by EVEN_EXP 623 Thus ODD p by ODD_ADD, EVEN_ODD 624 so p MOD 4 = 1 or p MOD 4 = 3 by mod_4_odd 625 but p MOD 4 <> 3 by mod_4_not_squares 626 ==> p MOD 4 = 1 627*) 628Theorem fermat_two_squares_iff: 629 !p. prime p ==> ((p MOD 4 = 1) <=> ?!u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2)) 630Proof 631 rw[EQ_IMP_THM] >- 632 rw[fermat_two_squares_thm] >> 633 `?u v. ODD u /\ EVEN v /\ (p = u ** 2 + v ** 2)` by metis_tac[EXISTS_UNIQUE_THM] >> 634 `ODD (u ** 2)` by rw[ODD_EXP] >> 635 `EVEN (v ** 2)` by rw[EVEN_EXP] >> 636 `ODD p` by fs[ODD_ADD, EVEN_ODD] >> 637 `(p MOD 4 = 1) \/ (p MOD 4 = 3)` by fs[mod_4_odd] >> 638 fs[mod_4_not_squares] 639QED 640 641(* ------------------------------------------------------------------------- *) 642(* Fermat Two-Squares Algorithm. *) 643(* ------------------------------------------------------------------------- *) 644 645 646(* Theorem: (iterate_period (zagier o flip) (1,1,z) = 1) <=> (z = 1) *) 647(* Proof: 648 By iterate_period_thm, this is to show: 649 (1) zagier (flip (1,1,z)) = (1,1,z) ==> z = 1 650 This is true by zagier_flip_1_1_z 651 (2) zagier (flip (1,1,1)) = (1,1,1) 652 This is true by zagier_flip_1_1_z 653*) 654Theorem zagier_flip_1_1_z_period: 655 !z. (iterate_period (zagier o flip) (1,1,z) = 1) <=> (z = 1) 656Proof 657 rw[iterate_period_thm, zagier_flip_1_1_z] 658QED 659 660 661(* Theorem: prime p /\ (p MOD 4 = 1) /\ (u = (1, 1, p DIV 4)) /\ 662 (n = iterate_period (zagier o flip) u) /\ (g = \(x,y,z). y <> z) ==> 663 !j. j < n DIV 2 ==> g (FUNPOW (zagier o flip) j u) /\ 664 ~g (FUNPOW (zagier o flip) (n DIV 2) u) *) 665(* Proof: 666 Let s = mills p, 667 f = zagier o flip, 668 h = n DIV 2, 669 v = FUNPOW f h u. 670 The goal is: ~g v /\ !j. j < h ==> g (FUNPOW f j u) 671 672 Note ~square p by prime_non_square 673 so FINITE s by mills_non_square_finite 674 and zagier involute s by zagier_involute_mills 675 and flip involute s by flip_involute_mills 676 so f PERMUTES s by involute_involute_permutes 677 Also u IN s by mills_element_trivial 678 and 0 < n by iterate_period_pos 679 so h < n by HALF_LESS, 0 < n 680 681 Without assuming the following: 682 Note fixes flip s = {v} by flip_fixes_prime 683 But use: 684 Note fixes zagier s = {u} by zagier_fixes_prime 685 so ODD n by involute_involute_fix_sing_period_odd 686 Then use involute_involute_fix_odd_period_fix: 687 !j. 0 < j /\ j < n ==> 688 (FUNPOW (zagier o flip) j u IN fixes flip s <=> (j = h)) 689 690 Case: n = 1, 691 Then h = 0 by HALF_EQ_0, n <> 0 692 so v = u by FUNPOW_0 693 = (1,1,1) by zagier_flip_1_1_z_period 694 or ~g v by definition of g 695 and the other case needs j < 0, which is irrelevant. 696 697 Case: n <> 1. 698 and h <> 0, so 0 < h by HALF_EQ_0 699 700 Claim: ~g v 701 Note 0 < h < n, 702 so v IN fixes flip s by involute_involute_fix_odd_period_fix 703 Now ?x y z. v = (x,y,z) by triple_parts 704 and y = z by fixes_element, flip_fix 705 or ~g v by definition of g 706 707 Claim: !j. j < h ==> g (FUNPOW f j u) 708 By contradiction, suppose ?j. j < h /\ ~g (FUNPOW f j u). 709 Let w = FUNPOW f j u. 710 Note ?x y z. w = (x,y,z) by triple_parts 711 so y = z by ~g w 712 Now w IN s by FUNPOW_closure 713 and flip w = w by flip_fix 714 ==> w IN (fixes flip u} by fixes_element 715 If j = 0, w = u by FUNPOW_0 716 Then z = 1, so n = 1 by zagier_flip_1_1_z_period, PAIR_EQ 717 which contradicts n <> 1. 718 If j <> 0, 719 Then j = h by involute_involute_fix_odd_period_fix 720 which contradicts j < h. 721*) 722Theorem flip_fixes_iterates_prime: 723 !p u n g. prime p /\ (p MOD 4 = 1) /\ (u = (1, 1, p DIV 4)) /\ 724 (n = iterate_period (zagier o flip) u) /\ (g = \(x,y,z). y <> z) ==> 725 (~g (FUNPOW (zagier o flip) (n DIV 2) u) /\ 726 (!j. j < n DIV 2 ==> g (FUNPOW (zagier o flip) j u))) 727Proof 728 ntac 5 strip_tac >> 729 qabbrev_tac `f = zagier o flip` >> 730 qabbrev_tac `s = mills p` >> 731 qabbrev_tac `h = n DIV 2` >> 732 qabbrev_tac `v = FUNPOW f h u` >> 733 `~square p` by rw[prime_non_square] >> 734 `FINITE s` by metis_tac[mills_non_square_finite] >> 735 `zagier involute s` by metis_tac[zagier_involute_mills, ONE_NOT_ZERO] >> 736 `flip involute s` by metis_tac[flip_involute_mills] >> 737 `f PERMUTES s` by fs[involute_involute_permutes, Abbr`f`] >> 738 `u IN s` by rw[mills_element_trivial, Abbr`s` ] >> 739 `0 < n` by metis_tac[iterate_period_pos] >> 740 `h < n` by rw[HALF_LESS, Abbr`h`] >> 741 `fixes zagier s = {u}` by metis_tac[zagier_fixes_prime] >> 742 drule_then strip_assume_tac (involute_involute_fix_sing_period_odd |> ISPEC ``zagier``) >> 743 last_x_assum (qspecl_then [`flip`, `n`, `u`] strip_assume_tac) >> 744 `ODD n` by rfs[] >> 745 Cases_on `n = 1` >| [ 746 `h = 0` by rw[Abbr`h`] >> 747 `v = u` by rw[Abbr`v`] >> 748 `_ = (1,1,1)` by metis_tac[zagier_flip_1_1_z_period] >> 749 rfs[], 750 `h <> 0` by rw[HALF_EQ_0, Abbr`h`] >> 751 `0 < h` by decide_tac >> 752 drule_then strip_assume_tac (involute_involute_fix_odd_period_fix |> ISPEC ``f: num # num # num -> num # num # num``) >> 753 `!j. 0 < j /\ j < n ==> 754 (FUNPOW (zagier o flip) j u IN fixes flip s <=> (j = HALF n))` by rfs[] >> 755 `~g v` by 756 (`FUNPOW (zagier o flip) h u IN fixes flip s <=> (h = HALF n)` by rfs[] >> 757 `v IN fixes flip s` by metis_tac[] >> 758 `?x y z. v = (x,y,z)` by rw[triple_parts] >> 759 `y = z` by metis_tac[fixes_element, flip_fix] >> 760 fs[]) >> 761 rpt strip_tac >> 762 spose_not_then strip_assume_tac >> 763 qabbrev_tac `w = FUNPOW f j u` >> 764 `?x y z. w = (x,y,z)` by rw[triple_parts] >> 765 `~g w <=> (y = z)` by fs[] >> 766 `y = z` by fs[] >> 767 `w IN s` by rw[FUNPOW_closure, Abbr`w`] >> 768 `w IN (fixes flip s)` by fs[flip_fix, fixes_element] >> 769 Cases_on `j = 0` >| [ 770 `w = u` by fs[Abbr`w`] >> 771 metis_tac[zagier_flip_1_1_z_period, pairTheory.PAIR_EQ], 772 `0 < j` by decide_tac >> 773 `FUNPOW (zagier o flip) j u IN fixes flip s <=> (j = HALF n)` by rfs[] >> 774 `j = h` by metis_tac[] >> 775 decide_tac 776 ] 777 ] 778QED 779 780(* This proof is using: 781 zagier_flip_1_1_z_period and involute_involute_fix_odd_period_fix. 782 In part4, there is another proof using: 783 flip_fixes_prime, which depends on fermat_two_squares_unique_odd_even. 784*) 785 786(* ------------------------------------------------------------------------- *) 787(* Computation by WHILE loop. *) 788(* ------------------------------------------------------------------------- *) 789 790(* Define the exit condition *) 791val found_def = Define` 792 found (x:num, y:num, z:num) <=> (y = z) 793`; 794 795 796(* Theorem: $~ o found = (\(x,y,z). y <> z) *) 797(* Proof: by found_def, FUN_EQ_THM. *) 798Theorem found_not: 799 $~ o found = (\(x,y,z). y <> z) 800Proof 801 rw[found_def, pairTheory.FORALL_PROD, FUN_EQ_THM] 802QED 803 804(* Idea: use WHILE for search. Develop theory in iterateCompute. *) 805 806(* Compute two squares of Fermat's theorem by WHILE loop. *) 807val two_sq_def = Define` 808 two_sq n = WHILE ($~ o found) (zagier o flip) (1,1,n DIV 4) 809`; 810 811(* 812> EVAL ``two_sq 5``; = (1,2): thm (1,1,1) 813> EVAL ``two_sq 13``; = (3,2): thm (3,1,1) 814> EVAL ``two_sq 17``; = (1,4): thm (1,2,2) 815> EVAL ``two_sq 29``; = (5,2): thm (5,1,1) 816> EVAL ``two_sq 97``; = (9,4): thm (9,2,2) 817*) 818 819 820(* Theorem: two_sq n = WHILE (\(x,y,z). y <> z) (zagier o flip) (1, 1, n DIV 4) *) 821(* Proof: by two_sq_def, found_not. *) 822Theorem two_sq_alt: 823 !n. two_sq n = WHILE (\(x,y,z). y <> z) (zagier o flip) (1, 1, n DIV 4) 824Proof 825 simp[two_sq_def, found_not] 826QED 827 828(* Using direct WHILE, no need for Hoare specification. *) 829 830(* Theorem: prime p /\ (p MOD 4 = 1) ==> two_sq p IN fixes flip (mills p) *) 831(* Proof: 832 Let s = mills p, 833 f = zagier o flip, 834 u = (1,1,p DIV 4), 835 n = iterate_period f u. 836 By two_sq_def, this is to show: WHILE ($~ o found) f u IN fixes flip s 837 838 Note (~) o found = (\t. flip t <> t) by found_def, flip_def, FUN_EQ_THM 839 Thus the goal is: WHILE (\t. flip t <> t) f u IN fixes flip s 840 841 Note ~square p by prime_non_square 842 so FINITE s by mills_non_square_finite 843 and zagier involute s by zagier_involute_mills, ~square p, p MOD 4 <> 0 844 and flip involute s by flip_involute_mills 845 also fixes zagier s = {u} by zagier_fixes_prime 846 so u IN fixes zagier s by IN_SING 847 and ODD n by involute_involute_fix_sing_period_odd 848 ==> WHILE (\y. flip y <> y) f u IN fixes flip s 849 by involute_involute_fix_odd_period_fix_while 850*) 851Theorem two_sq_thm: 852 !p. prime p /\ (p MOD 4 = 1) ==> two_sq p IN fixes flip (mills p) 853Proof 854 rw[two_sq_def] >> 855 qabbrev_tac `s = mills p` >> 856 qabbrev_tac `f = zagier o flip` >> 857 qabbrev_tac `u = (1,1,p DIV 4)` >> 858 qabbrev_tac `n = iterate_period f u` >> 859 `(~) o found = \t. flip t <> t` by 860 (rw[FUN_EQ_THM] >> 861 `?x y z. t = (x,y,z)` by rw[triple_parts] >> 862 rw[found_def, flip_def]) >> 863 simp[] >> 864 assume_tac (involute_involute_fix_odd_period_fix_while |> ISPEC ``zagier``) >> 865 last_x_assum (qspecl_then [`flip`, `s`, `n`, `u`] strip_assume_tac) >> 866 `~square p` by rw[prime_non_square] >> 867 `FINITE s` by rw[mills_non_square_finite, Abbr`s`] >> 868 `zagier involute s` by metis_tac[zagier_involute_mills, ONE_NOT_ZERO] >> 869 `flip involute s` by metis_tac[flip_involute_mills] >> 870 `fixes zagier s = {u}` by rw[zagier_fixes_prime, Abbr`s`, Abbr`u`] >> 871 drule_then strip_assume_tac involute_involute_fix_sing_period_odd >> 872 last_x_assum (qspecl_then [`zagier`, `flip`, `n`, `u`] strip_assume_tac) >> 873 rfs[Abbr`f`] 874QED 875 876(* Very good -- nice and simple! *) 877 878 879(* Theorem: prime p /\ (p MOD 4 = 1) ==> 880 HOARE_SPEC (fixes zagier (mills p)) 881 (WHILE ($~ o found) (zagier o flip)) 882 (fixes flip (mills p)) *) 883(* Proof: 884 Let s = mills p, 885 f = zagier o flip, 886 u = (1,1,p DIV 4), 887 n = iterate_period f u, 888 h = n DIV 2, 889 v = FUNPOW f h u, 890 g = $~ o found, use ~g to exit WHILE loop, 891 a = fixes zagier s, 892 b = fixes flip s. 893 The goal is: HOARE_SPEC a (WHILE g f) b 894 895 Note a = {u} by zagier_fixes_prime 896 and b = {v} by flip_fixes_prime 897 898 If n = 1, that is, period = 1 for prime p = 5. 899 then u = (1,1,1) by zagier_flip_1_1_z_period, triple_parts 900 so ~g u by definition of g, found_def 901 but u IN s by mills_element_trivial 902 so u IN b by fixes_element, flip_fix 903 ==> b = a by EQUAL_SING, IN_SING 904 Thus HOARE_SPEC {u} (WHILE g f) {u} 905 by iterate_while_hoare_0 906 or HOARE_SPEC a (WHILE g f) b. 907 908 If n <> 1, 909 Note ~square p by prime_non_square 910 so FINITE s by mills_non_square_finite 911 and zagier involute s by zagier_involute_mills, ~square p 912 and flip involute s by flip_involute_mills 913 so f PERMUTES s by involute_involute_permutes 914 Also 0 < n by iterate_period_pos, u IN s 915 and ODD n by involute_involute_fix_sing_period_odd 916 so n <> 2 by EVEN_2, ODD_EVEN 917 ==> 1 + h < n by half_add1_lt, 2 < n 918 or h < n - 1 by arithmetic 919 920 Now ~g v /\ (!j. j < h ==> g (FUNPOW f j u)) 921 by flip_fixes_iterates_prime, found_not 922 Thus HOARE_SPEC a (WHILE g f) b by iterate_while_hoare 923*) 924Theorem two_sq_while_hoare: 925 !p. prime p /\ (p MOD 4 = 1) ==> 926 HOARE_SPEC (fixes zagier (mills p)) 927 (WHILE ($~ o found) (zagier o flip)) 928 (fixes flip (mills p)) 929Proof 930 rpt strip_tac >> 931 qabbrev_tac `s = mills p` >> 932 qabbrev_tac `f = zagier o flip` >> 933 qabbrev_tac `u = (1,1,p DIV 4)` >> 934 qabbrev_tac `n = iterate_period f u` >> 935 qabbrev_tac `h = n DIV 2` >> 936 qabbrev_tac `v = FUNPOW f h u` >> 937 qabbrev_tac `g = $~ o found` >> 938 qabbrev_tac `a = fixes zagier s` >> 939 qabbrev_tac `b = fixes flip s` >> 940 `u IN s` by rw[mills_element_trivial, Abbr`u`, Abbr`s`] >> 941 `a = {u}` by rw[zagier_fixes_prime, Abbr`a`, Abbr`u`, Abbr`s`] >> 942 `b = {v}` by metis_tac[flip_fixes_prime, SING_DEF] >> 943 Cases_on `n = 1` >| [ 944 `u = (1,1,1)` by metis_tac[zagier_flip_1_1_z_period, triple_parts] >> 945 `~g u` by fs[found_def, Abbr`g`] >> 946 `u IN s` by rw[mills_element_trivial, Abbr`u`, Abbr`s`] >> 947 `u IN b` by rw[fixes_element, flip_fix, Abbr`b`] >> 948 `b = a` by metis_tac[EQUAL_SING, IN_SING] >> 949 rw[iterate_while_hoare_0], 950 `~square p` by rw[prime_non_square] >> 951 `FINITE s` by metis_tac[mills_non_square_finite] >> 952 `zagier involute s` by metis_tac[zagier_involute_mills, ONE_NOT_ZERO] >> 953 `flip involute s` by metis_tac[flip_involute_mills] >> 954 `f PERMUTES s` by fs[involute_involute_permutes, Abbr`f`] >> 955 `0 < n` by metis_tac[iterate_period_pos] >> 956 drule_then strip_assume_tac involute_involute_fix_sing_period_odd >> 957 last_x_assum (qspecl_then [`zagier`, `flip`, `n`, `u`] strip_assume_tac) >> 958 `ODD n` by rfs[] >> 959 `n <> 2` by metis_tac[EVEN_2, ODD_EVEN] >> 960 `1 + h < n` by rw[half_add1_lt, Abbr`h`] >> 961 `h < n - 1` by decide_tac >> 962 drule_then strip_assume_tac iterate_while_hoare >> 963 last_x_assum (qspecl_then [`u`, `f`, `n-1`, `n`, `g`, `h`] strip_assume_tac) >> 964 `~g v /\ (!j. j < h ==> g (FUNPOW f j u))` by metis_tac[found_not, flip_fixes_iterates_prime] >> 965 rfs[] 966 ] 967QED 968 969 970(* Theorem: prime p /\ (p MOD 4 = 1) ==> 971 let (x,y,z) = WHILE ($~ o found) (zagier o flip) (1,1,p DIV 4) 972 in p = x ** 2 + (y + z) ** 2 *) 973(* Proof: 974 Let s = mills p, 975 u = (1,1,p DIV 4), 976 a = fixes zagier s, 977 b = fixes flip s. 978 979 Note HOARE_SPEC a 980 (WHILE ($~ o found) (zagier o flip)) b 981 by two_sq_while_hoare 982 and a = {u} by zagier_fixes_prime 983 984 By HOARE_SPEC_DEF, this is to show: 985 !s. a s ==> b (WHILE ($~ o found) (zagier o flip) s) 986 Thus s = u by IN_SING, function as set 987 Let w = WHILE ($~ o found) (zagier o flip) u. 988 Then ?x y z. w = (x,y,z) by triple_parts 989 and w IN b by b w, function as set 990 so w IN s /\ y = z by fixes_element, flip_fix 991 Thus p 992 = windmill x y z by mills_element, w IN s 993 = windmill x y y by y = z 994 = x ** 2 + (2 * y) ** 2 by windmill_by_squares 995 = x ** 2 + (y + z) ** 2 by arithmetic, y = z 996*) 997Theorem two_sq_while_thm: 998 !p. prime p /\ (p MOD 4 = 1) ==> 999 let (x,y,z) = WHILE ($~ o found) (zagier o flip) (1,1,p DIV 4) 1000 in p = x ** 2 + (y + z) ** 2 1001Proof 1002 rw_tac bool_ss[] >> 1003 drule_then strip_assume_tac two_sq_while_hoare >> 1004 rfs[whileTheory.HOARE_SPEC_DEF] >> 1005 qabbrev_tac `u = (1,1,p DIV 4)` >> 1006 `fixes zagier (mills p) = {u}` by rw[zagier_fixes_prime, Abbr`u`] >> 1007 last_x_assum (qspecl_then [`u`] strip_assume_tac) >> 1008 `(x,y,z) IN fixes flip (mills p)` by rfs[SPECIFICATION] >> 1009 `(x,y,z) IN (mills p) /\ (y = z)` by fs[fixes_element, flip_fix] >> 1010 metis_tac[mills_element, windmill_by_squares, DECIDE``y + y = 2 * y``] 1011QED 1012 1013(* A beautiful theorem! *) 1014 1015(* Define the algorithm. *) 1016val two_squares_def = Define` 1017 two_squares n = let (x,y,z) = two_sq n in (x, y + z) 1018`; 1019 1020(* 1021> EVAL ``two_squares 5``; = (1,2) 1022> EVAL ``two_squares 13``; = (3,2) 1023> EVAL ``two_squares 17``; = (1,4) 1024> EVAL ``two_squares 29``; = (5,2) 1025> EVAL ``two_squares 97``; = (9,4) 1026> EVAL ``MAP two_squares [5; 13; 17; 29; 37; 41; 53; 61; 73; 89; 97]``; 1027= [(1,2); (3,2); (1,4); (5,2); (1,6); (5,4); (7,2); (5,6); (3,8); (5,8); (9,4)]: thm 1028*) 1029 1030 1031(* Theorem: prime p /\ (p MOD 4 = 1) ==> 1032 let (u,v) = two_squares p in (p = u ** 2 + v ** 2) *) 1033(* Proof: by two_squares_def, two_sq_def, two_sq_while_thm. *) 1034Theorem two_squares_thm: 1035 !p. prime p /\ (p MOD 4 = 1) ==> 1036 let (u,v) = two_squares p in (p = u ** 2 + v ** 2) 1037Proof 1038 rw[two_squares_def, two_sq_def] >> 1039 drule_then strip_assume_tac two_sq_while_thm >> 1040 qabbrev_tac `loop = WHILE ($~ o found) (zagier o flip) (1,1,p DIV 4)` >> 1041 `?x y z. loop = (x,y,z)` by fs[triple_parts] >> 1042 fs[] 1043QED 1044 1045(* Another proof of the same theorem, using two_sq_thm. *) 1046 1047(* Theorem: prime p /\ (p MOD 4 = 1) ==> 1048 let (u,v) = two_squares p in (p = u ** 2 + v ** 2) *) 1049(* Proof: 1050 Let t = two_sq p. 1051 Then t IN fixes flip (mills p) by two_sq_thm 1052 and ?x y z. t = (x,y,z) by triple_parts 1053 ==> (x,y,z) IN mills p /\ (y = z) by fixes_element, flip_fix 1054 so p = windmill x y y by mills_element 1055 = x ** 2 + (2 * y) ** 2 by windmill_by_squares 1056 = x ** 2 + (y + z) ** 2 by y = z 1057 = u ** 2 + v ** 2 by two_squares_def 1058*) 1059Theorem two_squares_thm: 1060 !p. prime p /\ (p MOD 4 = 1) ==> 1061 let (u,v) = two_squares p in (p = u ** 2 + v ** 2) 1062Proof 1063 rw[two_squares_def] >> 1064 qabbrev_tac `t = two_sq p` >> 1065 `t IN fixes flip (mills p)` by rw[two_sq_thm, Abbr`t`] >> 1066 `?x y z. t = (x,y,z)` by rw[triple_parts] >> 1067 `(x,y,z) IN mills p /\ (y = z)` by fs[fixes_element, flip_fix] >> 1068 `p = windmill x y y` by fs[mills_element] >> 1069 simp[windmill_by_squares] 1070QED 1071 1072(* ------------------------------------------------------------------------- *) 1073(* Fermat's Two Squares Theorem by Group Action. *) 1074(* ------------------------------------------------------------------------- *) 1075 1076(* ------------------------------------------------------------------------- *) 1077(* Relation to Fixes and Pairs. *) 1078(* ------------------------------------------------------------------------- *) 1079 1080(* Theorem: f involute X ==> fixed_points (FUNPOW f) Z2 X = fixes f X *) 1081(* Proof: 1082 By fixed_points_def, fixes_def, EXTENSION, this is to show: 1083 (1) !a. a < 2 ==> (FUNPOW f a x = x) ==> f x = x 1084 Note f x = FUNPOW f 1 x = x by FUNPOW_1, 1 < 2 1085 (2) f x = x ==> !a. a < 2 ==> FUNPOW f a x = x 1086 When a = 0, FUNPOW f 0 x = x by FUNPOW_0 1087 When a = 1, FUNPOW f 1 x = f x = x by FUNPOW_1, f x = x 1088*) 1089Theorem involute_fixed_points_eq_fixes: 1090 !f X. f involute X ==> fixed_points (FUNPOW f) Z2 X = fixes f X 1091Proof 1092 rw[fixed_points_def, fixes_def, EXTENSION, EQ_IMP_THM] >- 1093 metis_tac[FUNPOW_1, DECIDE``1 < 2``] >> 1094 (`a = 0 \/ a = 1` by decide_tac >> simp[]) 1095QED 1096 1097(* Theorem: f involute X ==> BIGUNION (multi_orbits (FUNPOW f) Z2 X) = pairs f X *) 1098(* Proof: 1099 By multi_orbits_def, pairs_def, BIGUNION, EXTENSION, this is to show: 1100 (1) x IN e /\ e IN orbits (FUNPOW f) Z2 X /\ ~SING e ==> x IN X /\ f x <> x 1101 Note x IN X by involute_orbits_element_element 1102 and e = orbit (FUNPOW f) Z2 x by involute_orbits_element_is_orbit 1103 By contradiction, suppose f x = x, 1104 Then e = {x} by involute_orbit_fixed 1105 with contradicts ~SING e by SING 1106 (2) x IN X /\ f x <> x ==> ?e. x IN e /\ e IN orbits (FUNPOW f) Z2 X /\ ~SING e 1107 Let e = {x; f x}. 1108 Then x IN e, and ~SING e by f x <> x 1109 The goal is to show: {x; f x} IN orbits (FUNPOW f) Z2 X 1110 1111 Note {x; f x} 1112 = orbit (FUNPOW f) Z2 x by involute_orbit_not_fixed 1113 which is IN orbits (FUNPOW f) Z2 X by funpow_orbit_in_orbits 1114*) 1115Theorem involute_multi_orbits_union_eq_pairs: 1116 !f X. f involute X ==> BIGUNION (multi_orbits (FUNPOW f) Z2 X) = pairs f X 1117Proof 1118 rw[multi_orbits_def, pairs_def, EXTENSION, EQ_IMP_THM] >- 1119 metis_tac[involute_orbits_element_element] >- 1120 (`x IN X` by metis_tac[involute_orbits_element_element] >> 1121 `s = orbit (FUNPOW f) Z2 x` by metis_tac[involute_orbits_element_is_orbit] >> 1122 metis_tac[involute_orbit_fixed, SING]) >> 1123 qexists_tac `{x; f x}` >> 1124 simp[] >> 1125 `{x; f x} = orbit (FUNPOW f) Z2 x` by metis_tac[involute_orbit_not_fixed] >> 1126 rw[funpow_orbit_in_orbits] 1127QED 1128 1129(* ------------------------------------------------------------------------- *) 1130(* Fermat's Two-Squares Theorem (Existence) *) 1131(* ------------------------------------------------------------------------- *) 1132 1133(* With zagier and flip both involutions on (mills p), 1134 and zagier with a unique fixed point, 1135 this implies flip has at least a fixed point, 1136 giving the existence of Fermat's two squares. 1137*) 1138 1139(* Proof based on involute_two_fixed_points_both_odd *) 1140 1141(* Theorem: prime p /\ p MOD 4 = 1 ==> 1142 fixed_points (FUNPOW zagier) Z2 (mills p) = {(1,1,p DIV 4)} *) 1143(* Proof: 1144 By fixes_def, mills_def, this is to show: 1145 (1) p = windmill x y z /\ zagier (x,y,z) = (x,y,z) ==> (x,y,z) = (1,1,p DIV 4) 1146 Note ~square p by prime_non_square 1147 and x <> 0 by mills_element, mills_triple_nonzero 1148 so x = y by zagier_fix 1149 ==> (x,y,z) = (1,1,p DIV 4) by windmill_trivial_prime 1150 (2) p MOD 4 = 1 ==> p = windmill 1 1 (p DIV 4) 1151 This is true by windmill_trivial_prime 1152 (3) a < 2 ==> FUNPOW zagier p (1,1,p DIV 4) = (1,1,p DIV 4) 1153 When a = 0, true by FUNPOW_0 1154 When a = 1, true by FUNPOW_1, zagier_fix 1155*) 1156Theorem zagier_fixed_points_sing: 1157 !p. prime p /\ p MOD 4 = 1 ==> 1158 fixed_points (FUNPOW zagier) Z2 (mills p) = {(1,1,p DIV 4)} 1159Proof 1160 rw[fixed_points_def, mills_def, Zadd_element, EXTENSION] >> 1161 simp[EQ_IMP_THM] >> 1162 rpt strip_tac >| [ 1163 rename1 ���x = (x',y,z)��� >> 1164 `~square p` by metis_tac[prime_non_square] >> 1165 `p MOD 4 <> 0` by decide_tac >> 1166 `zagier x = x` by metis_tac[FUNPOW_1, DECIDE``1 < 2``] >> 1167 `x' = y` by metis_tac[zagier_fix, mills_element, mills_triple_nonzero] >> 1168 metis_tac[windmill_trivial_prime], 1169 metis_tac[windmill_trivial_prime], 1170 (`a = 0 \/ a = 1` by decide_tac >> fs[zagier_fix]) 1171 ] 1172QED 1173 1174(* A better proof of the same theorem *) 1175 1176(* Theorem: prime p /\ p MOD 4 = 1 ==> 1177 fixed_points (FUNPOW zagier) Z2 (mills p) = {(1,1,p DIV 4)} *) 1178(* Proof: 1179 Let X = mills p, a = fixed_points (FUNPOW zagier) Z2 X. 1180 The goal is: a = {(1,1,p DIV 4)}. 1181 1182 Note ~square p by prime_non_square 1183 so zagier involute X by zagier_involute_mills 1184 1185 By EXTENSION, this is to show: 1186 (1) t IN X ==> t = (1,1,p DIV 4) 1187 Note t IN X by involute_fixed_points_element_element 1188 and zagier t = t by involute_fixed_points 1189 Now ?x y z. t = (x, y, z) by triple_parts 1190 and x <> 0 by mills_triple_nonzero, ~square p 1191 ==> x = y by zagier_fix 1192 Note p = windmill x y z by mills_element 1193 ==> (x,y,z) = (1,1,p DIV 4) by windmill_trivial_prime 1194 1195 (2) (1,1,p DIV 4) IN a 1196 Note (1,1,p DIV 4) IN X by mills_element_trivial 1197 and zagier (1,1,p DIV 4) 1198 = zagier (1,1,p DIV 4) by zagier_1_1_z 1199 so (1,1,p DIV 4) IN a by involute_fixed_points_iff 1200*) 1201Theorem zagier_fixed_points_sing: 1202 !p. prime p /\ p MOD 4 = 1 ==> 1203 fixed_points (FUNPOW zagier) Z2 (mills p) = {(1,1,p DIV 4)} 1204Proof 1205 rpt strip_tac >> 1206 qabbrev_tac `X = mills p` >> 1207 qabbrev_tac `a = fixed_points (FUNPOW zagier) Z2 X` >> 1208 `~square p` by metis_tac[prime_non_square] >> 1209 `p MOD 4 <> 0` by decide_tac >> 1210 `zagier involute X` by metis_tac[zagier_involute_mills] >> 1211 rw[EXTENSION, EQ_IMP_THM] >| [ 1212 `x IN X` by metis_tac[involute_fixed_points_element_element] >> 1213 `zagier x = x` by metis_tac[involute_fixed_points] >> 1214 `?x1 y z. x = (x1, y, z)` by rw[triple_parts] >> 1215 `x1 <> 0` by metis_tac[mills_triple_nonzero] >> 1216 `x1 = y` by fs[zagier_fix] >> 1217 metis_tac[windmill_trivial_prime, mills_element], 1218 `(1,1,p DIV 4) IN X` by rw[mills_element_trivial, Abbr`X`] >> 1219 metis_tac[zagier_1_1_z, involute_fixed_points_iff] 1220 ] 1221QED 1222 1223(* Theorem: prime p /\ p MOD 4 = 1 ==> ?u v. p = u ** 2 + v ** 2 *) 1224(* Proof: 1225 Let X = mills p, the solutions (x,y,z) of p = x ** 2 + 4 * y * z. 1226 Note ~square p by prime_non_square 1227 so FINITE X by mills_non_square_finite 1228 Now !x y z. (x,y,z) IN X ==> 1229 x <> 0 /\ y <> 0 /\ z <> 0 by mills_triple_nonzero 1230 and zagier involute X by zagier_involute_mills 1231 and flip involute X by flip_involute_mills 1232 1233 Let a = fixed_points (FUNPOW zagier) Z2 X, 1234 b = fixed_points (FUNPOW flip) Z2 X. 1235 Then ODD (CARD a) <=> ODD (CARD b) by involute_two_fixed_points_both_odd 1236 1237 The punch line: 1238 Let k = p DIV 4. 1239 Then a = {(1,1,k)} by zagier_fixed_points_sing 1240 Thus CARD a = 1 by CARD_SING 1241 so ODD (CARD b) by ODD_1, above 1242 ==> CARD b <> 0 by EVEN_0 1243 or b <> EMPTY by CARD_EMPTY 1244 thus ?x y z. (x,y,z) IN b by MEMBER_NOT_EMPTY, triple_parts 1245 so flip (x, y, z) = (x, y, z) by involute_fixed_points 1246 ==> y = z by flip_fix 1247 Note (x,y,z) IN X by fixed_points_element 1248 Put u = x, v = 2 * y. 1249 Then p = u ** 2 + v ** 2 by mills_element, windmill_by_squares 1250*) 1251Theorem fermat_two_squares_exists_alt: 1252 !p. prime p /\ p MOD 4 = 1 ==> ?u v. p = u ** 2 + v ** 2 1253Proof 1254 rpt strip_tac >> 1255 qabbrev_tac `X = mills p` >> 1256 `~square p` by metis_tac[prime_non_square] >> 1257 `FINITE X` by fs[mills_non_square_finite, Abbr`X`] >> 1258 `zagier involute X` by metis_tac[zagier_involute_mills, DECIDE``1 <> 0``] >> 1259 `flip involute X` by metis_tac[flip_involute_mills] >> 1260 qabbrev_tac `a = fixed_points (FUNPOW zagier) Z2 X` >> 1261 qabbrev_tac `b = fixed_points (FUNPOW flip) Z2 X` >> 1262 `ODD (CARD a) <=> ODD (CARD b)` by rw[involute_two_fixed_points_both_odd, Abbr`a`, Abbr`b`] >> 1263 qabbrev_tac `k = p DIV 4` >> 1264 `a = {(1,1,k)}` by rw[zagier_fixed_points_sing, Abbr`a`, Abbr`X`, Abbr`k`] >> 1265 `CARD a = 1` by rw[] >> 1266 `CARD b <> 0` by metis_tac[ODD_1, EVEN_0, ODD_EVEN] >> 1267 `b <> EMPTY` by metis_tac[CARD_EMPTY] >> 1268 `?x y z. (x,y,z) IN b` by metis_tac[MEMBER_NOT_EMPTY, triple_parts] >> 1269 `flip (x, y, z) = (x, y, z)` by metis_tac[involute_fixed_points] >> 1270 `y = z` by fs[flip_fix] >> 1271 `(x,y,z) IN X` by fs[fixed_points_element, Abbr`b`] >> 1272 metis_tac[mills_element, windmill_by_squares] 1273QED 1274 1275(* Very good! *) 1276 1277(* ------------------------------------------------------------------------- *) 1278(* References *) 1279(* ------------------------------------------------------------------------- *) 1280 1281(* 1282 1283Biscuits of Number Theory (Arthur T. BenjamIN snd Ezra Brown editors) 1284ISBN: 9780883853405 1285Part IV: Sum of Squares and Polygonal Numbers 1286p.143 1287A One-Sentence Proof that Every Prime p = 1 (mod 4) is a Sum of Two Squares 1288Don Zagier 1289 1290 1291An Algorithm to List All the Fixed-Point Free Involutions on a Finite Set 1292Cyril Prissette (20 Jun 2010) 1293https://hal.archives-ouvertes.fr/hal-00493605/document 1294 1295The One-Sentence Proof in Multiple Sentences 1296Marcel Moosbrugger (Feb 3, 2020) 1297https://medium.com/cantors-paradise/the-one-sentence-proof-in-multiple-sentences-ab2657efc576 1298 1299*) 1300 1301 1302(* ------------------------------------------------------------------------- *) 1303 1304(* export theory at end *) 1305val _ = export_theory(); 1306 1307(*===========================================================================*) 1308