<< problem 211 - Divisor Square Sum | Totient Chains - problem 214 >> |
Problem 213: Flea Circus
(see projecteuler.net/problem=213)
A 30x30 grid of squares contains 900 fleas, initially one flea per square.
When a bell is rung, each flea jumps to an adjacent square at random (usually 4 possibilities, except for fleas on the edge of the grid or at the corners).
What is the expected number of unoccupied squares after 50 rings of the bell? Give your answer rounded to six decimal places.
My Algorithm
Assume there is a 3x2 with two fleas A and B:A..
.B.
The probability land_{A,1}(x,y) of flea A to land on a certain square after one jump is:
00.50
0.500
Note: Flea A had 2 options when it started at square (1,2), therefore directions(1,2) = 100% / 2 = 50% = 0.5.
The probability land_{B,1}(x,y) of flea B to land on a certain square after one jump is:
00.3330
0.33300.333
Note: Flea B had 3 options when it started at square (2,1), therefore directions(2,1) = 100% / 3 = 33.3% = 0.333.
To find the probability of flea A to be NOT on a certain square I have to compute empty_{A,1}(x,y) = 1 - land_{A,1}(x,y):
10.51
0.511
And the same for flea B: empty_{B,1}(x,y) = 1 - land_{B,1}(x,y):
10.6661
0.66610.666
A square is empty is neither flea A nor flea B landed on it: empty_1(x,y) = empty_{A,1}(x,y) * empty_{B,1}(x,y):
10.3331
0.33310.666
So the expected total number of empty squares after one jump is sum{empty_1} = 1 + 0.333 + 1 + 0.333 + 1 + 0.666 = 4.333.
There were exactly four empty squares before the first jump but there are on average 4.333 empty squares after one jump (→ plus 0.333)
because it's possible that both fleas land on the same square, leaving five square empty instead of four.
In the second iteration, the probability land_{A,2}(x,y) of flea A depends on the probability of being on the square in the previous iteration, too.
The probabilities land_{A,2}(x,y) of flea A after two jumps are:
0.41600.166
00.416*0
Let's examine the square land_{A,2}(2,1) which is denoted by a star:
land_{A,2}(2,1) = land_{A,1}(2,2) * directions(2,2) + land_{A,1}(1,1) * directions(1,1) = 0.5 * 0.333 + 0.5 * 0.5 = 0.416.
Implementation
Initially all squares are empty, therefore
empty[x][y] = 1
.Then my program tracks each a single flea (such as flea A or B) over 50 rounds.
It has a grid
current
which is intially zero and contains a single 1 at the start position of the flea.After each jump, the probabilities of each square the flea can jump to are updated and stored in
next
.If all squares are processed,
next
is copied to current
.After 50 rounds, the "emptiness" of each square is
empty[x][y] *= 1 - current[x][y]
.This process is repeated for all 30x30 fleas after which the sum of all squares in
empty
is the result.
Alternative Approaches
Even though the algorithm is pretty simple, I struggled with a few details and couldn't find the correct answer for a while.
Therefore I wrote a Monte-Carlo simulation "to get a feeling" for the result.
Obviously, the Monte-Carlo simulation wasn't able to compute an exact result, even after millions of repetitions.
But it helped me to spot my errors and fix them.
Note
The grid offers several opportunies to exploit symmetries.
I mirror along the x- and/or y-axis if possible (see #define FAST_SYMMETRY
).
In the end I get a 3x speedup.
You could mirror along the diagonals, too, but then the code might get messy.
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho "3 2 1" | ./213
Output:
Note: the original problem's input 30 30 50
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.
#include <iostream>
#include <iomanip>
#include <vector>
// size of the flea circus
unsigned int width = 30;
unsigned int height = 30;
// number of jumps for each flea
unsigned int rounds = 50;
// a fixed 2D array
typedef std::vector<std::vector<double>> Grid;
// it looks weird, but the inner array describes the second index
// create a 2D (width x height) and set all elements to initialValue
Grid makeGrid(double initialValue)
{
Grid result(width);
for (auto& row : result)
row.resize(height, initialValue);
return result;
}
int main()
{
// read parameters
std::cin >> width >> height >> rounds;
// probability that squares are empty
// ==> no fleas yet, therefore "emptiness" probability is 1 (=100%) for each square
Grid empty = makeGrid(1);
// use symmetry to speed up the process
bool mirrorX = false;
bool mirrorY = false;
unsigned int maxX = width;
unsigned int maxY = height;
#define FAST_SYMMETRY
#ifdef FAST_SYMMETRY
if (width % 2 == 0)
{
mirrorX = true;
maxX = width / 2;
}
if (height % 2 == 0)
{
mirrorY = true;
maxY = height / 2;
}
#endif
// analyze each flea
for (unsigned int startX = 0; startX < maxX; startX++)
for (unsigned int startY = 0; startY < maxY; startY++)
{
// drop a single flea at (x,y):
// probability of seeing the flea is zero for each square ...
Grid current = makeGrid(0);
// ... except for the square where it starts
current[startX][startY] = 1;
// trace probability of landing at each square after a few rounds
for (unsigned int round = 0; round < rounds; round++)
{
// track probability of seeing the flea after one more jump
Grid next = makeGrid(0);
// look at each square from the previous round
for (unsigned int x = 0; x < width; x++)
for (unsigned int y = 0; y < height; y++)
{
// flea couldn't be there ? ==> skip
if (current[x][y] == 0)
continue;
// count number of possible destination squares
unsigned int directions = 4;
if (x == 0 || x == width - 1)
directions--;
if (y == 0 || y == height - 1)
directions--;
// probability of landing on a square is
// "probability of starting at current square" / "number of options"
double probability = current[x][y] / directions;
// add probability to each allowed destination square
if (x > 0)
next[x - 1][y] += probability;
if (x < width - 1)
next[x + 1][y] += probability;
if (y > 0)
next[x][y - 1] += probability;
if (y < height - 1)
next[x][y + 1] += probability;
}
current = std::move(next);
}
// 1. current[x][y] contains the probability that our flea is on square (x,y) after xx rounds
// 2. the probability of being empty is the opposite, it's: 1 - current[x][y]
// 3. and the probability that no flea finishes on that square (no matter where it started)
// is the product of the "emptiness" probabilities of all fleas
for (unsigned int x = 0; x < width; x++)
for (unsigned int y = 0; y < height; y++)
{
empty[x][y] *= 1 - current[x][y];
// use symmetries
if (mirrorX)
empty[x][y] *= 1 - current[width - 1 - x][y];
if (mirrorY)
empty[x][y] *= 1 - current[x][height - 1 - y];
if (mirrorX && mirrorY)
empty[x][y] *= 1 - current[width - 1 - x][height - 1 - y];
}
}
// sum of all probabilities
double result = 0;
for (unsigned int x = 0; x < width; x++)
for (unsigned int y = 0; y < height; y++)
result += empty[x][y];
// display result
std::cout << std::fixed << std::setprecision(6) << result << std::endl;
return 0;
}
This solution contains 17 empty lines, 28 comments and 6 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 0.05 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
July 5, 2017 submitted solution
July 5, 2017 added comments
Difficulty
Project Euler ranks this problem at 60% (out of 100%).
Links
projecteuler.net/thread=213 - the best forum on the subject (note: you have to submit the correct solution first)
Code in various languages:
Python github.com/hughdbrown/Project-Euler/blob/master/euler-213.py (written by Hugh Brown)
Python github.com/Meng-Gen/ProjectEuler/blob/master/213.py (written by Meng-Gen Tsai)
Python github.com/smacke/project-euler/blob/master/python/213.py (written by Stephen Macke)
C++ github.com/roosephu/project-euler/blob/master/213.cpp (written by Yuping Luo)
Java github.com/HaochenLiu/My-Project-Euler/blob/master/213.java (written by Haochen Liu)
Java github.com/thrap/project-euler/blob/master/src/Java/Problem213.java (written by Magnus Solheim Thrap)
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 |
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 |
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 211 - Divisor Square Sum | Totient Chains - problem 214 >> |