<< problem 139 - Pythagorean tiles | Perfect Square Collection - problem 142 >> |
Problem 141: Investigating progressive numbers, n, which are also square
(see projecteuler.net/problem=141)
A positive integer, n, is divided by d and the quotient and remainder are q and r respectively.
In addition d, q, and r are consecutive positive integer terms in a geometric sequence, but not necessarily in that order.
For example, 58 divided by 6 has quotient 9 and remainder 4. It can also be seen that 4, 6, 9 are consecutive terms in a geometric sequence (common ratio 3/2).
We will call such numbers, n, progressive.
Some progressive numbers, such as 9 and 10404 = 102^2, happen to also be perfect squares.
The sum of all progressive perfect squares below one hundred thousand is 124657.
Find the sum of all progressive perfect squares below one trillion (10^12).
My Algorithm
The problem asks for solution of
(1) n = d * q + r
I can assume that d <= q (otherwise I would just swap both)
(2) r < d <= q
The geometric progression is
(3) d / r = q / d
which is equivalent to
(4) d = sqrt{q * r} (remember: d must be an integer)
Now n is
(5) n = sqrt{q * r} * q + r
My simple bruteForce
algorithm finds all solutions for n < 100000 is less than 0.01 seconds. And would take forever to solve n < 10^12.
Nevertheless, that's what bruteForce
does:
- iterate over all r < q < 100000
- compute n using equation (5)
- make sure that n < 100000 and n is a perfect square and q * r is a perfect square
(6) k = d / r = q / d
(7) d = k * r
(8) q = k^2 * r
If q is an integer then q = k^2 * r must be an integer, too.
Since k is rational there are some values x and y
(9) k = dfrac{x}{y}
(10) k^2 = dfrac{x^2}{y^2}
(11) q = dfrac{x^2}{y^2} * r
In equation (11) the variable r must be a multiple of y^2 otherwise q wouldn't be an integer.
That means r is the product of some integer and a perfect square y^2.
My program iterates over all d and all factor and all y such that r = factor * y^2 and r < d. From (3) follows:
(12) q = d^2 * r
(13) n = dfrac{d^3}{r} + r
Knowing potential candidates for d and r I have to check whether frac{d^3}{r} is actually an integer (it's not for most pairs !).
If n < 10^12 and n is a perfect square then another solution was found. There are only 23 solution for n < 10^12.
Be careful: there are multiple combinations of
factor
and y*y
which evaluate to the same remainder
. For example 36 = 9 * 2*2 = 4 * 3*3
.Each remainder must only be used once for each divisor. I saw in my logs that is at most one solution per remainder and divisor.
I don't know why - so maybe my fast exit with the variable
found
is wrong.I'm not very proud that there is a heuristic in the most inner loop which aborts
if (n > 2*limit)
.There is no true science behind that line - I just know it works for 10^5 and 10^12.
Actually I can abort
if (n > limit)
for 10^12 but it misses one solution for 10^5.In the end, this heuristic shortens my execution time to a reasonable level but it could produce wrong answers for other inputs.
Alternative Approaches
After I solved this problem I took a look at other solutions. Most people came up with a much faster and reliable way which I missed completely (their code finishes in less than a second).
Note
It took me a long time find a solution (well, bruteForce
was easy but my final algorithm is completely different).
There are multiple "hacks" in my code - yes, I can solve the two cases 100000 and 10^12 but have no idea how many inputs fail.
So be extra careful when using my live test !
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho 100000 | ./141
Output:
Note: the original problem's input 1000000000000
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 <cmath>
// true if n is a perfect square
template <typename T>
bool isSquare(T n)
{
T root = sqrt(n);
return root * root == n;
}
// iterate over all remainders and quotients, it's okay for limit = 100000 but too slow for 10^12
void bruteForce(unsigned int limit)
{
unsigned int sum = 0;
for (unsigned int remainder = 1; remainder < limit; remainder++)
for (unsigned int quotient = remainder + 2; quotient < limit; quotient++)
{
// d / r = q / d, therefore d^2 = r * q
auto divisorSquared = remainder * quotient;
unsigned int divisor = sqrt(divisorSquared);
// find n
auto n = divisor * quotient + remainder;
// too large ?
if (n > limit)
break;
// not a perfect square ?
if (!isSquare(n))
continue;
// most d^2 are actually NOT a perfect square, therefore all computations above yield garbage
if (!isSquare(divisorSquared))
continue;
// yes, found another combination
//std::cout << n << " = " << divisor << " * " << quotient << " + " << remainder << std::endl;
sum += n;
}
std::cout << sum << std::endl;
}
int main()
{
unsigned long long limit = 1000000000000ULL;
std::cin >> limit;
//bruteForce(limit); return 0;
unsigned long long sum = 0;
// remainder < divisor <= quotient
for (unsigned long long divisor = 1; divisor * divisor <= limit; divisor++)
{
// it seems there is at most one solution per divisor
bool found = false;
// all possible perfect squares ...
for (unsigned long long y = sqrt(divisor); y >= 1 && !found; y--)
// ... with all possible factors ...
for (unsigned long long factor = 1; factor * y * y <= divisor; factor++)
{
// ... generate all possible remainders
auto remainder = factor * y * y;
// n = d^3 / r + r
auto n = divisor * divisor * divisor / remainder + remainder;
// heuristic (MIGHT FAIL !) to skip very large n
if (n > 2*limit)
break;
// ignore if beyond limit
if (n > limit)
continue;
// must be a perfect square (I moved this check in front of the next check because that turned out to be faster)
if (!isSquare(n))
continue;
// there's a high likelihood that d^3 isn't a multiple of r and therefore d^3/r isn't an integer
auto quotient = n / divisor;
// cross-check whether all variables produce a valid result
if (divisor * quotient + remainder != n)
continue;
sum += n;
// it seems there is at most one solution per divisor
found = true;
break;
// note: I kept track of various remainders (to avoid repetitions for the same divisor)
// aborting after a valid solution eliminates the need
}
}
// display result
std::cout << sum << std::endl;
return 0;
}
This solution contains 19 empty lines, 25 comments and 2 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 21.3 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
September 1, 2017 submitted solution
September 1, 2017 added comments
Difficulty
Project Euler ranks this problem at 60% (out of 100%).
Links
projecteuler.net/thread=141 - the best forum on the subject (note: you have to submit the correct solution first)
Code in various languages:
C# www.mathblog.dk/project-euler-141investigating-progressive-numbers-n-which-are-also-square/ (written by Kristian Edlund)
C# github.com/HaochenLiu/My-Project-Euler/blob/master/141.cs (written by Haochen Liu)
Python github.com/hughdbrown/Project-Euler/blob/master/euler-141.py (written by Hugh Brown)
Python github.com/Meng-Gen/ProjectEuler/blob/master/141.py (written by Meng-Gen Tsai)
Python github.com/steve98654/ProjectEuler/blob/master/141.py
C++ github.com/roosephu/project-euler/blob/master/141.cpp (written by Yuping Luo)
C++ github.com/smacke/project-euler/blob/master/cpp/141.cpp (written by Stephen Macke)
C++ github.com/steve98654/ProjectEuler/blob/master/141.cpp
Java github.com/thrap/project-euler/blob/master/src/Java/Problem141.java (written by Magnus Solheim Thrap)
Go github.com/frrad/project-euler/blob/master/golang/Problem141.go (written by Frederick Robinson)
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 |
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 139 - Pythagorean tiles | Perfect Square Collection - problem 142 >> |