<< problem 510 - Tangent Circles | Prime triples and geometric sequences - problem 518 >> |
Problem 516: 5-smooth totients
(see projecteuler.net/problem=516)
5-smooth numbers are numbers whose largest prime factor doesn't exceed 5.
5-smooth numbers are also called Hamming numbers.
Let S(L) be the sum of the numbers n not exceeding L such that Euler's totient function phi(n) is a Hamming number.
S(100)=3728.
Find S(10^12). Give your answer modulo 2^32.
My Algorithm
The three fundamental properties of phi are (see en.wikipedia.org/wiki/Euler's_totient_function):
(1) phi(p) = p - 1 if p is prime
(2) phi(p^k) = p^{k-1} * (p - 1) if p is prime
(3) phi(xy) = phi(x) * phi(y) if x and y are coprime
A 5-smooth Hamming number can be written as 2^a 3^b 5^c.
From (1) follows that if the successor of a Hamming number is prime, then that prime has a totient which is a Hamming number:
(4) phi(2^a 3^b 5^c + 1) = 2^a 3^b 5^c if 2^a 3^b 5^c + 1 is prime
From (2) and (3) follows that the totient of each Hamming number is a Hamming number, too:
(5) phi(2^a 3^b 5^c) = phi(2^a) * phi(3^b) * phi(5^c)
(6) = 2^{a-1} * (2-1) * 3^{b-1} * (3-1) * 5^{c-1} * (5-1)
(7) = 2^{a-1} * 3^{b-1} * 2 * 5^{c-1} * 4
(8) = 2^{a+2} * 3^{b-1} * 5^{c-1}
And from (3) and (4) follows that if the totient of two primes p_1 and p_2 is a Hamming numbers each, then their product is a Hamming number as well:
(9) phi(p_1 p_2) = phi(p_1) * phi(p_2)
(10) = 2^{a_1} 3^{b_1} 5^{c_1} * 2^{a_2} 3^{b_2} 5^{c_2}
(11) = 2^{a_1 + a_2} 3^{b_1 + b_2} 5^{c_1 + c_2}
This applies not only to the product of two primes but also to the product of any number of (distinct) primes with a Hamming totient.
To avoid further confusion: let's assume that all those primes are bigger than 5.
Multiplying each such product by 2^k produces another Hamming number because:
(12) phi(2^k p_1 p_2) = phi(2^k) * phi(p) = 2^{k-1} * phi(p)
Multiplying each such product by 3^k produces another Hamming number because:
(13) phi(3^k p) = phi(3^k) * phi(p) = 3^{k-1} * 2 * phi(p)
Multiplying each such product by 5^k produces another Hamming number because:
(14) phi(5^k p) = phi(5^k) * phi(p) = 5^{k-1} * 2^2 * phi(p)
In general: since 2, 3 and 5 are coprime, multiplying each such product by 2^a 3^b 5^c produces another Hamming number (see (3) ).
The algorithm now looks as follows:
- find all Hamming numbers up to 10^12 and store them in
hamming
(3428 numbers) - find all prime numbers that are the successor of a Hamming number and store them in
primes
(545 numbers) - multiply all values of
primes
in every possible way, but each prime must appear at most once per product (and obviously the product must not exceed 10^12). - multiply each such product by all Hamming numbers (again: not exceeding 10^12)
Note
I wasn't sure whether I could have missed a few numbers whose totient is a Hamming number, too.
My program found the correct value for S(100) so I hoped for the best - and indeed, my result for S(10^12) was correct.
The main loop processes a queue of all prime number products.
Initially I was keen on keeping the queue sorted but soon discovered that it doesn't matter.
Even more, it doesn't matter which numbers you process first: processing the last (= largest) numbers first is faster and saves memory, too,
because the queue doesn't grow too much.
Almost all of my solutions which need a prime test either use a simple prime sieve or my implementation of the Miller-Rabin test.
Since the Miller-Rabin test is more than 100 lines long and not easily understandable, I often hesitate to include it.
This time I decided to give my wheel-based prime test a chance (see en.wikipedia.org/wiki/Wheel_factorization , it's part of my toolbox, too).
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho 100 | ./516
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. Or just jump to my GitHub repository.
#include <iostream>
#include <vector>
#include <algorithm>
// ---------- finally I found a problem for my wheel-based prime test (waiting in my toolbox for ages ...) ----------
// wheel-based prime test
bool isPrime(unsigned long long x)
{
// prime test for 2, 3 and 5 and their multiples
if (x % 2 == 0 || x % 3 == 0 || x % 5 == 0)
return x == 2 || x == 3 || x == 5;
// wheel with size 30 (=2*3*5):
// test against 30m+1, 30m+7, 30m+11, 30m+13, 30m+17, 30m+19, 30m+23, 30m+29
// their deltas/increments are:
const unsigned int Delta[] = { 6, 4, 2, 4, 2, 4, 6, 2 };
// start with 7, which is 30*0+7
unsigned long long i = 7;
// 7 belongs to the second test group
auto pos = 1;
// check numbers up to sqrt(x)
while (i*i <= x)
{
// not prime ?
if (x % i == 0)
return false;
// skip forward to next test divisor
i += Delta[pos];
// next delta/increment
pos = (pos + 1) & 7;
}
// passed all tests, must be a prime number
return x > 1;
}
// ---------- problem-specific code ----------
// store all 3428 Hamming number smaller than 10^12
// and all 545 primes that are a successor of a Hamming number
std::vector<unsigned long long> hamming;
std::vector<unsigned long long> primes;
int main()
{
unsigned long long limit = 1000000000000; // 10^12
std::cin >> limit;
// iterate over all numbers 2^a * 3^b * 5^c
for (unsigned long long two = 1; two <= limit; two *= 2)
for (unsigned long long three = 1; two*three <= limit; three *= 3)
for (unsigned long long five = 1; two*three*five <= limit; five *= 5)
{
auto current = two * three * five;
// including 1 which is not divisible by 2, 3, or 5 but still 5-smooth
hamming.push_back(current);
if (isPrime(current + 1) && current > 5)
primes.push_back(current + 1);
}
// in ascending order ...
std::sort(hamming.begin(), hamming.end());
std::sort(primes.begin(), primes.end());
// store all numbers where phi(x) is 5-smooth
// first=number, second=largest prime factor of that number
typedef std::pair<unsigned long long, unsigned long long> Candidate;
std::vector<Candidate> todo;
// "seed value"
todo.push_back({ 1, 1 });
unsigned long long sum = 0;
while (!todo.empty())
{
// pick next candidate
// note: it's much faster to pick larger candidates first to keep the queue short
auto number = todo.back().first;
auto largestPrime = todo.back().second;
todo.pop_back();
// multiply with all Hamming numbers, each result's totient is a Hamming number, too
for (auto x : hamming)
{
auto next = x * number;
if (next > limit)
break;
sum += next;
}
// multiply with all primes that are larger than the currently largest prime
auto nextPrime = std::upper_bound(primes.begin(), primes.end(), largestPrime);
for (; nextPrime != primes.end(); nextPrime++)
{
auto next = *nextPrime * number;
if (next > limit)
break;
todo.push_back({ next, *nextPrime });
}
}
// same as modulo 2^32
sum &= 0xFFFFFFFF; // alternatively, I could have cast it to a 32-bit integer
std::cout << sum << std::endl;
return 0;
}
This solution contains 20 empty lines, 27 comments and 3 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 0.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
November 15, 2017 submitted solution
November 15, 2017 added comments
Difficulty
Project Euler ranks this problem at 20% (out of 100%).
Links
projecteuler.net/thread=516 - 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/516.py (written by Meng-Gen Tsai)
C++ github.com/evilmucedin/project-euler/blob/master/euler516/516.cpp (written by Den Raskovalov)
C++ github.com/roosephu/project-euler/blob/master/516.cpp (written by Yuping Luo)
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 |
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 510 - Tangent Circles | Prime triples and geometric sequences - problem 518 >> |