<< problem 501 - Eight Divisors | Tangent Circles - problem 510 >> |
Problem 504: Square on the Inside
(see projecteuler.net/problem=504)
Let ABCD be a quadrilateral whose vertices are lattice points lying on the coordinate axes as follows:
A(a, 0), B(0, b), C(-c, 0), D(0, -d), where 1 <= a, b, c, d <= m and a, b, c, d, m are integers.
It can be shown that for m = 4 there are exactly 256 valid ways to construct ABCD.
Of these 256 quadrilaterals, 42 of them strictly contain a square number of lattice points.
How many quadrilaterals ABCD strictly contain a square number of lattice points for m = 100?
My Algorithm
Pick's theorem states A = i + b/2 - 1 (see en.wikipedia.org/wiki/Pick's_theorem):
if a simple polygon has i interior integer points and b boundary integer points, then A must be a perfect square to match the conditions of the problem.
So all I have to do is finding i and b - that's what my function countLatticePoints
is for.
The bounding box that enclosed the quadrilateral and is parallel to the x- and y-axis has an area of
ab + bc + cd + ad = a(b + d) + c(b + d) = (a + c)(b + d)
The four triangles cover only half the area of the bounding box. My variable inside
equals i of the formula:
i = dfrac{1}{2} (a + c)(b + d)
The boundary integer points can be subdivided into two groups:
- 4 corners A(a, 0), B(0, b), C(-c, 0), D(0, -d)
- the edges between A and B, B and C, C and D, D and A
And there are gcd(x,y) + 1 such integer lattice points: the origin plus all multiples of k * ( dfrac{x}{gcd(x,y)}, dfrac{y}{gcd(x,y)} ).
However, each corner appears twice: A in AB and DA, B in AB and BC, ...
Using just gcd(x,y) (without minus one) solves this problem.
Unfortunately, my results of
countLatticePoints
were off by 4:the corners A, B, C and D are not part of the original bounding-box and must be subtracted, too.
Alternative Approaches
There are many symmetries. For example I could assume a <= b and reduce the search space accordingly.
Note
I set up a few caches to speed up the computations considerably:
- the same GCDs have to by computed quite many times, therefore I store them in
pointsOnEdge
sqrt
is pretty slow and a simple containerisSquare
stores atrue
flag for each square up to the largest area of an quadrilateral (100,100,100,100) → about 20000.
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho 4 | ./504
Output:
Note: the original problem's input 100
cannot be entered
because just copying results is a soft skill reserved for idiots.
(this interactive test is still under development, computations will be aborted after one second)
My code
… was written in C++11 and can be compiled with G++, Clang++, Visual C++. You can download it, too. Or just jump to my GitHub repository.
#include <iostream>
#include <vector>
// cache gcd(a,b)
std::vector<std::vector<unsigned int>> pointsOnEdge;
// greatest common divisor
template <typename T>
T gcd(T a, T b)
{
while (a != 0)
{
T c = a;
a = b % a;
b = c;
}
return b;
}
// count lattice points with Pick's theorem
// A(a, 0), B(b, 0), C(-c, 0), D(-d, 0)
unsigned int countLatticePoints(unsigned int a, unsigned int b, unsigned int c, unsigned int d)
{
// Pick's theorem: A = i + b/2 - 1
// where i is the area inside the quadrilateral
// and b is the number of lattice points on the boundary
// https://en.wikipedia.org/wiki/Pick%27s_theorem
// total area of the bounding box
auto rectangle = (a + c) * (b + d);
// only half of the area belongs to the quadrilateral
auto inside = rectangle / 2;
// points on the edges of the quadrilateral
auto boundary = pointsOnEdge[a][b] + pointsOnEdge[b][c] + pointsOnEdge[c][d] + pointsOnEdge[a][d];
// the four corners are not inside the bounding box (they are on its edges)
boundary -= 4;
// Pick's formula
return inside - (boundary / 2) - 1;
}
int main()
{
unsigned int limit = 100;
std::cin >> limit;
// precompute number of lattice points lying exactly on the sides/edges
pointsOnEdge.resize(limit + 1);
for (auto& x : pointsOnEdge)
x.resize(limit + 1);
for (unsigned int a = 1; a <= limit; a++)
for (unsigned int b = 1; b <= limit; b++)
// note: only one endpoint is included for each edge
// e.g. a=4 b=2 gcd(4,2)=2
// => has actually 3 lattice points:
// both endpoints plus point (3,1)
// but each endpoint is part of two edges therefore it would be counted twice
// that means that keeping only one endpoint is okay and much simpler
pointsOnEdge[a][b] = pointsOnEdge[b][a] = gcd(a, b);
// precompute square numbers
// note: std::vector<bool> would be sufficient but std::vector<char> is about 20% faster
std::vector<unsigned char> isSquare(countLatticePoints(limit, limit, limit, limit) + 1, false);
for (size_t i = 1; i*i < isSquare.size(); i++)
isSquare[i*i] = true;
// plain enumeration of all combinations, 100^4 = 10^8 iterations
auto count = 0;
for (unsigned int a = 1; a <= limit; a++)
for (unsigned int b = 1; b <= limit; b++)
for (unsigned int c = 1; c <= limit; c++)
for (unsigned int d = 1; d <= limit; d++)
// count only perfect squares
if (isSquare[countLatticePoints(a, b, c, d)])
count++;
// display result
std::cout << count << std::endl;
return 0;
}
This solution contains 11 empty lines, 25 comments and 2 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 0.22 seconds on an Intel® Core™ i7-2600K CPU @ 3.40GHz.
(compiled for x86_64 / Linux, GCC flags: -O3 -march=native -fno-exceptions -fno-rtti -std=gnu++11 -DORIGINAL
)
See here for a comparison of all solutions.
Note: interactive tests run on a weaker (=slower) computer. Some interactive tests are compiled without -DORIGINAL
.
Changelog
August 11, 2017 submitted solution
August 11, 2017 added comments
Difficulty
Project Euler ranks this problem at 15% (out of 100%).
Links
projecteuler.net/thread=504 - the best forum on the subject (note: you have to submit the correct solution first)
Code in various languages:
Python github.com/Meng-Gen/ProjectEuler/blob/master/504.py (written by Meng-Gen Tsai)
C++ github.com/evilmucedin/project-euler/blob/master/euler504/504.cpp (written by Den Raskovalov)
C++ github.com/roosephu/project-euler/blob/master/504.cpp (written by Yuping Luo)
C++ github.com/smacke/project-euler/blob/master/cpp/504.cpp (written by Stephen Macke)
Java github.com/thrap/project-euler/blob/master/src/Java/Problem504.java (written by Magnus Solheim Thrap)
Perl github.com/shlomif/project-euler/blob/master/project-euler/504/euler-504-v1.pl (written by Shlomi Fish)
Those links are just an unordered selection of source code I found with a semi-automatic search script on Google/Bing/GitHub/whatever.
You will probably stumble upon better solutions when searching on your own.
Maybe not all linked resources produce the correct result and/or exceed time/memory limits.
Heatmap
Please click on a problem's number to open my solution to that problem:
green | solutions solve the original Project Euler problem and have a perfect score of 100% at Hackerrank, too | |
yellow | solutions score less than 100% at Hackerrank (but still solve the original problem easily) | |
gray | problems are already solved but I haven't published my solution yet | |
blue | solutions are relevant for Project Euler only: there wasn't a Hackerrank version of it (at the time I solved it) or it differed too much | |
orange | problems are solved but exceed the time limit of one minute or the memory limit of 256 MByte | |
red | problems are not solved yet but I wrote a simulation to approximate the result or verified at least the given example - usually I sketched a few ideas, too | |
black | problems are solved but access to the solution is blocked for a few days until the next problem is published | |
[new] | the flashing problem is the one I solved most recently |
I stopped working on Project Euler problems around the time they released 617.
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 |
26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 |
51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 |
76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | 91 | 92 | 93 | 94 | 95 | 96 | 97 | 98 | 99 | 100 |
101 | 102 | 103 | 104 | 105 | 106 | 107 | 108 | 109 | 110 | 111 | 112 | 113 | 114 | 115 | 116 | 117 | 118 | 119 | 120 | 121 | 122 | 123 | 124 | 125 |
126 | 127 | 128 | 129 | 130 | 131 | 132 | 133 | 134 | 135 | 136 | 137 | 138 | 139 | 140 | 141 | 142 | 143 | 144 | 145 | 146 | 147 | 148 | 149 | 150 |
151 | 152 | 153 | 154 | 155 | 156 | 157 | 158 | 159 | 160 | 161 | 162 | 163 | 164 | 165 | 166 | 167 | 168 | 169 | 170 | 171 | 172 | 173 | 174 | 175 |
176 | 177 | 178 | 179 | 180 | 181 | 182 | 183 | 184 | 185 | 186 | 187 | 188 | 189 | 190 | 191 | 192 | 193 | 194 | 195 | 196 | 197 | 198 | 199 | 200 |
201 | 202 | 203 | 204 | 205 | 206 | 207 | 208 | 209 | 210 | 211 | 212 | 213 | 214 | 215 | 216 | 217 | 218 | 219 | 220 | 221 | 222 | 223 | 224 | 225 |
226 | 227 | 228 | 229 | 230 | 231 | 232 | 233 | 234 | 235 | 236 | 237 | 238 | 239 | 240 | 241 | 242 | 243 | 244 | 245 | 246 | 247 | 248 | 249 | 250 |
251 | 252 | 253 | 254 | 255 | 256 | 257 | 258 | 259 | 260 | 261 | 262 | 263 | 264 | 265 | 266 | 267 | 268 | 269 | 270 | 271 | 272 | 273 | 274 | 275 |
276 | 277 | 278 | 279 | 280 | 281 | 282 | 283 | 284 | 285 | 286 | 287 | 288 | 289 | 290 | 291 | 292 | 293 | 294 | 295 | 296 | 297 | 298 | 299 | 300 |
301 | 302 | 303 | 304 | 305 | 306 | 307 | 308 | 309 | 310 | 311 | 312 | 313 | 314 | 315 | 316 | 317 | 318 | 319 | 320 | 321 | 322 | 323 | 324 | 325 |
326 | 327 | 328 | 329 | 330 | 331 | 332 | 333 | 334 | 335 | 336 | 337 | 338 | 339 | 340 | 341 | 342 | 343 | 344 | 345 | 346 | 347 | 348 | 349 | 350 |
351 | 352 | 353 | 354 | 355 | 356 | 357 | 358 | 359 | 360 | 361 | 362 | 363 | 364 | 365 | 366 | 367 | 368 | 369 | 370 | 371 | 372 | 373 | 374 | 375 |
376 | 377 | 378 | 379 | 380 | 381 | 382 | 383 | 384 | 385 | 386 | 387 | 388 | 389 | 390 | 391 | 392 | 393 | 394 | 395 | 396 | 397 | 398 | 399 | 400 |
401 | 402 | 403 | 404 | 405 | 406 | 407 | 408 | 409 | 410 | 411 | 412 | 413 | 414 | 415 | 416 | 417 | 418 | 419 | 420 | 421 | 422 | 423 | 424 | 425 |
426 | 427 | 428 | 429 | 430 | 431 | 432 | 433 | 434 | 435 | 436 | 437 | 438 | 439 | 440 | 441 | 442 | 443 | 444 | 445 | 446 | 447 | 448 | 449 | 450 |
451 | 452 | 453 | 454 | 455 | 456 | 457 | 458 | 459 | 460 | 461 | 462 | 463 | 464 | 465 | 466 | 467 | 468 | 469 | 470 | 471 | 472 | 473 | 474 | 475 |
476 | 477 | 478 | 479 | 480 | 481 | 482 | 483 | 484 | 485 | 486 | 487 | 488 | 489 | 490 | 491 | 492 | 493 | 494 | 495 | 496 | 497 | 498 | 499 | 500 |
501 | 502 | 503 | 504 | 505 | 506 | 507 | 508 | 509 | 510 | 511 | 512 | 513 | 514 | 515 | 516 | 517 | 518 | 519 | 520 | 521 | 522 | 523 | 524 | 525 |
526 | 527 | 528 | 529 | 530 | 531 | 532 | 533 | 534 | 535 | 536 | 537 | 538 | 539 | 540 | 541 | 542 | 543 | 544 | 545 | 546 | 547 | 548 | 549 | 550 |
551 | 552 | 553 | 554 | 555 | 556 | 557 | 558 | 559 | 560 | 561 | 562 | 563 | 564 | 565 | 566 | 567 | 568 | 569 | 570 | 571 | 572 | 573 | 574 | 575 |
576 | 577 | 578 | 579 | 580 | 581 | 582 | 583 | 584 | 585 | 586 | 587 | 588 | 589 | 590 | 591 | 592 | 593 | 594 | 595 | 596 | 597 | 598 | 599 | 600 |
601 | 602 | 603 | 604 | 605 | 606 | 607 | 608 | 609 | 610 | 611 | 612 | 613 | 614 | 615 | 616 | 617 | 618 | 619 | 620 | 621 | 622 | 623 | 624 | 625 |
626 | 627 | 628 | 629 | 630 | 631 | 632 | 633 | 634 | 635 | 636 | 637 | 638 | 639 | 640 | 641 | 642 | 643 | 644 | 645 | 646 | 647 | 648 | 649 | 650 |
651 | 652 | 653 | 654 | 655 | 656 | 657 | 658 | 659 | 660 | 661 | 662 | 663 | 664 | 665 | 666 | 667 | 668 | 669 | 670 | 671 | 672 | 673 | 674 | 675 |
676 | 677 | 678 | 679 | 680 | 681 | 682 | 683 | 684 | 685 | 686 | 687 | 688 | 689 | 690 | 691 | 692 | 693 | 694 | 695 | 696 | 697 | 698 | 699 | 700 |
701 | 702 | 703 | 704 | 705 | 706 | 707 | 708 | 709 | 710 | 711 |
I scored 13526 points (out of 15700 possible points, top rank was 17 out of ≈60000 in August 2017) at Hackerrank's Project Euler+.
My username at Project Euler is stephanbrumme while it's stbrumme at Hackerrank.
Look at my progress and performance pages to get more details.
Copyright
I hope you enjoy my code and learn something - or give me feedback how I can improve my solutions.
All of my solutions can be used for any purpose and I am in no way liable for any damages caused.
You can even remove my name and claim it's yours. But then you shall burn in hell.
The problems and most of the problems' images were created by Project Euler.
Thanks for all their endless effort !!!
<< problem 501 - Eight Divisors | Tangent Circles - problem 510 >> |