1(* ------------------------------------------------------------------------- *) 2(* Windmills and Involutions. *) 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 "windmill"; 12 13(* ------------------------------------------------------------------------- *) 14 15 16(* open dependent theories *) 17(* arithmeticTheory -- load by default *) 18(* val _ = load "helperTheory"; *) 19open helperTheory; 20open helperNumTheory; 21open helperSetTheory; 22open helperFunctionTheory; 23open arithmeticTheory pred_setTheory; 24(* val _ = load "dividesTheory"; *) 25open dividesTheory; 26(* val _ = load "gcdTheory"; *) 27open gcdTheory; (* for P_EUCLIDES *) 28 29open pairTheory; (* for FORALL_PROD, PAIR_EQ *) 30 31 32(* ------------------------------------------------------------------------- *) 33(* Windmills and Involutions Documentation *) 34(* ------------------------------------------------------------------------- *) 35(* 36 37 Helper Theorem: 38 triple_parts |- !t. ?x y z. t = (x,y,z) 39 40 Windmill Theory: 41 windmill_def |- !x y z. windmill x y z = x ** 2 + 4 * y * z 42 windmill_eq_0 |- !x y z. windmill x y z = 0 <=> x = 0 /\ y * z = 0 43 windmill_comm |- !x y z. windmill x y z = windmill x z y 44 windmill_trivial |- !k. windmill 1 1 k = 4 * k + 1 45 windmill_by_squares |- !x y. windmill x y y = x ** 2 + (2 * y) ** 2 46 windmill_x_y_y |- !n x y. n = windmill x y y ==> 47 n = x ** 2 + (2 * y) ** 2 /\ (ODD n <=> ODD x) 48 windmill_trivial_prime 49 |- !p. prime p /\ p MOD 4 = 1 ==> 50 !x z. p = windmill x x z <=> x = 1 /\ z = p DIV 4 51 52 Mills Theory: 53 mills_def |- !n. mills n = {(x,y,z) | n = windmill x y z} 54 mills_element |- !n x y z. (x,y,z) IN mills n <=> n = windmill x y z 55 mills_element_flip |- !n x y z. (x,y,z) IN mills n ==> (x,z,y) IN mills n 56 mills_element_triple |- !n t. t IN mills n <=> 57 ?x y z. t = (x,y,z) /\ n = windmill x y z 58 mills_element_trivial |- !n. n MOD 4 = 1 ==> (1,1,n DIV 4) IN mills n 59 mills_0 |- mills 0 = {(0,0,n) | n IN univ(:num)} UNION 60 {(0,n,0) | n IN univ(:num)} 61 mills_0_infinite |- INFINITE (mills 0) 62 mills_1 |- mills 1 = {1} CROSS ({0} CROSS univ(:num)) UNION 63 {1} CROSS (univ(:num) CROSS {0}) 64 mills_1_infinite |- INFINITE (mills 1) 65 mills_non_square_bound |- !n x y z. ~square n /\ 66 (x,y,z) IN mills n ==> x < n /\ y < n /\ z < n 67 mills_non_square_subset |- !n. ~square n ==> 68 mills n SUBSET count n CROSS (count n CROSS count n) 69 mills_non_square_card_upper 70 |- !n. ~square n ==> CARD (mills n) < n ** 3 71 mills_non_square_finite |- !n. ~square n ==> FINITE (mills n) 72 mills_finite_non_square |- !n. FINITE (mills n) <=> ~square n 73 74 mills_with_no_mind |- !n. (?y z. (0,y,z) IN mills n) <=> n MOD 4 = 0 75 mills_with_all_mind |- !n. n MOD 4 <> 0 <=> !x y z. (x,y,z) IN mills n ==> x <> 0 76 mills_with_no_arms |- !n. (?x z. (x,0,z) IN mills n) \/ 77 (?x y. (x,y,0) IN mills n) <=> square n 78 mills_with_all_arms |- !n. ~square n <=> !x y z. (x,y,z) IN mills n ==> y <> 0 /\ z <> 0 79 mills_quad_suc_non_empty|- !n. n MOD 4 = 1 ==> mills n <> {} 80 mills_triple_nonzero |- !n. ~square n /\ n MOD 4 <> 0 ==> 81 !x y z. (x,y,z) IN mills n ==> x <> 0 /\ y <> 0 /\ z <> 0 82 mills_prime_triple_nonzero 83 |- !p x y z. prime p /\ (x,y,z) IN mills p ==> 84 x <> 0 /\ y <> 0 /\ z <> 0 85 86 Flip involution: 87 flip_def |- !x y z. flip (x,y,z) = (x,z,y) 88 flip_fix |- !x y z. flip (x,y,z) = (x,y,z) <=> y = z 89 flip_closure |- !n x y z. (x,y,z) IN mills n ==> flip (x,y,z) IN mills n 90 flip_closure_alt |- !n t. t IN mills n ==> flip t IN mills n 91 flip_involute |- !x y z. flip (flip (x,y,z)) = (x,y,z) 92 flip_involute_alt |- !t. flip (flip t) = t 93 flip_involute_mills |- !n. flip involute (mills n) 94 95 Zagier's involution: 96 zagier_def |- !x y z. zagier (x,y,z) = 97 if x < y - z then (x + 2 * z,z,y - z - x) 98 else if x < 2 * y then (2 * y - x,y,x + z - y) 99 else (x - 2 * y,x + z - y,y) 100 zagier_0_y_z |- !y z. zagier (0,y,z) = 101 if z < y then (2 * z,z,y - z) 102 else if 0 < y then (2 * y,y,z - y) 103 else (0,z,0) 104 zagier_x_0_z |- !x z. zagier (x,0,z) = (x,x + z,0) 105 zagier_x_y_0 |- !x y. zagier (x,y,0) = 106 if x < y then (x,0,y - x) 107 else if x < 2 * y then (2 * y - x,y,x - y) 108 else (x - 2 * y,x - y,y) 109 zagier_1_y_1 |- !y. zagier (1,y,1) = 110 if y = 0 then (1,2,0) 111 else if y = 1 then (1,1,1) 112 else if y = 2 then (3,2,0) 113 else (3,1,y - 2) 114 zagier_1_1_z |- !z. zagier (1,1,z) = (1,1,z) 115 zagier_x_0_0 |- !x. zagier (x,0,0) = (x,x,0) 116 zagier_0_y_0 |- !y. zagier (0,y,0) = (0,0,y) 117 zagier_0_0_z |- !z. zagier (0,0,z) = (0,z,0) 118 zagier_0_0_0 |- zagier (0,0,0) = (0,0,0) 119 zagier_x_y_y |- !x y. zagier (x,y,y) = 120 if x < 2 * y then (2 * y - x,y,x) else (x - 2 * y,x,y) 121 zagier_fix |- !x y z. x <> 0 ==> (zagier (x,y,z) = (x,y,z) <=> x = y) 122 zagier_x_x_z |- !x z. x <> 0 ==> zagier (x,x,z) = (x,x,z) 123 zagier_closure |- !n x y z. (x,y,z) IN mills n ==> zagier (x,y,z) IN mills n 124 zagier_closure_alt |- !n t. t IN mills n ==> zagier t IN mills n 125 zagier_involute |- !x y z. x <> 0 /\ y <> 0 /\ z <> 0 ==> 126 zagier (zagier (x,y,z)) = (x,y,z) 127 zagier_involute_xyz |- !x y z. x * y * z <> 0 ==> zagier (zagier (x,y,z)) = (x,y,z) 128 zagier_involute_thm |- !t. FST t <> 0 /\ FST (SND t) <> 0 /\ SND (SND t) <> 0 ==> 129 zagier (zagier t) = t 130 cuboid_def |- !x y z. cuboid (x,y,z) = x * y * z 131 cuboid_eq_0 |- !x y z. (cuboid (x,y,z) = 0) <=> (x = 0) \/ (y = 0) \/ (z = 0) 132 zagier_involute_alt |- !t. cuboid t <> 0 ==> zagier (zagier t) = t 133 zagier_involute_mills 134 |- !n. ~square n /\ n MOD 4 <> 0 ==> zagier involute (mills n) 135 136 Mind of a windmill: 137 mind_def |- !x y z. mind (x,y,z) = 138 if x < y - z then x + 2 * z 139 else if x < y then 2 * y - x 140 else x 141 mind_zagier_eqn |- !x y z. mind (zagier (x,y,z)) = mind (x,y,z) 142 mind_zagier_thm |- !t. mind (zagier t) = mind t 143 mind_flip_eqn |- !x y z. mind (flip (x,y,z)) = 144 if x < z - y then x + 2 * y 145 else if x < z then 2 * z - x 146 else x 147 mind_flip_1_1_z |- !z. mind (flip (1,1,z)) = if z < 2 then 1 else 3 148 mind_flip_x_x_z |- !x z. mind (flip (x,x,z)) = 149 if x < z - x then 3 * x else if x < z then 2 * z - x else x 150 mind_flip_x_y_y |- !x y. mind (flip (x,y,y)) = if x < y then 2 * y - x else x 151 152 Gap of a Windmill: 153 gap_def |- !x y z. gap (x,y,z) = if y < z then z - y else y - z 154 gap_flip_eqn |- !x y z. gap (flip (x,y,z)) = gap (x,y,z) 155 gap_flip_thm |- !t. gap (flip t) = gap t 156 157 Zagier and Flip: 158 zagier_flip_1_1_z |- !z. (zagier o flip) (1,1,z) = 159 if z = 0 then (1,2,0) 160 else if z = 1 then (1,1,1) 161 else if z = 2 then (3,2,0) 162 else (3,1,z - 2) 163 164 Computation of (mills n): 165 tuples_helper_def |- (!k. tuples_helper k 0 = []) /\ 166 !k n. tuples_helper k (SUC n) = 167 ZIP (GENLIST (K (SUC n)) k, 168 GENLIST SUC k) ++ tuples_helper k n 169 tuples_def |- !k. tuples k = tuples_helper k k 170 triples_helper_def |- (!k. triples_helper k 0 = []) /\ 171 !k n. triples_helper k (SUC n) = 172 ZIP (GENLIST (K (SUC n)) (k * k), 173 tuples k) ++ triples_helper k n 174 triples_def |- !k. triples k = triples_helper k k 175 mills_of_def |- !n. mills_of n = 176 FILTER (\(x,y,z). n = windmill x y z) (triples n) 177*) 178 179(* ------------------------------------------------------------------------- *) 180(* Helper Theorems *) 181(* ------------------------------------------------------------------------- *) 182 183(* Theorem: ?x y z. t = (x, y, z) *) 184(* Proof: by x = FST t, FST (SND t), and SND (SND t). *) 185Theorem triple_parts: 186 !(t :num # num # num). ?x y z. t = (x, y, z) 187Proof 188 rpt strip_tac >> 189 qabbrev_tac `x = FST t` >> 190 qabbrev_tac `y = FST (SND t)` >> 191 qabbrev_tac `z = SND (SND t)` >> 192 `t = (x, y, z)` by rw[Abbr`x`, Abbr`y`, Abbr`z`] >> 193 metis_tac[] 194QED 195 196(* ------------------------------------------------------------------------- *) 197(* Windmill Theory *) 198(* ------------------------------------------------------------------------- *) 199 200(* 201 202 +---+ 203 | | 204 | | 205 | | 206 | | 207 | | y 208 | +---------------------+ 209 | | | z 210 | +-----+---+-----------+ 211 | | x | | 212 | | | | 213 +-----------+---+-----+ | 214 | | | 215 +---------------------+ | 216 | | 217 | | 218 | | 219 | | 220 | | 221 +---+ 222 223 Note: need n not a square, so that y * z <> 0. 224 Note: need n not a multiple of 4, so that x <> 0. 225 226 Note: need x * y * z <> 0 for the windmill picture: 227 228 Algorithm: 229 1. draw the square of side x. 230 2. from each square vertex, put the line y alongside, in clockwise manner. 231 3. complete the 4 rectangles y * z, around the central square. 232 233 234 The 5 types of windmills: 235 ------------------------- 236 237 Type 1: x < y - z, so x < y. 238 example: zagier (3,6,1) = (5,1,2) 239 240 * * * * 241 * * * * 242 * * * * * * * * * * * * * * * * 243 * * * * * * * * * * * * 244 * * * * * * 245 * * * * --> * * 246 * * * * * * * * * * * * 247 * * * * * * * * * * * * * * * * 248 * * * * 249 * * * * 250 251 square: 3x3 -> 5x5, mind = 5x5. 252 x' = x + left z + right z = x + 2z 253 y' = z 254 z' = y - (x + z) = y - x - z, because x + z < y. 255 256 Type 2: y - z < x, and x < y, so x < 2y. 257 example: zagier (3,6,4) = (9,6,1) 258 259 * * * * * * * * * * * * * * 260 * * * * * * * * * * * * * * * * * * * * * * 261 * * * * * * * * * * * * * * 262 * * * * * * * * * * * * * * 263 * * * * * * * * * * * * * * * * 264 * * * * * * * * * * * * * * 265 * * * * * * * * * * --> * * * * 266 * * * * * * * * * * * * * * * * 267 * * * * * * * * * * * * * * 268 * * * * * * * * * * * * * * 269 * * * * * * * * * * * * * * * * * * * * * * 270 * * * * * * * * * * * * * * 271 272 square: 3x3 -> 9x9, mind = 9x9. 273 x' = x + left (y - x) + right (y - x) = 2y - x 274 y' = y 275 z' = (x + z) - y = x + z - y, because y < x + z. 276 277 Type 3: y = x, so y - z < x, and x < 2y. 278 example: zagier (4,4,2) = (4,4,2) 279 280 * * * * * * * * * * 281 * * * * * * * * * * 282 * * * * * * * * * * * * * * * * * * 283 * * * * * * * * * * * * 284 * * * * * * --> * * * * * * 285 * * * * * * * * * * * * 286 * * * * * * * * * * * * * * * * * * 287 * * * * * * * * * * 288 * * * * * * * * * * 289 290 square: 4x4 -> 4x4, mind = 4x4. 291 x' = x 292 y' = y 293 z' = z, all unchanged. 294 295 Type 4: y < x, but x < 2y. 296 example: zagier (6,4,2) = (2,4,4) 297 298 * * * * * * * * * * 299 * * * * * * * * * * 300 * * * * * * * * * * * * * * * * * * 301 * * * * * * * * * * * * * 302 * * * * * * * * * * * * * * * * * 303 * * * * * * --> * * * * * * * * * * 304 * * * * * * * * * * * * * * * * * 305 * * * * * * * * * * * * * 306 * * * * * * * * * * * * * * * * * * 307 * * * * * * * * * * * 308 * * * * * * * * * * * 309 310 square: 6x6 -> 2x2, mind = 6x6. 311 x' = x - left (x - y) - right (x - y) = 2y - x 312 y' = y 313 z' = (x + z) - y = x + z - y, because y < x + z. 314 315 Type 5: 2y < x, so y < x. 316 example: zagier (8,3,2) = (2,7,3) 317 318 * * * * * * * * 319 * * * * * * * * 320 * * * * * * * * * * * * * * * * * * * * * * 321 * * * * * * * * * * * * * * * 322 * * * * * * * * * * * * * * * 323 * * * * * * * * * * * * * * * 324 * * --> * * * * * * * * 325 * * * * * * * * * * * * * * * 326 * * * * * * * * * * * * * * * 327 * * * * * * * * * * * * * * * 328 * * * * * * * * * * * * * * * * * * * * * * 329 * * * * * * * * 330 * * * * * * * * 331 332 square: 8x8 -> 2x2, mind = 6x6. 333 x' = x - left y - right y = x - 2y 334 y' = (x + z) - y = x + z - y, because y < x + z. 335 z' = y. 336 337 "mind" is the maximum central square: 338 if x < y - z, central square expands to x + left z + right z = x + 2z. 339 else if x < y, central square expands to x + left (y - x) + right (y - x) = 2y - x. 340 else central square is at maximum, keeps as x. 341 342*) 343 344(* ------------------------------------------------------------------------- *) 345(* Part 1: A windmill *) 346(* ------------------------------------------------------------------------- *) 347 348(* Define windmill of three numbers *) 349Definition windmill_def: 350 windmill x y z = x ** 2 + 4 * y * z 351End 352 353(* Theorem: windmill x y z = 0 <=> x = 0 /\ y * z = 0 *) 354(* Proof: 355 windmill x y z = 0 356 <=> x ** 2 + 4 * y * z = 0 by windmill_def 357 <=> x * x + 4 * y * z = 0 by EXP_2 358 <=> (x * x = 0) /\ (4 * y * z = 0) by ADD_EQ_0 359 <=> (x = 0) /\ (y * z = 0) by MULT_EQ_0 360*) 361Theorem windmill_eq_0: 362 !x y z. windmill x y z = 0 <=> x = 0 /\ y * z = 0 363Proof 364 simp[windmill_def] 365QED 366 367(* Theorem: windmill x y z = windmill x z y *) 368(* Proof: 369 windmill x y z 370 = x ** 2 + 4 * y * z by windmill_def 371 = x ** 2 + 4 * z * y by MULT_COMM 372 = windmill x z y by windmill_def 373*) 374Theorem windmill_comm: 375 !x y z. windmill x y z = windmill x z y 376Proof 377 simp[windmill_def] 378QED 379 380(* Theorem: windmill 1 1 k = 4 * k + 1 *) 381(* Proof: 382 windmill 1 1 k 383 = 1 ** 2 + 4 * 1 * k by windmill_def 384 = 1 + 4 * k by arithmetic 385 = 4 * k + 1 by arithmetic 386*) 387Theorem windmill_trivial: 388 !k. windmill 1 1 k = 4 * k + 1 389Proof 390 simp[windmill_def] 391QED 392 393(* Theorem: windmill x y y = x ** 2 + (2 * y) ** 2 *) 394(* Proof: 395 windmill x y y 396 = x ** 2 + 4 * y * y by windmill_def 397 = x ** 2 + (2 * y) * (2 * y) by arithmetic 398 = x ** 2 + (2 * y) ** 2 by EXP_2 399*) 400Theorem windmill_by_squares: 401 !x y. windmill x y y = x ** 2 + (2 * y) ** 2 402Proof 403 simp[windmill_def, EXP_BASE_MULT] 404QED 405 406(* Theorem: n = windmill x y y ==> 407 n = x ** 2 + (2 * y) ** 2 /\ (ODD n <=> ODD x) *) 408(* Proof: 409 Note n = x ** 2 + (2 * y) ** 2 by windmill_by_squares 410 and EVEN (2 * y) by EVEN_DOUBLE 411 so EVEN (2 * y) ** 2 by EVEN_SQ 412 ODD n 413 <=> ODD (x ** 2) by ODD_ADD, ODD_EVEN 414 <=> ODD x by ODD_SQ 415*) 416Theorem windmill_x_y_y: 417 !n x y. n = windmill x y y ==> 418 n = x ** 2 + (2 * y) ** 2 /\ (ODD n <=> ODD x) 419Proof 420 ntac 4 strip_tac >> 421 `n = x ** 2 + (2 * y) ** 2` by rw[windmill_by_squares] >> 422 `EVEN (2 * y)` by rw[EVEN_DOUBLE] >> 423 `EVEN ((2 * y) ** 2)` by rw[EVEN_SQ] >> 424 `ODD n <=> ODD (x ** 2)` by metis_tac[ODD_ADD, ODD_EVEN] >> 425 simp[ODD_SQ] 426QED 427 428(* Theorem: prime p /\ p MOD 4 = 1 ==> 429 !x z. p = windmill x x z <=> x = 1 /\ z = p DIV 4 *) 430(* Proof: 431 Let k = p DIV 4, 432 Then p = 4 * k + 1 by DIVISION 433 434 If part: p = windmill x x z ==> ((x = 1) /\ (z = k)) 435 Note p = windmill x x z 436 = x ** 2 + 4 * x * z by windmill_def 437 = x * (x + 4 * z) by arithmetic 438 Thus x divides p by divides_def 439 so x = 1 or x = p by prime_def 440 If x = 1, 441 p = 4 * k + 1 442 = 1 + 4 * z by above 443 so z = k by arithmetic 444 If x = p, 445 then p + 4 * z = 1 by EQ_MULT_LCANCEL 446 but 1 < p by ONE_LT_PRIME 447 and this is impossible. 448 Only-if part: p = windmill 1 1 k 449 Note p = 4 * k + 1 450 = windmill 1 1 k by windmill_trivial 451*) 452Theorem windmill_trivial_prime: 453 !p. prime p /\ p MOD 4 = 1 ==> 454 !x z. p = windmill x x z <=> x = 1 /\ z = p DIV 4 455Proof 456 rpt strip_tac >> 457 qabbrev_tac `k = p DIV 4` >> 458 `p = k * 4 + 1` by metis_tac[DIVISION, DECIDE``0 < 4``] >> 459 `_ = 4 * k + 1` by simp[] >> 460 simp[EQ_IMP_THM] >> 461 ntac 2 strip_tac >| [ 462 `p = x ** 2 + 4 * x * z` by rw[windmill_def] >> 463 `_ = x * x + 4 * z * x` by simp[] >> 464 `_ = (x + 4 * z) * x` by decide_tac >> 465 `x divides p` by metis_tac[divides_def] >> 466 `(x = 1) \/ (x = p)` by metis_tac[prime_def] >- 467 fs[] >> 468 fs[] >> 469 `p + 4 * z = 1` by metis_tac[MULT_RIGHT_1, NOT_PRIME_0, EQ_MULT_LCANCEL] >> 470 `1 < p` by rw[ONE_LT_PRIME] >> 471 decide_tac, 472 rw[windmill_trivial] 473 ] 474QED 475 476(* ------------------------------------------------------------------------- *) 477(* Part 2: Set of windmills *) 478(* ------------------------------------------------------------------------- *) 479 480 481(* The set of windmills of a number *) 482val mills_def = zDefine` 483 mills n = {(x,y,z) | n = windmill x y z} 484`; 485(* use zDefine as this is not effective. *) 486 487(* Theorem: (x, y, z) IN mills n <=> n = windmill x y z *) 488(* Proof: by mills_def. *) 489Theorem mills_element: 490 !n x y z. (x, y, z) IN mills n <=> n = windmill x y z 491Proof 492 simp[mills_def] 493QED 494 495(* Theorem: (x, y, z) IN mills n ==> (x, z, y) IN mills n *) 496(* Proof: 497 (x, y, z) IN mills n 498 <=> n = windmill x y z by mills_def 499 <=> n = x ** 2 + 4 * y * z by windmill_def 500 <=> n = x ** 2 + 4 * z * y by MULT_COMM 501 <=> n = windmill x z y by windmill_def 502 <=> (x, z, y) IN mills n by mills_def 503*) 504Theorem mills_element_flip: 505 !n x y z. (x, y, z) IN mills n ==> (x, z, y) IN mills n 506Proof 507 simp[mills_def, windmill_def] 508QED 509 510(* Theorem: t IN mills n <=> ?x y z. (t = (x, y, z)) /\ n = windmill x y z *) 511(* Proof: by mills_def. *) 512Theorem mills_element_triple: 513 !n t. t IN mills n <=> ?x y z. (t = (x, y, z)) /\ n = windmill x y z 514Proof 515 simp[mills_def, FORALL_PROD] 516QED 517 518(* Theorem: n MOD 4 = 1 ==> (1, 1, n DIV 4) IN mills n *) 519(* Proof: 520 Note n = (n DIV 4) * 4 + n MOD 4 by DIVISION 521 = 4 * (n DIV 4) + 1 by arithmetic 522 = windmill 1 1 (n DIV 4) by windmill_trivial 523 Thus (1, 1, n DIV 4) IN (mills n) by mills_element 524*) 525Theorem mills_element_trivial: 526 !n. n MOD 4 = 1 ==> (1, 1, n DIV 4) IN mills n 527Proof 528 rpt strip_tac >> 529 `n = (n DIV 4) * 4 + n MOD 4` by rw[DIVISION] >> 530 `_ = windmill 1 1 (n DIV 4)` by rw[windmill_trivial] >> 531 fs[mills_element] 532QED 533 534(* Theorem: mills 0 = {(0,0,n) | n IN univ(:num)} UNION 535 {(0,n,0) | n IN univ(:num)} *) 536(* Proof: 537 Let (x,y,z) IN mills 0 538 Then 0 = windmill x y z by mills_def 539 = x ** 2 + 4 * y * z by windmill_def 540 ==> x ** 2 = 0 /\ 4 * y * z = 0 by ADD_EQ_0 541 ==> x * x = 0 /\ y * z = 0 by EXP_2 542 ==> x = 0 /\ (y = 0 \/ z = 0) by MULT_EQ_0 543*) 544Theorem mills_0: 545 mills 0 = {(0,0,n) | n IN univ(:num)} UNION 546 {(0,n,0) | n IN univ(:num)} 547Proof 548 rw[mills_def, windmill_def, FORALL_PROD, EXTENSION] 549QED 550 551(* Theorem: INFINITE (mills 0) *) 552(* Proof: 553 Let f = (\n. (0,0,n)). 554 Then INJ f univ(:num) (mills 0) by INJ_DEF, mills_0 555 Note INFINITE univ(:num) by INFINITE_NUM_UNIV 556 ==> INFINITE (mills 0) by INFINITE_INJ 557*) 558Theorem mills_0_infinite: 559 INFINITE (mills 0) 560Proof 561 qabbrev_tac `f = \n:num. (0,0,n)` >> 562 `INJ f univ(:num) (mills 0)` by rw[INJ_DEF, mills_0, Abbr`f`] >> 563 `INFINITE univ(:num)` by rw[] >> 564 metis_tac[INFINITE_INJ] 565QED 566 567(* Theorem: mills 1 = {1} CROSS ({0} CROSS univ(:num)) UNION 568 {1} CROSS (univ(:num) CROSS {0}) *) 569(* Proof: 570 Let (x,y,z) IN mills 1 571 Note 4 * y * z <> 1 by MULT_EQ_1 572 Then 1 = windmill x y z by mills_def 573 = x ** 2 + 4 * y * z by windmill_def 574 <=> x ** 2 = 1 /\ 4 * y * z = 0 by ADD_EQ_1 575 <=> x * x = 1 /\ y * z = 0 by EXP_2 576 <=> x = 1 /\ (y = 0 \/ z = 0) by MULT_EQ_0 577 <=> (x,y,z) IN (1,0,n) where n IN univ(:num) 578 or (x,y,z) IN (1,n,0) where n IN univ(:num) 579 <=> (x,y,z) IN {1} CROSS ({0} CROSS univ(:num)) 580 or (x,y,z) IN {1} CROSS (univ(:num) CROSS {0}) 581*) 582Theorem mills_1: 583 mills 1 = {1} CROSS ({0} CROSS univ(:num)) UNION 584 {1} CROSS (univ(:num) CROSS {0}) 585Proof 586 rw[mills_def, windmill_def, FORALL_PROD, EXTENSION, ADD_EQ_1] 587QED 588 589(* Theorem: INFINITE (mills 1) *) 590(* Proof: 591 Let f = (\n. (1,0,n)). 592 Then INJ f univ(:num) (mills 0) by INJ_DEF, mills_1 593 Note INFINITE univ(:num) by INFINITE_NUM_UNIV 594 ==> INFINITE (mills 0) by INFINITE_INJ 595*) 596Theorem mills_1_infinite: 597 INFINITE (mills 1) 598Proof 599 qabbrev_tac `f = \n:num. (1,0,n)` >> 600 `INJ f univ(:num) (mills 1)` by rw[INJ_DEF, mills_1, Abbr`f`] >> 601 `INFINITE univ(:num)` by rw[] >> 602 metis_tac[INFINITE_INJ] 603QED 604 605(* Theorem: ~square n /\ (x, y, z) IN mills n ==> x < n /\ y < n /\ z < n *) 606(* Proof: 607 Expand by square_def, mills_def, windmill_def, this is to show, 608 given !k. 4 * (y * z) + x ** 2 <> k ** 2: 609 (1) x < 4 * (y * z) + x ** 2 610 Note y * z <> 0 by given condition 611 so y <> 0 /\ z <> 0 by MULT_EQ_0 612 or 4 * (y * z) <> 0 by MULT_EQ_0 613 Note x <= x ** 2 by X_LE_X_SQUARED 614 Hence x < 4 * (y * z) + x ** 2 by adding nonzero term to RHS. 615 (2) y < 4 * (y * z) + x ** 2 616 Note y * z <> 0 by given condition 617 so y <> 0 /\ z <> 0 by MULT_EQ_0 618 thus 4 * z <> 0 by MULT_EQ_0, z <> 0 619 and 4 * z <> 1 by MULT_EQ_1, 4 <> 1 620 so y < y * (4 * z) by LT_MULT_CANCEL_LBARE 621 Hence y < 4 * (y * z) + x ** 2 by adding to RHS. 622 (3) y <> 0 /\ z <> 0 ==> z < 4 * (y * z) + x ** 2 623 Note y * z <> 0 by given condition 624 so y <> 0 /\ z <> 0 by MULT_EQ_0 625 thus 4 * y <> 0 by MULT_EQ_0, y <> 0 626 and 4 * y <> 1 by MULT_EQ_1, 4 <> 1 627 so z < z * (4 * y) by LT_MULT_CANCEL_LBARE 628 Hence z < 4 * (y * z) + x ** 2 by adding to RHS. 629*) 630Theorem mills_non_square_bound: 631 !n x y z. ~square n /\ (x, y, z) IN mills n ==> x < n /\ y < n /\ z < n 632Proof 633 rw[square_def, mills_def, windmill_def] >| [ 634 `y * z <> 0` by metis_tac[MULT_0, ADD] >> 635 `4 * (y * z) <> 0` by metis_tac[MULT_EQ_0, DECIDE``4 <> 0``] >> 636 `x <= x ** 2` by rw[X_LE_X_SQUARED] >> 637 decide_tac, 638 `y * z <> 0` by metis_tac[MULT_0, ADD] >> 639 `y <> 0 /\ z <> 0` by metis_tac[MULT_EQ_0] >> 640 `4 * z <> 0` by metis_tac[MULT_EQ_0, DECIDE``4 <> 0``] >> 641 `4 * z <> 1` by metis_tac[MULT_EQ_1, DECIDE``4 <> 1``] >> 642 `y < y * (4 * z)` by rw[LT_MULT_CANCEL_LBARE] >> 643 decide_tac, 644 `y * z <> 0` by metis_tac[MULT_0, ADD] >> 645 `y <> 0 /\ z <> 0` by metis_tac[MULT_EQ_0] >> 646 `4 * y <> 0` by metis_tac[MULT_EQ_0, DECIDE``4 <> 0``] >> 647 `4 * y <> 1` by metis_tac[MULT_EQ_1, DECIDE``4 <> 1``] >> 648 `z < z * (4 * y)` by rw[LT_MULT_CANCEL_LBARE] >> 649 decide_tac 650 ] 651QED 652 653(* Theorem: ~square n ==> 654 mills n SUBSET (count n) CROSS (count n CROSS (count n)) *) 655(* Proof: by SUBSET_DEF, mills_non_square_bound, IN_COUNT *) 656Theorem mills_non_square_subset: 657 !n. ~square n ==> 658 mills n SUBSET (count n) CROSS (count n CROSS (count n)) 659Proof 660 (rw[SUBSET_DEF, FORALL_PROD] >> metis_tac[mills_non_square_bound]) 661QED 662 663(* Theorem: ~square n ==> CARD (mills n) < n ** 3 *) 664(* Proof: 665 Let s = count n CROSS (count n CROSS count n), 666 t = mills n. 667 Then t SUBSET s by mills_non_square_subset 668 and FINITE s by FINITE_CROSS, FINITE_COUNT 669 Note n <> 0 by square_0, ~square n. 670 and n <> 1 by square_1, ~square n. 671 so 1 < n by arithmetic 672 Thus (1,1,0) IN s by IN_CROSS, IN_COUNT, 1 < n 673 but windmill 1 1 0 674 = 1 <> n by windmill_trivial 675 so (1,1,0) NOTIN t by mills_element 676 ==> t <> s by EXTENSION 677 ==> t PSUBSET s by PSUBSET_DEF 678 Now CARD s = n ** 3 by CARD_CROSS, CARD_COUNT 679 so CARD t < n ** 3 by CARD_PSUBSET 680*) 681Theorem mills_non_square_card_upper: 682 !n. ~square n ==> CARD (mills n) < n ** 3 683Proof 684 rpt strip_tac >> 685 qabbrev_tac `s = count n CROSS (count n CROSS count n)` >> 686 qabbrev_tac `t = mills n` >> 687 `t SUBSET s` by fs[mills_non_square_subset, Abbr`t`, Abbr`s`] >> 688 `FINITE s` by rw[Abbr`s`] >> 689 (Cases_on `n = 0` >> fs[square_def]) >> 690 (Cases_on `n = 1` >> fs[square_def]) >> 691 `1 < n` by decide_tac >> 692 `(1,1,0) IN s` by rw[Abbr`s`] >> 693 `(1,1,0) NOTIN t` by rw[mills_element, windmill_trivial, Abbr`t`] >> 694 `t <> s` by metis_tac[EXTENSION] >> 695 `t PSUBSET s` by fs[PSUBSET_DEF] >> 696 `CARD s = n ** 3` by rw[CARD_CROSS, Abbr`s`] >> 697 metis_tac[CARD_PSUBSET] 698QED 699 700(* Theorem: ~square n ==> FINITE (mills n) *) 701(* Proof: 702 Let c = count n. 703 Note FINITE c by FINITE_COUNT 704 so FINITE (c CROSS c CROSS c) by FINITE_CROSS 705 Now mills n SUBSET (c CROSS c CROSS c) by mills_non_square_subset 706 so FINITE (mills n) by SUBSET_FINITE 707*) 708Theorem mills_non_square_finite: 709 !n. ~square n ==> FINITE (mills n) 710Proof 711 metis_tac[mills_non_square_subset, SUBSET_FINITE, FINITE_CROSS, FINITE_COUNT] 712QED 713 714(* Theorem: FINITE (mills n) <=> ~square n *) 715(* Proof: 716 If part: FINITE (mills n) ==> ~square n 717 By contradiction, suppose square n. 718 Then ?k. n = k * k = k ** 2 by square_def 719 Then n = k ** 2 + 4 * 0 * t for any t, 720 = windmill k 0 t by windmill_def 721 so (k,0,t) IN (mills n) 722 Let f = \t:num. (k,0,t)`); 723 Then INJ f univ(:num) (mills n) by INJ_DEF 724 But INFINITE univ(:num) by INFINITE_NUM_UNIV 725 ==> INFINITE (mills n) by INFINITE_INJ 726 This contradicts FINITE (mills n). 727 Only-if part: ~square n ==> FINITE (mills n) 728 This is true by mills_non_square_finite 729*) 730Theorem mills_finite_non_square: 731 !n. FINITE (mills n) <=> ~square n 732Proof 733 rw[EQ_IMP_THM] >| [ 734 spose_not_then strip_assume_tac >> 735 fs[square_def] >> 736 qabbrev_tac `f = \t:num. (k, 0, t)` >> 737 `INJ f univ(:num) (mills n)` by 738 (rw[INJ_DEF, mills_def, windmill_def] >| [ 739 qexists_tac `k` >> 740 qexists_tac `0` >> 741 qexists_tac `x` >> 742 simp[], 743 fs[PAIR_EQ, Abbr`f`] 744 ]) >> 745 metis_tac[INFINITE_NUM_UNIV, INFINITE_INJ], 746 rw[mills_non_square_finite] 747 ] 748QED 749 750(* Theorem: (?y z. (0, y, z) IN mills n) <=> n MOD 4 = 0 *) 751(* Proof: 752 If part: (?y z. (0, y, z) IN mills n) ==> n MOD 4 = 0 753 (0, y, z) IN mills n 754 ==> n = windmill 0 y z by mills_def 755 ==> n = 0 ** 2 + 4 * y * z by windmill_def 756 ==> n = 4 * y * z by arithmetic 757 ==> n MOD 4 = 0 by MOD_EQ_0 758 Only-if part: n MOD 4 = 0 ==> (?y z. (0, y, z) IN mills n) 759 Note n = n DIV 4 * 4 + n MOD 4 by DIVISION 760 Let y = n DIV 4, z = 1. 761 Then n = 0 ** 2 + 4 * n DIV 4 by n MOD 4 = 0 762 ==> (0, y, z) IN mills n 763*) 764Theorem mills_with_no_mind: 765 !n. (?y z. (0, y, z) IN mills n) <=> n MOD 4 = 0 766Proof 767 rw[mills_def, windmill_def] >> 768 (rw[EQ_IMP_THM] >> simp[]) >> 769 `n = n DIV 4 * 4 + n MOD 4` by rw[DIVISION] >> 770 qexists_tac `n DIV 4` >> 771 qexists_tac `1` >> 772 simp[] 773QED 774 775(* Theorem: n MOD 4 <> 0 <=> !x y z. (x, y, z) IN mills n ==> x <> 0 *) 776(* Proof: by mills_with_no_mind *) 777Theorem mills_with_all_mind: 778 !n. n MOD 4 <> 0 <=> !x y z. (x, y, z) IN mills n ==> x <> 0 779Proof 780 metis_tac[mills_with_no_mind] 781QED 782 783(* Theorem: (?x z. (x, 0, z) IN mills n) \/ (?x y. (x, y, 0) IN mills n) 784 <=> square n *) 785(* Proof: 786 If part: (?x z. (x, 0, z) IN mills n) \/ (?x y. (x, y, 0) IN mills n) 787 ==> square n 788 (x, 0, z) IN mills n 789 ==> n = windmill x 0 z by mills_def 790 ==> n = x ** 2 + 4 * 0 * z by windmill_def 791 ==> n = x ** 2 by arithmetic 792 ==> square n by square_def, EXP_2 793 (x, y, 0) IN mills n 794 ==> n = windmill x y 0 by mills_def 795 ==> n = x ** 2 + 4 * y * 0 by windmill_def 796 ==> n = x ** 2 by arithmetic 797 ==> square n by square_def, EXP_2 798 Only-if part: square n ==> (?x z. (x, 0, z) IN mills n) \/ (?x y. (x, y, 0) IN mills n) 799 Note ?k. n = k * k by square_def 800 Let x = k, any z in first case, any y in second case, 801 then n = x ** 2 + 4 * 0 * z, so (x, 0, z) IN mills n. 802 or n = x ** 2 + 4 * y * 0, so (x, y, 0) IN mills n. 803*) 804Theorem mills_with_no_arms: 805 !n. (?x z. (x, 0, z) IN mills n) \/ (?x y. (x, y, 0) IN mills n) <=> square n 806Proof 807 rw[mills_def, windmill_def, square_def] 808QED 809 810(* Theorem: ~square n <=> !x y z. (x, y, z) IN mills n ==> y <> 0 /\ z <> 0 *) 811(* Proof: by mills_with_no_arms *) 812Theorem mills_with_all_arms: 813 !n. ~square n <=> !x y z. (x, y, z) IN mills n ==> y <> 0 /\ z <> 0 814Proof 815 metis_tac[mills_with_no_arms] 816QED 817 818(* Theorem: n MOD 4 = 1 ==> mills n <> {} *) 819(* Proof: 820 By contradiction, suppose (mills n = EMPTY). 821 Now ?k. n = k * 4 + 1 by DIVISION 822 so n = 1 ** 2 + 4 * 1 * k by arithmetic 823 = windmill 1 1 k by windmill_def 824 Thus (1, 1, k) IN mills n by mills_def 825 This contradicts mills n = EMPTY by MEMBER_NOT_EMPTY 826*) 827Theorem mills_quad_suc_non_empty: 828 !n. n MOD 4 = 1 ==> mills n <> {} 829Proof 830 rpt strip_tac >> 831 `?k. n = k * 4 + 1` by metis_tac[DIVISION, DECIDE``0 < 4``] >> 832 `_ = windmill 1 1 k` by rw[windmill_def] >> 833 `(1, 1, k) IN (mills n)` by metis_tac[mills_element] >> 834 metis_tac[MEMBER_NOT_EMPTY] 835QED 836 837(* Theorem: ~square n /\ n MOD 4 <> 0 ==> 838 !x y z. (x,y,z) IN (mills n) ==> x <> 0 /\ y <> 0 /\ z <> 0 *) 839(* Proof: 840 Note n = windmill x y z by mills_def 841 = x ** 2 + 4 * y * z by windmill_def 842 By contradiction, suppose x = 0, or y = 0, or z = 0. 843 If x = 0, then n MOD 4 = 0 by MOD_EQ_0 844 which contradicts n MOD 4 <> 0. 845 If y = 0 or z = 0, then n = x ** 2 846 which contradicts ~square n 847 by square_def, EXP_2 848*) 849Theorem mills_triple_nonzero: 850 !n. ~square n /\ n MOD 4 <> 0 ==> 851 !x y z. (x,y,z) IN (mills n) ==> x <> 0 /\ y <> 0 /\ z <> 0 852Proof 853 spose_not_then strip_assume_tac >> 854 rfs[square_def, mills_def, windmill_def] >> 855 `y <> 0 /\ z <> 0` by metis_tac[MULT_EQ_0, ADD] >> 856 `x <> 0` by metis_tac[SQ_0, ADD_0, MULT_COMM, MOD_EQ_0, DECIDE``0 < 4``] >> 857 decide_tac 858QED 859 860(* Theorem: prime p /\ (x,y,z) IN (mills p) ==> x <> 0 /\ y <> 0 /\ z <> 0 *) 861(* Proof: 862 Note ~square p by prime_non_square 863 and p MOD 4 <> 0 by prime_mod_eq_0, NOT_PRIME_4 864 so (x,y,z) IN (mills p) 865 ==> x <> 0 /\ y <> 0 /\ z <> 0 by mills_triple_nonzero 866*) 867Theorem mills_prime_triple_nonzero: 868 !p x y z. prime p /\ (x,y,z) IN (mills p) ==> x <> 0 /\ y <> 0 /\ z <> 0 869Proof 870 ntac 5 strip_tac >> 871 `~square p` by metis_tac[prime_non_square] >> 872 `p MOD 4 <> 0` by metis_tac[prime_mod_eq_0, NOT_PRIME_4, DECIDE``1 < 4``] >> 873 metis_tac[mills_triple_nonzero] 874QED 875 876(* ------------------------------------------------------------------------- *) 877(* Flip involution. *) 878(* ------------------------------------------------------------------------- *) 879 880(* Use temporary overload, leave proper overload in involute.hol *) 881val _ = temp_overload_on("involute", ``\f s. !x. x IN s ==> f x IN s /\ (f (f x) = x)``); 882val _ = set_fixity "involute" (Infix(NONASSOC, 450)); (* same as relation *) 883 884(* Define the flip function for last two elements of a triple. *) 885Definition flip_def: 886 flip (x:num, y:num, z:num) = (x, z, y) 887End 888 889(* Theorem: flip (x,y,z) = (x,y,z) <=> y = z *) 890(* Proof: 891 flip (x,y,z) = (x,y,z) 892 <=> (x,z,y) = (x,y,z) by flip_def 893 <=> x = x, z = y, y = z. 894*) 895Theorem flip_fix: 896 !x y z. flip (x,y,z) = (x,y,z) <=> y = z 897Proof 898 simp[flip_def] 899QED 900 901(* Theorem: (x,y,z) IN mills n ==> flip (x,y,z) IN mills n *) 902(* Proof: by flip_def, mills_element_flip. *) 903Theorem flip_closure: 904 !n x y z. (x,y,z) IN mills n ==> flip (x,y,z) IN mills n 905Proof 906 rw[flip_def, mills_element_flip] 907QED 908 909(* Theorem: t IN mills n ==> flip t IN mills n *) 910(* Proof: by flip_closure. *) 911Theorem flip_closure_alt: 912 !n t. t IN mills n ==> flip t IN mills n 913Proof 914 metis_tac[triple_parts, flip_closure] 915QED 916 917(* Theorem: flip (flip (x,y,z)) = (x,y,z) *) 918(* Proof: by flip_def *) 919Theorem flip_involute: 920 !x y z. flip (flip (x,y,z)) = (x,y,z) 921Proof 922 rw[flip_def] 923QED 924 925(* Theorem: flip (flip t) = t *) 926(* Proof: by flip_involute. *) 927Theorem flip_involute_alt: 928 !t. flip (flip t) = t 929Proof 930 metis_tac[triple_parts, flip_involute] 931QED 932 933(* Theorem: flip involute (mills n) *) 934(* Proof: 935 Let t = (x,y,z). by triple_parts 936 Then flip (x,y,z) IN mills n by flip_closure 937 and flip (flip (x,y,z)) = (x,y,z) by flip_involute 938 so flip involute (mills n) by involution 939*) 940Theorem flip_involute_mills: 941 !n. flip involute (mills n) 942Proof 943 metis_tac[flip_closure, flip_involute, triple_parts] 944QED 945 946(* flip_involute_mills |> SPEC ``n:num``; 947val it = |- flip involute mills n: thm 948*) 949 950(* ------------------------------------------------------------------------- *) 951(* Zagier's involution. *) 952(* *) 953(* A One-Sentence Proof That *) 954(* Every Prime p = 1 (mod 4) Is a Sum of Two Squares (Don Zagier) *) 955(* http://igor-kortchemski.perso.math.cnrs.fr/mathclub/Zagier.pdf *) 956(* The American Mathematical Monthly, Vol. 97, No. 2 (Feb., 1990), p. 144 *) 957(* ------------------------------------------------------------------------- *) 958 959(* Define the Zagier's involution: a self-inverse bijection *) 960Definition zagier_def: 961 zagier (x, y, z) = 962 if x < y - z then (x + 2 * z, z, y - z - x) 963 else if x < 2 * y then (2 * y - x, y, x + z - y) (* 2 * y - x = y + y - x *) 964 else (x - 2 * y, x + z - y, y) (* here y - z <= x /\ y <= x *) 965End 966(* 967At the two boundaries: 968x = y - z, windmill x y z = (y - z) ** 2 + 4 * y * z = (y + z) ** 2, not a windmill. 969x = 2 * y, windmill x y z = (2 * y) ** 2 + 4 * y * z = 4 * y * (y + z), not for a prime. 970*) 971 972 973(* 974For p = 41 = 4 * 10 + 1, k = 10. 975 976EVAL ``zagier (3,8,1)``; -> (5,1,4) -> (3,8,1) 977EVAL ``zagier (3,4,2)``; -> (5,4,1) -> (3,4,2) 978EVAL ``zagier (1,10,1)``; -> (3,1,8) -> (1,10,1) 979EVAL ``zagier (1,5,2)``; -> (5,2,2) -> (1,5,2) 980EVAL ``zagier (1,2,5)``; -> (3,2,4) -> (1,2,5) 981EVAL ``zagier (1,1,10)``; -> (1,1,10) -> (1,1,10) 982 983EVAL ``MAP zagier [(3,8,1);(3,4,2);(1,10,1);(1,5,2);(1,2,5);(1,1,10)]``; 984-> [(5,1,4); (5,4,1); (3,1,8); (5,2,2); (3,2,4); (1,1,10)]: thm 985EVAL ``MAP zagier [(5,1,4);(5,4,1);(3,1,8);(5,2,2);(3,2,4);(1,1,10)]``; 986-> [(3,8,1); (3,4,2); (1,10,1); (1,5,2); (1,2,5); (1,1,10)]: thm 987*) 988 989 990(* Theorem: zagier (0, y, z) = 991 if z < y then (2 * z,z,y - z) 992 else if 0 < y then (2 * y,y,z - y) else (0,z,0) *) 993(* Proof: 994 If x = 0, for 0 < y - z, need z < y. 995 In this case, 0 < y, so 0 < 2 * y. 996 Otherwise, y <= z. If ~(0 < y), then y = 0 997 zagier (0, y, z) 998 = if z < y then (2 * z,z,y - z) 999 else if 0 < y then (2 * y,y,z - y) 1000 else (0,z,0) 1001*) 1002Theorem zagier_0_y_z: 1003 !y z. zagier (0, y, z) = 1004 if z < y then (2 * z,z,y - z) 1005 else if 0 < y then (2 * y,y,z - y) 1006 else (0,z,0) 1007Proof 1008 rw[zagier_def] 1009QED 1010 1011(* Theorem: zagier (x, 0, z) = (x, x + z, 0) *) 1012(* Proof: 1013 If y = 0, y - z = 0, so x < y - z <=> F. 1014 Also x < 2 * y <=> F. Thus this is the third case: 1015 zagier (x, 0, z) 1016 = (x - 2 * 0,x + z - 0,0) 1017 = (x, x + z, 0) 1018*) 1019Theorem zagier_x_0_z: 1020 !x z. zagier (x, 0, z) = (x, x + z, 0) 1021Proof 1022 simp[zagier_def] 1023QED 1024 1025(* Theorem: zagier (x, y, 0) = 1026 if x < y then (x, 0, y - x) 1027 else if x < 2 * y then (2 * y - x, y, x - y) 1028 else (x - 2 * y, x - y, y) *) 1029(* Proof: 1030 If z = 0, y - z = y, so x < y - z <=> x < y. 1031 zagier (x, y, 0) 1032 = if x < y then (x + 2 * 0,0,y - 0 - x) 1033 if x < 2 * y then (2 * y - x,y,x + 0 - y) 1034 else (x - 2 * y,x + 0 - y,y) 1035 = if x < y then (x, 0, y - x) 1036 else if x < 2 * y then (2 * y - x, y, x - y) 1037 else (x - 2 * y, x - y, y) 1038*) 1039Theorem zagier_x_y_0: 1040 !x y. zagier (x, y, 0) = 1041 if x < y then (x, 0, y - x) 1042 else if x < 2 * y then (2 * y - x, y, x - y) 1043 else (x - 2 * y, x - y, y) 1044Proof 1045 simp[zagier_def] 1046QED 1047 1048(* Theorem: zagier (1, 1, z) = (1, 1, z) *) 1049(* Proof: 1050 If x = 1 and y = 1, then x < y - z <=> F, but x < 2 * y <=> T. 1051 zagier (1, 1, z) 1052 = (2 * 1 - 1,1,1 + z - 1) 1053 = (1,1,z) 1054*) 1055Theorem zagier_1_1_z: 1056 !z. zagier (1, 1, z) = (1, 1, z) 1057Proof 1058 simp[zagier_def] 1059QED 1060 1061(* Theorem: zagier (1, y, 1) = 1062 if y = 0 then (1, 2, 0) 1063 else if y = 1 then (1, 1, 1) 1064 else if y = 2 then (3, 2, 0) 1065 else (3, 1, y - 2) *) 1066(* Proof: 1067 If x = 1 and z = 1, then x < y - z <=> 2 < y, 1068 and zagier (1, y, 1) = (1 + 2 * 1,1,y - 1 - 1) = (3, 1, y - 2) 1069 Otherwise, y <= 2, so y = 0, 1, 2. 1070 For x < 2 * y <=> 1 < 2 * y, y = 1 or 2 1071 and zagier (1, y, 1) = (2 * y - 1,y,1 + 1 - y) = (2 * y - 1, y, 2 - y) 1072 when y = 1, zagier (1, 1, 1) = (1, 1, 1) 1073 when y = 2, zagier (1, 2, 1) = (3, 2, 0) 1074 Otherwise y = 0, zagier (1, 0, 1) = (1 - 2 * 0,1 + 1 - 0,0) = (1, 2, 0) 1075*) 1076Theorem zagier_1_y_1: 1077 !y. zagier (1, y, 1) = 1078 if y = 0 then (1, 2, 0) 1079 else if y = 1 then (1, 1, 1) 1080 else if y = 2 then (3, 2, 0) 1081 else (3, 1, y - 2) 1082Proof 1083 rw[zagier_def] 1084QED 1085 1086(* 1087g `!n x y z. (x, y, z) IN mills n ==> my_zagier (x, y, z) IN mills n`; 1088e (rw[mills_def, windmill_def, my_zagier_def]); (* >> 4 *) 1089 1090g `!n x y z. (x, y, z) IN mills n ==> me_zagier (x, y, z) IN mills n`; 1091e (rw[mills_def, windmill_def, me_zagier_def]); (* >> 3 *) 1092but me_zagier is wrong! 1093 1094g `!n x y z. (x, y, z) IN mills n ==> mi_zagier (x, y, z) IN mills n`; 1095e (rw[mills_def, windmill_def, mi_zagier_def]); (* >> 4 *) 1096*) 1097 1098(* Theorem: zagier (x,0,0) = (x,x,0) *) 1099(* Proof: 1100 If y = 0 and z = 0, x < y - z becomes x < 0, which is false. 1101 Next x < 2 * y becomes x < 2 * 0 = 0 is false, too. 1102 Thus zagier (x,0,0) = (x - 2 * y,x + z - y,y) = (x,x,0). 1103*) 1104Theorem zagier_x_0_0: 1105 !x. zagier (x,0,0) = (x,x,0) 1106Proof 1107 simp[zagier_def] 1108QED 1109 1110(* Theorem: zagier (0,y,0) = (0,0,y) *) 1111(* Proof: 1112 If x = 0 and z = 0, x < y - z becomes 0 < y. 1113 If 0 < y, 1114 then zagier (0,y,0) = (x + 2 * z,z,y - z - x) = (0,0,y) 1115 Otherwise y = 0, and x < 2 * y becomes 0 < 2 * 0 = 0 is false. 1116 Thus zagier (0,0,0) = (x - 2 * y,x + z - y,y) = (0,0,0). 1117*) 1118Theorem zagier_0_y_0: 1119 !y. zagier (0,y,0) = (0,0,y) 1120Proof 1121 rw[zagier_def] 1122QED 1123 1124(* Theorem: zagier (0,0,z) = (0,z,0) *) 1125(* Proof: 1126 If x = 0 and y = 0, x < y - z becomes 0 < 0 - z = 0 is false. 1127 Next x < 2 * y becomes 0 < 2 * 0 = 0 is false, too. 1128 Thus zagier (0,0,z) = (x - 2 * y,x + z - y,y) = (0,z,0). 1129*) 1130Theorem zagier_0_0_z: 1131 !z. zagier (0,0,z) = (0,z,0) 1132Proof 1133 simp[zagier_def] 1134QED 1135 1136(* Theorem: zagier (0,0,0) = (0,0,0) *) 1137(* Proof: by zagier_x_0_0, or zagier_0_0_z or zagier_0_y_0 *) 1138Theorem zagier_0_0_0: 1139 zagier (0,0,0) = (0,0,0) 1140Proof 1141 simp[zagier_def] 1142QED 1143 1144(* Theorem: zagier (x, y, y) = 1145 if x < 2 * y then (2 * y - x,y,x) else (x - 2 * y,x,y) *) 1146(* Proof: by zagier_def, x < y - y = 0 is false. *) 1147Theorem zagier_x_y_y: 1148 !x y. zagier (x, y, y) = 1149 if x < 2 * y then (2 * y - x,y,x) else (x - 2 * y,x,y) 1150Proof 1151 rw[zagier_def] 1152QED 1153 1154(* Theorem: x <> 0 ==> (zagier (x,y,z) = (x,y,z) <=> x = y) *) 1155(* Proof: 1156 By zagier_def, 1157 If x < y - z, then 0 < y by x <> 0. 1158 (x + 2 * z,z,y - z - x) = (x, y, z) 1159 <=> x + 2 * z = x, z = y, y - z - x = z 1160 <=> x + 2 * y = x, z = y, 0 - x = y 1161 <=> x + 2 * y = x, z = y, y = 0 1162 This contradicts 0 < y. 1163 Next, if x < 2 * y, 1164 (2 * y - x,y,x + z - y) = (x, y, z) 1165 <=> 2 * y - x = x, y = y, x + z - y = z 1166 <=> 2 * y - x = x, y = y, x - y = 0 1167 <=> x = y 1168 Otherwise, 1169 (x - 2 * y,x + z - y,y) = (x, y, z) 1170 <=> x - 2 * y = x, x + z - y = y, y = z 1171 <=> x - 2 * y = x, x = y, y = z 1172 <=> x - 2 * x = x, x = y, y = z 1173 <=> x = 0, y = 0, z = 0 1174 This contradicts x <> 0. 1175*) 1176Theorem zagier_fix: 1177 !x y z. x <> 0 ==> (zagier (x,y,z) = (x,y,z) <=> x = y) 1178Proof 1179 rw[zagier_def] 1180QED 1181 1182(* Theorem: x <> 0 ==> zagier (x, x, z) = (x, x, z) *) 1183(* Proof: by zagier_fix. *) 1184Theorem zagier_x_x_z: 1185 !x z. x <> 0 ==> zagier (x, x, z) = (x, x, z) 1186Proof 1187 simp[zagier_fix] 1188QED 1189 1190(* Theorem: (x, y, z) IN mills n ==> zagier (x, y, z) IN mills n *) 1191(* Proof: 1192 By mills_def, windmill_def, zagier_def, this is to show: 1193 (1) x < y - z ==> 1194 4 * (y * z) + x ** 2 = 4 * (z * (y - (x + z))) + (x + 2 * z) ** 2 1195 or windmill x y z = windmill (x + 2 * z) z (y - (x + z)) 1196 From x < y - z 1197 so x + z < y by LESS_SUB_ADD_LESS, or in detail: 1198 Note x < 0 is impossible, so ~(y <= z), or z < y, implies z <= y. 1199 so x + z < y - z + z = y by SUB_ADD, z <= y 1200 Thus x + z <= y by x + z < y 1201 1202 windmill (x + 2 * z) z (y - (x + z)) 1203 = (x + 2 * z) ** 2 + 4 * z * (y - (x + z)) 1204 = x ** 2 + (2 * z) ** 2 + 4 * x * z + 4 * z * (y - (x + z)) by binomial_add 1205 = x ** 2 + 4 * z * z + 4 * z * x + 4 * z (y - (x + z)) by EXP_2 1206 = x ** 2 + 4 * z * (x + z) + 4 * z (y - (x + z)) by LEFT_ADD_DISTRIB 1207 = x ** 2 + 4 * z * ((x + z) + (y - (x + z))) by LEFT_ADD_DISTRIB 1208 = x ** 2 + 4 * z * ((x + z) + y - (x + z)) by LESS_EQ_ADD_SUB, x + z <= y 1209 = x ** 2 + 4 * z * (y + (x + z) - (x + z)) 1210 = x ** 2 + 4 * z * y by ADD_SUB 1211 = x ** 2 + 4 * y * z 1212 = windmill x y z 1213 1214 (2) ~(x < y - z) /\ x < 2 * y ==> 1215 4 * (y * z) + x ** 2 = 4 * (y * (x + z - y)) + (2 * y - x) ** 2 1216 or windmill x y z = windmill (2 * y - x) y (x + z - y) 1217 Note y - z <= x by ~(x < y - z) 1218 so y <= x + z by SUB_RIGHT_LESS_EQ 1219 1220 windmill (2 * y - x) y (x + z - y) 1221 = (2 * y - x) ** 2 + 4 * y * (x + z - y) 1222 = (2 * y - x) ** 2 + 8 * y * x - 8 * y * x + 4 * y * (x + z - y) by ADD_SUB 1223 = (2 * y - x) ** 2 + 8 * y * x + 4 * y * (x + z - y) - 8 * y * x by SUB_RIGHT_ADD, 1224 since 8 * y * x <= (2 * y - x) ** 2 + 8 * y * x 1225 = (2 * y + x) ** 2 + 4 * y * (x + z - y) - 8 * y * x by binomial_sub_add, x < 2 * y 1226 = (2 * y) ** 2 + x ** 2 + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x 1227 by binomial_add 1228 = x ** 2 + 4 * y * y + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x 1229 by EXP_2 1230 = x ** 2 + 4 * y * (y + x) + 4 * y * (x + z - y) - 8 * y * x 1231 by LEFT_ADD_DISTRIB 1232 = x ** 2 + 4 * y * (y + x + (x + z - y)) - 8 * y * x by LEFT_ADD_DISTRIB 1233 = x ** 2 + 4 * y * (y + x + x + z - y) - 8 * y * x by LESS_EQ_ADD_SUB, y <= x + z 1234 = x ** 2 + 4 * y * (2 * x + z) - 4 * y * (2 * x) by arithmetic 1235 = x ** 2 + 4 * y * (2 * x + z - 2 * x) by LEFT_SUB_DISTRIB 1236 = x ** 2 + 4 * y * z by ADD_SUB 1237 = windmill x y z 1238 1239 (3) ~(x < y - z) /\ ~(x < 2 * y) ==> 1240 4 * (y * z) + x ** 2 = 4 * (y * (x + z - y)) + (x - 2 * y) ** 2 1241 or windmill x y z = windmill (x - 2 * y) y (x + z - y) 1242 Note y - z <= x by ~(x < y - z) 1243 so y <= x + z by SUB_RIGHT_LESS_EQ 1244 Also 2 * y <= x by ~(x < 2 * y) 1245 1246 windmill (x - 2 * y) y (x + z - y) 1247 = (x - 2 * y) ** 2 + 4 * y * (x + z - y) 1248 = (x - 2 * y) ** 2 + 8 * y * x - 8 * y * x + 4 * y * (x + z - y) by ADD_SUB 1249 = (x - 2 * y) ** 2 + 8 * y * x + 4 * y * (x + z - y) - 8 * y * x by SUB_RIGHT_ADD, 1250 since 8 * y * x <= (2 * y - x) ** 2 + 8 * y * x 1251 = (x + 2 * y) ** 2 + 4 * y * (x + z - y) - 8 * y * x by binomial_sub_add, 2 * y <= x 1252 = x ** 2 + (2 * y) ** 2 + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x 1253 by binomial_add 1254 = x ** 2 + 4 * y * y + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x 1255 by EXP_2 1256 = x ** 2 + 4 * y * (y + x) + 4 * y * (x + z - y) - 8 * y * x 1257 by LEFT_ADD_DISTRIB 1258 = x ** 2 + 4 * y * (y + x + (x + z - y)) - 8 * y * x by LEFT_ADD_DISTRIB 1259 = x ** 2 + 4 * y * (y + x + x + z - y) - 8 * y * x by LESS_EQ_ADD_SUB, y <= x + z 1260 = x ** 2 + 4 * y * (2 * x + z) - 4 * y * (2 * x) by arithmetic 1261 = x ** 2 + 4 * y * (2 * x + z - 2 * x) by LEFT_SUB_DISTRIB 1262 = x ** 2 + 4 * y * z by ADD_SUB 1263 = windmill x y z 1264*) 1265Theorem zagier_closure: 1266 !n x y z. (x, y, z) IN mills n ==> zagier (x, y, z) IN mills n 1267Proof 1268 rw[mills_def, windmill_def, zagier_def] >| [ 1269 `x + z < y` by decide_tac >> 1270 `(x + 2 * z) ** 2 + 4 * (z * (y - (x + z))) = 1271 x ** 2 + (2 * z) ** 2 + 4 * x * z + 4 * z * (y - (x + z))` by simp[binomial_add] >> 1272 `_ = x ** 2 + (2 * z) ** 2 + 4 * z * x + 4 * z * (y - (x + z))` by decide_tac >> 1273 `_ = x ** 2 + (2 * z) * (2 * z) + 4 * z * x + 4 * z * (y - (x + z))` by metis_tac[EXP_2] >> 1274 `_ = x ** 2 + 4 * z * z + 4 * z * x + 4 * z * (y - (x + z))` by decide_tac >> 1275 `_ = x ** 2 + 4 * z * (z + x) + 4 * z * (y - (x + z))` by decide_tac >> 1276 `_ = x ** 2 + 4 * z * ((x + z) + (y - (x + z)))` by rw[LEFT_ADD_DISTRIB] >> 1277 `_ = x ** 2 + 4 * z * ((x + z) + y - (x + z))` by rw[] >> 1278 `_ = x ** 2 + 4 * z * y` by decide_tac >> 1279 simp[], 1280 `y <= x + z` by decide_tac >> 1281 `(2 * y - x) ** 2 + 4 * (y * (x + z - y)) = 1282 (2 * y - x) ** 2 + 4 * y * (x + z - y)` by decide_tac >> 1283 `_ = (2 * y - x) ** 2 + 8 * y * x - 8 * y * x + 4 * y * (x + z - y)` by decide_tac >> 1284 `_ = (2 * y - x) ** 2 + 4 * (2 * y) * x + 4 * y * (x + z - y) - 8 * y * x` by fs[] >> 1285 `_ = (2 * y + x) ** 2 + 4 * y * (x + z - y) - 8 * y * x` by simp[binomial_sub_add] >> 1286 `_ = (2 * y) ** 2 + x ** 2 + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x` by simp[binomial_add] >> 1287 `_ = (2 * y) * (2 * y) + x ** 2 + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x` by metis_tac[EXP_2] >> 1288 `_ = x ** 2 + 4 * y * y + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x` by fs[] >> 1289 `_ = x ** 2 + 4 * y * (y + x) + 4 * y * (x + z - y) - 8 * y * x` by rw[LEFT_ADD_DISTRIB] >> 1290 `_ = x ** 2 + 4 * y * (y + x + (x + z - y)) - 8 * y * x` by rw[LEFT_ADD_DISTRIB] >> 1291 `_ = x ** 2 + 4 * y * (y + x + x + z - y) - 8 * y * x` by rw[LESS_EQ_ADD_SUB] >> 1292 `_ = x ** 2 + 4 * y * (x + x + z) - 4 * y * (x + x)` by decide_tac >> 1293 `_ = x ** 2 + 4 * y * (x + x + z - (x + x))` by decide_tac >> 1294 `_ = x ** 2 + 4 * y * z` by decide_tac >> 1295 simp[], 1296 `y <= x + z` by decide_tac >> 1297 `2 * y <= x` by decide_tac >> 1298 `(x - 2 * y) ** 2 + 4 * (y * (x + z - y)) = 1299 (x - 2 * y) ** 2 + 4 * y * (x + z - y)` by decide_tac >> 1300 `_ = (x - 2 * y) ** 2 + 8 * y * x - 8 * y * x + 4 * y * (x + z - y)` by decide_tac >> 1301 `_ = (x - 2 * y) ** 2 + 4 * x * (2 * y) + 4 * y * (x + z - y) - 8 * y * x` by fs[] >> 1302 `_ = (x + 2 * y) ** 2 + 4 * y * (x + z - y) - 8 * y * x` by simp[binomial_sub_add] >> 1303 `_ = x ** 2 + (2 * y) ** 2 + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x` by simp[binomial_add] >> 1304 `_ = x ** 2 + (2 * y) * (2 * y) + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x` by metis_tac[EXP_2] >> 1305 `_ = x ** 2 + 4 * y * y + 4 * y * x + 4 * y * (x + z - y) - 8 * y * x` by fs[] >> 1306 `_ = x ** 2 + 4 * y * (y + x) + 4 * y * (x + z - y) - 8 * y * x` by rw[LEFT_ADD_DISTRIB] >> 1307 `_ = x ** 2 + 4 * y * (y + x + (x + z - y)) - 8 * y * x` by rw[LEFT_ADD_DISTRIB] >> 1308 `_ = x ** 2 + 4 * y * (y + x + x + z - y) - 8 * y * x` by rw[LESS_EQ_ADD_SUB] >> 1309 `_ = x ** 2 + 4 * y * (x + x + z) - 4 * y * (x + x)` by decide_tac >> 1310 `_ = x ** 2 + 4 * y * (x + x + z - (x + x))` by decide_tac >> 1311 `_ = x ** 2 + 4 * y * z` by decide_tac >> 1312 simp[] 1313 ] 1314QED 1315 1316 1317(* Theorem: t IN mills n ==> zagier t IN mills n *) 1318(* Proof: by zagier_closure *) 1319Theorem zagier_closure_alt: 1320 !n t. t IN mills n ==> zagier t IN mills n 1321Proof 1322 metis_tac[triple_parts, zagier_closure] 1323QED 1324 1325 1326(* Theorem: x <> 0 /\ y <> 0 /\ z <> 0 ==> zagier (zagier (x, y, z)) = (x, y, z) *) 1327(* Proof: 1328 Put the 3 branches of zagier_def into 5 cases for a triple (x,y,z): 1329 Case 1: x < y and x < y - z for branch 1 1330 Case 2: x < y and y - z < x for branch 2 1331 Case 3: x = y for branch 2 1332 Case 4: y < x and x < 2 * y for branch 2 1333 Case 5: y < x and 2 * y < x for branch 3 1334 1335 Note x <> 0 /\ y <> 0 /\ z <> 0 implies proper windmills. 1336 Boundary cases correspond to improper windmills, hence not considered: 1337 For x < y, x = y - z gives 1338 n = (y - z) ** 2 + 4 * y * z = (y + z) ** 2, a windmill without arms. 1339 For y < x, x = 2 * y gives 1340 n = (2 * y) ** 2 + 4 * y * z = 4 * y * (y + z), a windmill with only four arms. 1341 1342 By zagier_def, 1343 If x < y - z, 1344 then (x,y,z) will be in case 1. 1345 zagier (x,y,z) = (x + 2 * z,z,y - z - x) = (x',y',z') by first branch 1346 so x' = x + 2 * z = z + x + z, expand inner square, 1347 y' = z, 1348 z' = y - z - x, excess to fit the mind. 1349 Thus 2 * y' = 2 * z < x + 2 * z = x' for x <> 0, 1350 which is case 5, so 1351 zagier (x',y',z') 1352 = (x' - 2 * y',x' + z' - y',y') by third branch 1353 = (x + 2 * z - 2 * z, x + 2 * z + (y - z - x) - z, z) 1354 = (x, y, z) 1355 1356 If ~(x < y - z), but x < 2 * y, 1357 then (x,y,z) will be in case 2 or case 3 or case 4. 1358 zagier (x,y,z) = (2 * y - x,y,x + z - y) = (x',y',z') by second branch 1359 so x' = 2 * y - x, 1360 y' = y, 1361 z' = x + z - y. 1362 1363 If x < y, this is case 2. 1364 x' = 2 * y - x = x + 2 * (y - x), expand inner square, 1365 y' = y, 1366 z' = x + z - y, excess to fit the mind. 1367 Thus x' = 2 * y - x < 2 * y = 2 * y' for x <> 0, 1368 which is case 4, so 1369 zagier (x',y',z') 1370 = (2 * y' - x',y',x' + z' - y') by second branch 1371 = (2 * y - (2 * y - x),y,(2 * y - x) + (x + z - y) - y) 1372 = (x, y, z) 1373 1374 If x = y, this is case 3. 1375 x' = 2 * y - x = 2 * y - y = x, no change to inner square, 1376 y' = y, 1377 z' = x + z - y = y + z - y = z, no excess, fit the mind. 1378 Thus x' = y', which is case 3, so 1379 zagier (x',y',z') 1380 = (2 * y' - x',y',x' + z' - y') by second branch 1381 = (2 * y - (2 * y - x),y,(2 * y - x) + (x + z - y) - y) 1382 = (x, y, z) 1383 1384 If y < x, this is case 4. 1385 x' = 2 * y - x = x - 2 * (x - y), shrink inner square, 1386 y' = y, 1387 z' = x + z - y = z + (x - y), excess to fit the mind. 1388 Thus y' - z' = y - (x + z - y) < 2 * y - x = x' for z <> 0, 1389 which is case 2, so 1390 zagier (x',y',z') 1391 = (2 * y' - x',y',x' + z' - y') by second branch 1392 = (2 * y - (2 * y - x),y,(2 * y - x) + (x + z - y) - y) 1393 = (x, y, z) 1394 1395 If ~(x < y - z), and ~(x < 2 * y), 1396 then (x,y,z) will be in case 5. 1397 zagier (x,y,z) = (x - 2 * y,x + z - y,y) = (x',y',z') by third branch 1398 so x' = x - 2 * y, shrink inner square 1399 y' = x + z - y, 1400 z' = y. 1401 Thus x' = x - 2 * y < (x + z - y) - y = y' - z' for z <> 0, 1402 which is case 1, so 1403 zagier (x',y',z') 1404 = (x' + 2 * z',z',y' - z' - x') by first branch 1405 = (x - 2 * y + 2 * y, y, (x + z - y) - y - (x - 2 * y)) 1406 = (x, y, z) 1407*) 1408Theorem zagier_involute: 1409 !x y z. x <> 0 /\ y <> 0 /\ z <> 0 ==> 1410 zagier (zagier (x, y, z)) = (x, y, z) 1411Proof 1412 rw[zagier_def] 1413QED 1414 1415(* Theorem: x * y * z <> 0 ==> zagier (zagier (x, y, z)) = (x, y, z) *) 1416(* Proof: by MULT3_EQ_0, zagier_involute. *) 1417Theorem zagier_involute_xyz: 1418 !x y z. x * y * z <> 0 ==> zagier (zagier (x, y, z)) = (x, y, z) 1419Proof 1420 rw[zagier_involute] 1421QED 1422 1423(* Theorem: FST t <> 0 /\ FST (SND t) <> 0 /\ SND (SND t) <> 0 ==> 1424 zagier (zagier t) = t *) 1425(* Proof: by zagier_involute *) 1426Theorem zagier_involute_thm: 1427 !t. FST t <> 0 /\ FST (SND t) <> 0 /\ SND (SND t) <> 0 ==> 1428 zagier (zagier t) = t 1429Proof 1430 simp[zagier_involute, FORALL_PROD] 1431QED 1432 1433(* Define cuboid of a triple *) 1434Definition cuboid_def: 1435 cuboid (x, y, z) = x * y * z 1436End 1437 1438(* Theorem: cuboid (x, y, z) = 0 <=> x = 0 \/ y = 0 \/ z = 0 *) 1439(* Proof: by cuboid_def, MULT_EQ_0 *) 1440Theorem cuboid_eq_0: 1441 !x y z. cuboid (x, y, z) = 0 <=> x = 0 \/ y = 0 \/ z = 0 1442Proof 1443 simp[cuboid_def] 1444QED 1445 1446(* Theorem: cuboid t <> 0 ==> zagier (zagier t) = t *) 1447(* Proof: 1448 Let t = (x, y, z). 1449 Then x <> 0 /\ y <> 0 /\ z <> 0 by cuboid_eq_0 1450 Thus zagier (zagier t) = t by zagier_involute 1451*) 1452Theorem zagier_involute_alt: 1453 !t. cuboid t <> 0 ==> zagier (zagier t) = t 1454Proof 1455 metis_tac[triple_parts, cuboid_eq_0, zagier_involute] 1456QED 1457 1458(* Theorem: ~square n /\ n MOD 4 <> 0 ==> zagier involute (mills n) *) 1459(* Proof: 1460 Let t = (x,y,z). by triple_parts 1461 Then zagier (x,y,z) IN mills n by zagier_closure 1462 Note !x y z. (x,y,z) IN mills n ==> 1463 x <> 0 /\ y <> 0 /\ z <> 0 by mills_triple_nonzero 1464 so zagier (zagier (x,y,z)) = (x,y,z) by zagier_involute 1465 so zagier involute (mills n) by involution 1466*) 1467Theorem zagier_involute_mills: 1468 !n. ~square n /\ n MOD 4 <> 0 ==> zagier involute (mills n) 1469Proof 1470 metis_tac[mills_triple_nonzero, zagier_closure, zagier_involute, triple_parts] 1471QED 1472 1473(* ------------------------------------------------------------------------- *) 1474(* Part 3: Mind of a windmill *) 1475(* ------------------------------------------------------------------------- *) 1476 1477(* The mind of a windmill: 1478 1479 +----+ 1480 | | y 1481 | +-------------+ 1482 | | |z 1483 | +-----+----+--+ 1484 | | x | z | 1485 | | | | 1486 +--+----+-----+ | 1487 | | | 1488 +-------------+ | 1489 | | 1490 +----+ 1491 Algorithm: 1492 1. draw the square of side x. 1493 2. from each square vertex, put the line y alongside, in clockwise manner. 1494 3. complete the 4 rectangles y * z, around the central square. 1495 1496 "mind" is the side of central square: 1497 from the picture above, it is evident that, 1498 assume y > z, then (y - z) is the slack for inner square x to fit in. 1499 If x < y - z, the inner square can grow to: 1500 x' = x + 2 * z, y' = z, z' = y - x - z. (result looks like the third diagram) 1501 1502 y 1503 +--------+ 1504 +----| |z 1505 | +-----+--+-+ 1506 | | x | z | 1507 | | | | 1508 +-+--+-----+ | 1509 | |----+ 1510 +--------+ 1511 1512 If y < z, assume y > x, the inner square can only grow to: 1513 x' = x + 2 * (y - x) = 2 * y - x, y' = y, z' = x + z - y. 1514 this only works with 0 < x' for x < 2 * y. 1515 1516 Otherwise, y < x, the mind has to shrink: 1517 x' = x - 2 * (x - y) = 2 * y - x, y' = x + z - y, z' = y. 1518 1519 y 1520 +---+ 1521 | |z 1522 +---+-+---+ 1523 +---| x | z | 1524 | | |---+ 1525 +---+-+---+ 1526 | | 1527 +---+ 1528 1529*) 1530(* Define the mind of a triple, for a windmill. *) 1531Definition mind_def: 1532 mind (x,y,z) = if x < y - z then x + 2 * z 1533 else if x < y then 2 * y - x 1534 else x (* maximal central square *) 1535End 1536 1537(* Investigation: 1538 1539zagier (3,8,2) = (7,2,3) 1540zagier (3,6,4) = (9,6,1) 1541zagier (7,3,2) = (1,6,3) 1542 1543EVAL ``mind (3,8,2)``; = 7 1544EVAL ``mind (zagier (3,8,2))``; = 7 1545 1546EVAL ``mind (3,6,4)``; = 9 1547EVAL ``mind (zagier (3,6,4))``; = 9 1548EVAL ``mind (9,6,1)``; = 9 1549 1550EVAL ``mind (7,3,2)``; = 7 1551EVAL ``mind (zagier (7,3,2))``; = 7 1552 1553*) 1554 1555 1556(* Theorem: mind (zagier (x, y, z)) = mind (x, y, z) *) 1557(* Proof: 1558 If x < y - z, 1559 then since 2 * z < z is false, 1560 so x + 2 * z < z - (y - z - x) is false 1561 and x + 2 * z < z is also false, so 1562 mind (zagier (x, y, z)) 1563 = mind (x + 2 * z,z,y - z - x) by zagier_def 1564 = x + 2 * z by mind_def, third case 1565 = mind (x, y, z) by mind_def 1566 Otherwise, if x < y, 1567 then x < 2 * y, 1568 mind (zagier (x, y, z)) 1569 = mind (2 * y - x,y,x + z - y) by zagier_def 1570 = 2 * y - x by mind_def, third case 1571 = mind (x, y, z) by mind_def 1572 Otherwise, if x < 2 * y, 1573 then if 2 * y - x < y 1574 mind (zagier (x, y, z)) 1575 = mind (2 * y - x,y,x + z - y) by zagier_def 1576 = 2 * y - (2 * y - x) by mind_def, second case 1577 = x 1578 = mind (x, y, z) by mind_def 1579 else ~(2 * y - x < y) 1580 so y = x, 1581 mind (zagier (x, y, z)) 1582 = mind (2 * y - x,y,x + z - y) by zagier_def 1583 = 2 * y - x by mind_def, third case 1584 = x by y = x 1585 = mind (x, y, z) by mind_def 1586 Otherwise, 1587 then if x - 2 * y < (x + z - y) - y 1588 mind (zagier (x, y, z)) 1589 = mind (x - 2 * y,x + z - y,y) by zagier_def 1590 = (x - 2 * y) + 2 * y by mind_def, first case 1591 = x 1592 = mind (x, y, z) by mind_def 1593 else if x - 2 * y < x + z - y 1594 so z = 0, by arithmetic 1595 mind (zagier (x, y, z)) 1596 = mind (x - 2 * y,x + z - y,y) by zagier_def 1597 = 2 * (x + z - y) - (x - 2 * y) by mind_def, second case 1598 = 2 * x + 2 * z - x by arithmetic 1599 = x by z = 0 1600 = mind (x, y, z) by mind_def 1601 else ~(x - 2 * y < x + z - y) 1602 so y = 0, by arithmetic 1603 mind (zagier (x, y, z)) 1604 = mind (x - 2 * y,x + z - y,y) by zagier_def 1605 = x - 2 * y by mind_def, third case 1606 = x by y = 0 1607 = mind (x, y, z) by mind_def 1608*) 1609Theorem mind_zagier_eqn: 1610 !x y z. mind (zagier (x, y, z)) = mind (x, y, z) 1611Proof 1612 rw[mind_def, zagier_def] 1613QED 1614 1615(* Theorem: mind (zagier t) = mind t *) 1616(* Proof: by triple_parts, mind_zagier_eqn. *) 1617Theorem mind_zagier_thm: 1618 !t. mind (zagier t) = mind t 1619Proof 1620 metis_tac[triple_parts, mind_zagier_eqn] 1621QED 1622 1623(* What is the mind of a flip? *) 1624 1625(* Theorem: mind (flip (x,y,z)) = 1626 if x < z - y then x + 2 * y 1627 else if x < z then 2 * z - x 1628 else x *) 1629(* Proof: by mind_def, flip_def *) 1630Theorem mind_flip_eqn: 1631 !x y z. mind (flip (x,y,z)) = 1632 if x < z - y then x + 2 * y 1633 else if x < z then 2 * z - x 1634 else x 1635Proof 1636 simp[mind_def, flip_def] 1637QED 1638 1639(* 1640> mind_flip_eqn |> SPEC ``1`` |> SPEC ``1`` |> SPEC ``k:num``; 1641val it = |- mind (flip (1,1,k)) = 1642if 1 < k - 1 then 1 + 2 * 1 else if 1 < k then 2 * k - 1 else 1: thm 1643*) 1644 1645(* Theorem: mind (flip (1,1,z)) = if z < 2 then 1 else 3 *) 1646(* Proof: by mind_flip_eqn. *) 1647Theorem mind_flip_1_1_z: 1648 !z. mind (flip (1,1,z)) = if z < 2 then 1 else 3 1649Proof 1650 simp[mind_flip_eqn] 1651QED 1652 1653(* 1654> mind_flip_eqn |> SPEC ``x:num`` |> SPEC ``x:num`` |> SPEC ``z:num``; 1655val it = |- mind (flip (x,x,z)) = 1656if x < z - x then x + 2 * x else if x < z then 2 * z - x else x: thm 1657*) 1658 1659(* Theorem: mind (flip (x,y,y)) = if x < y then 2 * y - x else x *) 1660(* Proof: by mind_flip_eqn. *) 1661Theorem mind_flip_x_x_z: 1662 !x z. mind (flip (x,x,z)) = 1663 if x < z - x then 3 * x else if x < z then 2 * z - x else x 1664Proof 1665 simp[mind_flip_eqn] 1666QED 1667 1668(* 1669> mind_flip_eqn |> SPEC ``x:num`` |> SPEC ``y:num`` |> SPEC ``y:num``; 1670val it = |- mind (flip (x,y,y)) = 1671if x < y - y then x + 2 * y else if x < y then 2 * y - x else x: thm 1672*) 1673 1674(* Theorem: mind (flip (x,y,y)) = if x < y then 2 * y - x else x *) 1675(* Proof: by mind_flip_eqn. *) 1676Theorem mind_flip_x_y_y: 1677 !x y. mind (flip (x,y,y)) = if x < y then 2 * y - x else x 1678Proof 1679 simp[mind_flip_eqn] 1680QED 1681 1682(* ------------------------------------------------------------------------- *) 1683(* Gap of a Windmill. *) 1684(* ------------------------------------------------------------------------- *) 1685 1686(* Define the gap, the absolute difference of y and z in a triple. *) 1687Definition gap_def: 1688 gap (x,y,z) = if y < z then z - y else y - z 1689End 1690 1691(* Theorem: gap (flip (x, y, z)) = gap (x, y, z) *) 1692(* Proof: by gap_def, flip_def. *) 1693Theorem gap_flip_eqn: 1694 !x y z. gap (flip (x, y, z)) = gap (x, y, z) 1695Proof 1696 simp[gap_def, flip_def] 1697QED 1698 1699(* Theorem: gap (flip t) = gap t *) 1700(* Proof: by triple_parts, gap_flip_eqn. *) 1701Theorem gap_flip_thm: 1702 !t. gap (flip t) = gap t 1703Proof 1704 metis_tac[triple_parts, gap_flip_eqn] 1705QED 1706 1707 1708(* ------------------------------------------------------------------------- *) 1709(* Zagier and Flip. *) 1710(* ------------------------------------------------------------------------- *) 1711 1712(* Theorem: (zagier o flip) (1,1,z) = 1713 if z = 0 then (1,2,0) 1714 else if z = 1 then (1,1,1) 1715 else if z = 2 then (3,2,0) 1716 else (3,1,z - 2) *) 1717(* Proof: 1718 (zagier o flip) (1,1,z) 1719 = zagier (flip (1,1,z)) by o_THM 1720 = zagier (1,z,1) by flip_def 1721 = if z = 0 then (1,2,0) 1722 else if z = 1 then (1,1,1) 1723 else if z = 2 then (3,2,0) 1724 else (3,1,z - 2) by zagier_def 1725 1726> zagier_def |> ISPEC ``1`` |> ISPEC ``z:num`` |> ISPEC ``1`` |> SIMP_RULE arith_ss []; 1727|- zagier (1,z,1) = 1728 if 1 < z - 1 then (3,1,z - 2) 1729 else if 1 < 2 * z then (2 * z - 1,z,2 - z) 1730 else (1 - 2 * z,2 - z,z): thm 1731zagier_def |> ISPEC ``x:num`` |> ISPEC ``z:num`` |> ISPEC ``x:num`` |> SIMP_RULE arith_ss []; 1732*) 1733Theorem zagier_flip_1_1_z: 1734 !z. (zagier o flip) (1,1,z) = 1735 if z = 0 then (1,2,0) 1736 else if z = 1 then (1,1,1) 1737 else if z = 2 then (3,2,0) 1738 else (3,1,z - 2) 1739Proof 1740 rw[zagier_def, flip_def] 1741QED 1742 1743(* ------------------------------------------------------------------------- *) 1744(* Computation of (mills n) *) 1745(* ------------------------------------------------------------------------- *) 1746 1747(* Generate the tuples *) 1748Definition tuples_helper_def: 1749 tuples_helper k 0 = [] /\ 1750 tuples_helper k (SUC n) = 1751 ZIP ((GENLIST (K (SUC n)) k), (GENLIST SUC k)) ++ tuples_helper k n 1752End 1753 1754Definition tuples_def: 1755 tuples k = tuples_helper k k 1756End 1757(* 1758EVAL ``tuples_helper 3 3``; 1759EVAL ``tuples 3``; 1760[(3,1); (3,2); (3,3); (2,1); (2,2); (2,3); (1,1); (1,2); (1,3)] 1761*) 1762 1763(* Generate the triples *) 1764Definition triples_helper_def: 1765 triples_helper k 0 = [] /\ 1766 triples_helper k (SUC n) = 1767 ZIP ((GENLIST (K (SUC n)) (k * k)), tuples k) ++ triples_helper k n 1768End 1769 1770Definition triples_def: 1771 triples k = triples_helper k k 1772End 1773 1774(* 1775EVAL ``triples_helper 3 3``; 1776EVAL ``triples 3``; 1777*) 1778 1779(* 1780EVAL ``FILTER (\(x, y, z). 5 = x * x + 4 * y * z) (triples 5)``; -> [(1,1,1)] 1781EVAL ``FILTER (\(x, y, z). 13 = x * x + 4 * y * z) (triples 13)``; -> [(3,1,1); (1,3,1); (1,1,3)] 1782*) 1783 1784(* Compute (mills n) *) 1785Definition mills_of_def: 1786 mills_of n = FILTER (\(x, y, z). n = windmill x y z) (triples n) 1787End 1788 1789(* 1790EVAL ``mills_of 5``; [(1,1,1)] 1791EVAL ``MAP zagier (mills_of 5)``; [(1,1,1)] 1792EVAL ``mills_of 13``; [(3,1,1); (1,3,1); (1,1,3)] 1793EVAL ``MAP zagier (mills_of 13)``; [(1,3,1); (3,1,1); (1,1,3)] 1794*) 1795 1796(* 1797 1798EVAL ``MAP2 (\x y. (x,y)) (GENLIST (K 1) 5) (GENLIST SUC 5)``; 1799EVAL ``ZIP ((GENLIST (K 1) 3), (GENLIST SUC 3))``; 1800 1801EVAL ``(count 13) CROSS (count 13) CROSS (count 13)``; 1802EVAL ``13 * 13 * 13``; -> 2197 1803 1804EVAL ``count 5 CROSS (count 5 CROSS count 5)``; 1805EVAL ``ZIP (GENLIST SUC 5, GENLIST SUC 5)``; 1806EVAL ``ZIP (GENLIST SUC 5, ZIP (GENLIST SUC 5, GENLIST SUC 5))``; 1807 1808EVAL ``FILTER (\(x, y, z). 5 = x * x + 4 * y * z) (ZIP (GENLIST SUC 5, ZIP (GENLIST SUC 5, GENLIST SUC 5)))``; 1809 1810(* Generate the triples *) 1811val triples_def = Define` 1812 triples n = FILTER (\(x, y, z). n = x * x + 4 * y * z) 1813 (ZIP (GENLIST SUC n, ZIP (GENLIST SUC n, GENLIST SUC n))) 1814`; 1815 1816EVAL ``triples 5``; -> [(1,1,1)] 1817EVAL ``triples 13``; the list is wrong! 1818 1819*) 1820 1821(* ------------------------------------------------------------------------- *) 1822 1823(* export theory at end *) 1824val _ = export_theory(); 1825 1826(*===========================================================================*) 1827