<< problem 461 - Almost Pi | Maximum number of divisors - problem 485 >> |
Problem 473: Phigital number base
(see projecteuler.net/problem=473)
Let phi be the golden ratio: phi = 1 + 5 sqrt{2}.
Remarkably it is possible to write every positive integer as a sum of powers of phi
even if we require that every power of phi is used at most once in this sum.
Even then this representation is not unique.
We can make it unique by requiring that no powers with consecutive exponents are used and that the representation is finite.
E.g: 2 = phi + phi^-2 and 3 = phi^2 + phi^-2
To represent this sum of powers of phi we use a string of 0's and 1's with a point to indicate where the negative exponents start.
We call this the representation in the phigital numberbase.
So 1=1_{phi}, 2=10.01_{phi}, 3=100.01_{phi} and 14=100100.001001_{phi}
The strings representing 1, 2 and 14 in the phigital number base are palindromic, while the string representing 3 is not.
(the phigital point is not the middle character).
The sum of the positive integers not exceeding 1000 whose phigital representation is palindromic is 4345.
Find the sum of the positive integers not exceeding 10^10 whose phigital representation is palindromic.
My Algorithm
According to the problem statement, two 1s must be separated by at least one zero.
When I looked at the found numbers, there were always at least TWO zeros between two 1s.
In search()
I generate only one half of the palindrome: the right-hand side. Obviously the left-hand side must be the mirrored version.
If current
is the sum of a few phi^n then I add phi^{n+2}, then phi^{n+3}, phi^{n+4}, ... until the sum gets too large.
Unfortunately, I can't work around precision issues and need GCC's __float128
extension.
Computing phi^n on-the-fly caused even more precision problems - so I precomputed more accurate values with Wolfram Alpha (see large table in phipow()
.
Alternative Approaches
I didn't realize that the golden ratio can be found in Fibonacci numbers, too (see en.wikipedia.org/wiki/Fibonacci_number):
F_n = dfrac{phi^n - (-phi)^{-n}}{2 phi - 1}
This way all my precision problem can be avoided - and even more important, a solution can be found in a few milliseconds.
Note
Working with GCC's __float128
extension requires the std=gnu++11
compiler switch.
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho 1000 | ./473
Output:
Note: the original problem's input 10000000000
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++. You can download it, too. Or just jump to my GitHub repository.
#include <iostream>
#include <vector>
#include <cmath>
typedef __float128 Number; // G++ only !
// acceptable error
const Number Epsilon = 1e-15Q;
// ln(10^10)/ln((1+sqrt(5))/2) is about 47.85, so at most 48 exponents
unsigned int maxExponent = 48;
// return phi^exponent
Number phipow(int exponent)
{
// golden ratio is (1 + sqrt(5)) / 2
//const Number Phi = 1.618033988749894848204586834365638117720309179805762862135Q;
// precompute values Phi^0 ... Phi^48
static const Number precomputed[] =
{
1,
1.618033988749894848204586834365638117720309179805762862135Q,
2.618033988749894848204586834365638117720309179805762862135Q,
4.236067977499789696409173668731276235440618359611525724270Q,
6.854101966249684544613760503096914353160927539417288586406Q,
11.09016994374947424102293417182819058860154589902881431067Q,
17.94427190999915878563669467492510494176247343844610289708Q,
29.03444185374863302665962884675329553036401933747491720776Q,
46.97871376374779181229632352167840047212649277592102010484Q,
76.01315561749642483895595236843169600249051211339593731260Q,
122.9918693812442166512522758901100964746170048893169574174Q,
199.0050249987406414902082282585417924771075170027128947300Q,
321.9968943799848581414605041486518889517245218920298521475Q,
521.0019193787254996316687324071936814288320388947427468775Q,
842.9988137587103577731292365558455703805565607867725990250Q,
1364.000733137435857404797968963039251809388599681515345902Q,
2206.999546896146215177927205518884822189945160468287944927Q,
3571.000280033582072582725174481924073999333760149803290830Q,
5777.999826929728287760652380000808896189278920618091235757Q,
9349.000106963310360343377554482732970188612680767894526588Q,
15126.99993389303864810402993448354186637789160138598576234Q,
24476.00004085634900844740748896627483656650428215388028893Q,
39602.99997474938765655143742344981670294439588353986605128Q,
64079.00001560573666499884491241609153951090016569374634021Q,
103681.9999903551243215502823358659082424552960492336123914Q,
167761.0000059608609865491272482819997819661962149273587317Q,
271442.9999963159853080994095841479080244214922641609711232Q,
439204.0000022768462946485368324299078063876884790883298549Q,
710646.9999985928316027479464165778158308091807432493009781Q,
1149851.00000086967789739648324900772363719686922233763Q,
1860497.99999946250950014442966558553946800604996558693Q,
3010349.00000033218739754091291459326310520291918792456Q,
4870846.99999979469689768534258017880257320896915351149Q,
7881196.00000012688429522625549477206567841188834143605Q,
12752042.9999999215811929115980749508682516208574949475Q,
20633239.0000000484654881378535697229339300327458363836Q,
33385281.9999999700466810494516446738021816536033313311Q,
54018521.0000000185121691873052143967361116863491677147Q,
87403802.9999999885588502367568590705382933399524990459Q,
141422324.000000007071019424062073467274405026301666760Q,
228826126.999999995629869660818932537812698366254165806Q,
370248451.000000002700889084881006005087103392555832567Q,
599074577.999999998330758745699938542899801758809998373Q,
969323029.000000001031647830580944547986905151365830941Q,
1568397606.99999999936240657628088309088670691017582931Q,
2537720636.00000000039405440686182763887361206154166025Q,
4106118242.99999999975646098314271072976031897171748957Q,
6643838879.00000000015051539000453836863393103325914982Q,
10749957121.999999999906976373147249098394250004976639Q // already too big: > 10^10
};
// lookup cached values
if (exponent >= 0)
return precomputed[+exponent];
else
return 1.0Q / precomputed[-exponent];
}
// return phi^pos + phi^-(pos+1)
Number phipowBoth(unsigned int pos)
{
// precompute phi^i + phi^-(i+1)
// => those are the pairs in the palindromic representation
static std::vector<Number> cache;
if (cache.empty())
{
cache.push_back(phipow(0));
for (int i = 1; i < 49; i++)
cache.push_back(phipow(i) + phipow(-(i+1)));
}
return cache[pos];
}
// generate all strings where there are at least two zeros between two ones
unsigned long long search(Number limit, unsigned int exponent = 0, Number current = 0)
{
// "out of bounds" ?
if (current > limit)
return 0;
unsigned long long result = 0;
// an integer ?
auto rounded = round(current); // avoid rounding issues, add a small "buffer"
auto diff = current - rounded;
if (diff > -Epsilon && diff < +Epsilon)
result += rounded;
// it turns out at least two zeros between a pair of ones
// go deeper by appending "001", "0001", "00001", etc.
// => at least two zeros and a single one
// except next to the dot: a zero on one side generates a zero on the other side, too (palindrome)
if (exponent == 0)
exponent++;
else
exponent += 2;
// until phi^48 > 10^10
while (exponent <= maxExponent)
{
// add phi^exponent (and its mirrored "twin" phi^-(exponent+1)
result += search(limit, exponent+1, current + phipowBoth(exponent));
exponent++;
}
return result;
}
int main()
{
// 10^10
unsigned long long limit = 10000000000ULL;
std::cin >> limit;
// find maximum length of one side of the palindrome
maxExponent = 0;
while (phipow(maxExponent) <= limit)
maxExponent++;
// only one phigital number can have a 1 next to the dot => that's the degenerated string "1"
// all other strings must have zeros next to the dot bceause "no powers with consecutive exponents"
std::cout << 1 + search(limit) << std::endl;
return 0;
}
This solution contains 16 empty lines, 23 comments and 3 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 19.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 30, 2017 submitted solution
September 30, 2017 added comments
Difficulty
Project Euler ranks this problem at 35% (out of 100%).
Links
projecteuler.net/thread=473 - 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/473.py (written by Meng-Gen Tsai)
C++ github.com/roosephu/project-euler/blob/master/473.cpp (written by Yuping Luo)
Go github.com/frrad/project-euler/blob/master/golang/Problem473.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 | 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 | 712 | 713 | 714 | 715 | 716 | 717 | 718 | 719 | 720 | 721 | 722 | 723 | 724 | 725 |
726 | 727 | 728 | 729 | 730 | 731 | 732 | 733 | 734 | 735 | 736 | 737 | 738 | 739 | 740 | 741 | 742 | 743 | 744 | 745 | 746 | 747 | 748 | 749 | 750 |
751 | 752 | 753 | 754 | 755 | 756 | 757 | 758 | 759 | 760 | 761 | 762 | 763 | 764 | 765 | 766 | 767 | 768 | 769 | 770 | 771 | 772 | 773 | 774 | 775 |
776 | 777 | 778 | 779 | 780 | 781 | 782 | 783 | 784 | 785 | 786 | 787 | 788 | 789 | 790 | 791 | 792 | 793 | 794 | 795 | 796 | 797 | 798 | 799 | 800 |
801 | 802 | 803 | 804 | 805 | 806 | 807 | 808 | 809 | 810 | 811 | 812 | 813 | 814 | 815 | 816 | 817 | 818 | 819 |
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 461 - Almost Pi | Maximum number of divisors - problem 485 >> |