<< problem 612 - Friend numbers | The millionth number with at least one million ... - problem 615 >> |
Problem 613: Pythagorean Ant
(see projecteuler.net/problem=613)
Dave is doing his homework on the balcony and, preparing a presentation about Pythagorean triangles,
has just cut out a triangle with side lengths 30cm, 40cm and 50cm from some cardboard, when a gust of wind blows the triangle down into the garden.
Another gust blows a small ant straight onto this triangle.
The poor ant is completely disoriented and starts to crawl straight ahead in random direction in order to get back into the grass.
Assuming that all possible positions of the ant within the triangle and all possible directions of moving on are equiprobable,
what is the probability that the ant leaves the triangle along its longest side?
Give your answer rounded to 10 digits after the decimal point.
My Algorithm
The side lengths 30, 40 and 50 sound familiar: yes, it's a triangle with a right angle.
Since the size of the triangle doesn't affect the way the ant can move, I replaced it by a simpler 3,4,5 triangle.
Furthermore, let's assume that the triangles landed in such a way that
- the side a=3 aligned with the x-axis of a virtual coordinate system.
- the side b=4 aligned with the y-axis of a virtual coordinate system.
- the ant is at point P
If I know the coordinates of P then I can compute the probability the the ant exits along the side c:
if the ant chooses any direction within the angle gamma then that probabiltiy is gamma / 2 pi
The angles alpha and beta can be derived from the coordinates P(x,y):
(1) tan(alpha) = dfrac{x}{b - y}
(2) tan(beta) = dfrac{y}{a - x}
Solving for alpha and beta:
(3) alpha = arctan(dfrac{x}{b - y})
(4) beta = arctan(dfrac{y}{a - x})
That's exactly what the functions
getAlpha
and getBeta
do.The third angle gamma is based solely on alpha and beta:
- the sum of all angles around the point P is always 2 pi (=360^{\circ})
- all angles in any triangle sum to pi
- a right angle equals frac{pi}{2}
- so the angle "below" P is pi - frac{pi}{2} - alpha = frac{pi}{2} - alpha
- so the angle "on the left side" of P is pi - frac{pi}{2} - beta = frac{pi}{2} - beta
(5) 2 pi = gamma + (frac{pi}{2} - alpha) + (frac{pi}{2} - beta) + frac{pi}{2}
(6) gamma = 2 pi - (frac{pi}{2} - alpha) - (frac{pi}{2} - beta) - frac{pi}{2}
(7) gamma = frac{pi}{2} + alpha + beta
And the probability is (as mentioned before):
(8) f(x,y) = dfrac{gamma}{2 pi}
(9) f(x,y) = dfrac{1}{2 pi} * (dfrac{pi}{2} + arctan(dfrac{x}{b - y}) + arctan(dfrac{y}{a - x}))
The average value f(x,y) of all points (x,y) on the triangle's surface will be the answer to the problem.
However, I couldn't find a fast and reliable way to do it in C++ (without external libraries).
The code you see below is good for the first four digits - but I need to provide ten accurate digits.
That's when I started playing around with Wolfram Alpha:
- the integral \int{\int{f(x,y)}dx}dy divided by the triangle's area dfrac{3 * 4}{2} = 6 is the exact answer
- and Wolfram Alpha simplified that to 0.5 + log((2^16 3^4 sqrt{3/5})/5^12) / 12 pi - I would have never found that formula !
- evaluating the equation gives the correct result, even in C++ ...
Note
I always strive to put all tools needed to solve a problem into a single C++ file.
This time my C++ code is basically useless because an external tool (Wolfram Alpha) condensed the whole solution to a single equation.
Therefore I don't like this problem at all - I found the right approach to the problem and just couldn't take the last step because
it required wizardry that is wildly outside my capabilities.
Reading the forums it appears that 95% of all solvers walked the same path and relied on Mathematica or similar tools.
Interactive test
This feature is not available for the current problem.
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 <cmath>
// triangle (same shape as 30,40,50 !)
const double a = 3; // shortest leg
const double b = 4; // longest leg
const double c = 5; // hypotenuse
// |\.
// | \.
// | \.
// b c
// | \.
// | \.
// +--a-\.
const double PI = 3.1415926535897932384626433832795028;
// compute angle in upper left corner
double getAlpha(double x, double y)
{
// tan(alpha) = opposite / adjacent = X / Y
// alpha = arctan(X / Y)
// where X = x and Y = b - y
return atan(x / (b - y));
}
// compute angle in lower right corner
double getBeta(double x, double y)
{
// similar to getAlpha:
// tan(beta) = opposite / adjacent = Y / X
// beta = arctan(Y / X)
// where X = a - x and Y = y
return atan(y / (a - x));
}
// return probability of leaving the triangle
double getProbability(double x, double y)
{
// the of all angles in a triangle is pi:
// pi = otherAlpha + otherBeta + gamma
// but the big triangle must fulfil the same condition, so
// pi = alpha + otherAlpha + beta + otherBeta + pi/2
// therefore the difference is
// 0 = gamma - alpha - beta - pi/2
// gamme = alpha + beta + pi/2
auto gamma = PI/2 + getAlpha(x, y) + getBeta(x, y);
// each direction is equally likely: "exit probability" = gamma / full circle
// and a full circle has 360 degrees => 2pi
return gamma / (2 * PI);
}
int main()
{
// display 10 digits
std::cout << std::fixed << std::setprecision(10);
double epsilon = 0.01;
//std::cin >> epsilon;
// add all percentages
double sum = 0;
unsigned long long numPoints = 0;
for (double x = epsilon; x < a; x += epsilon)
{
auto height = b - x * b/a;
for (double y = epsilon; y < height; y += epsilon)
{
auto current = getProbability(x, y);
sum += current;
numPoints++;
//if (numPoints % 10000000 == 0)
//std::cout << numPoints << std::endl;
}
}
// compute average percentage
auto result = sum / numPoints;
// ==> uncomment the following line to see my approximate result
//std::cout << result << std::endl;
// I reduced the triangle's dimensions from 30,40,50 to 3,4,5 (doesn't affect the probability)
// and Wolfram Alpha double-integrated to
// 0.5 + log((2^16 3^4 sqrt(3/5))/5^12) / 12pi
// (actually Wolfram Alpha's output contained large numbers: 5308416 = 2^16 3^4 and 244140625 = 5^12)
// I have no clue how Wolfram Alpha found that strange equation
// but the resulting value matched my approximation and it's indeed correct
result = 0.5 + log((pow(2, 16) * pow(3, 4) * sqrt(3.0 / 5)) / pow(5, 12)) / (12*PI);
// and here comes the accurate result
std::cout << result << std::endl;
return 0;
}
This solution contains 12 empty lines, 42 comments and 3 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in less than 0.01 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
November 21, 2017 submitted solution
November 21, 2017 added comments
Links
projecteuler.net/thread=613 - the best forum on the subject (note: you have to submit the correct solution first)
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 |
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 612 - Friend numbers | The millionth number with at least one million ... - problem 615 >> |