<< 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 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.

# 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.

# Interactive test

You can submit your own input to my program and it will be instantly processed at my server:

This is equivalent to`echo "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)*

# Benchmark

The correct solution to the original Project Euler problem was found in 0.05 seconds on a Intel® Core™ i7-2600K CPU @ 3.40GHz.

(compiled for x86_64 / Linux, GCC flags: `-O3 -march=native -fno-exceptions -fno-rtti -std=c++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)

# Heatmap

green problems solve the original Project Euler problem and have a perfect score of 100% at Hackerrank, too.

yellow problems score less than 100% at Hackerrank (but still solve the original problem).

gray problems are already solved but I haven't published my solution yet.

blue problems are solved and there wasn't a Hackerrank version of it at the time I solved it or I didn't care about it because it differed too much.

red problems are solved but exceed the time limit of one minute or the memory limit of 256 MByte.

*Please click on a problem's number to open my solution to that problem:*

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 |

I scored 13,183 points (out of 15300 possible points, top rank was 17 out of ≈60000 in August 2017) at Hackerrank's Project Euler+.

Look at my progress and performance pages to get more details.

My username at Project Euler is

**stephanbrumme**while it's stbrumme at Hackerrank.

# 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 >> |