<< problem 500 - Problem 500!!! | Square on the Inside - problem 504 >> |
Problem 501: Eight Divisors
(see projecteuler.net/problem=501)
The eight divisors of 24 are 1, 2, 3, 4, 6, 8, 12 and 24.
The ten numbers not exceeding 100 having exactly eight divisors are 24, 30, 40, 42, 54, 56, 66, 70, 78 and 88.
Let f(n) be the count of numbers not exceeding n with exactly eight divisors.
You are given f(100) = 10, f(1000) = 180 and f(10^6) = 224427.
Find f(10^12).
My Algorithm
When a number n is split into its prime factors p_1, p_2, ... and each prime factors is found q_1, q_2, ... times then
n = p_1^{q_1} * p_2^{q_2} * ...
The number of divisors is (see en.wikipedia.org/wiki/Divisor_function):
d = (q_1 + 1) * (q_2 + 1) * ...
For d = 8 there are only three possible combinations:
8 = 2 * 2 * 2 = 4 * 2 = 8 * 1
The first case means that n has three different prime factors p_1, p_2 and p_3 where:
q_1 = q_2 = q_3 = 1 when n = p_1 * p_2 * p_3
To fulfil the second case, n needs exactly two different prime factor p_1 and p_2 where:
q_1 = 3 and q_2 = 1 when n = p_1^3 * p_2
And the third case applies n has only one prime factor p_1:
q_1 = 7 when n = p_1^7
24 has eight prime factors because of the second case: 24 = 2^3 * 3.
Case 1 needs the most time: two nested loops iterate over all prime numbers p_1 and p_2 where p_1 < p_2.
Adding a third nested loop for p_3 (with p_2 < p_3) would work but take forever.
That's where the countPrimes(n)
function comes into play: it returns the number of prime up to n
.
To find the number of primes p_3 for p_2 < p_3 <= n: \#(p_3) = countPrimes(n / p_1 p_2) - countPrimes(p_2)
The second case is built on the same idea but needs only a single loops.
There are just a few numbers for the third case and a simple loop without countPrimes
suffices.
I didn't write the countPrimes
function - I found almost working code on am-just-a-nobody.blogspot.de/2015/11/c-code-for-primepi-function.html
and improved it to match my needs. See my code comments for further details.
In other programming languages it is often called PrimePi
.
Note
You can still find my initial brute-force code (see hasEightDivisors
). It can easily solve for n = 10^6 but is much too slow for n = 10^12.
There are probably faster algorithm to count a number of primes. My program finishes in just under a minute but needs only 15 MByte RAM.
This Wikipedia page lists a few alternatives: en.wikipedia.org/wiki/Prime-counting_function
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho 1000000 | ./501
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>
#include <cmath>
// ---------- brute-force code for small numbers, not used anymore ----------
// simple trial-division algorithm, too slow for large numbers
bool hasEightDivisors(unsigned long long n)
{
unsigned int count = 2; // always at least two divisors: 1 and n
for (unsigned long long i = 2; i*i <= n; i++)
if (n % i == 0) // found a divisor ?
{
if (i*i == n)
count++; // perfect square
else
count += 2; // i is a divisor and n/i as well
// abort early
if (count > 8)
return false;
}
// exactly 8 divisors ?
return count == 8;
}
// ---------- and now the currently used code ----------
// return the number of primes up to n
// similar to PrimePi in popular math software
std::vector<unsigned int> primes; // will contain all primes up to sqrt(n)
unsigned long long countPrimes(unsigned long long n)
{
// based on http://am-just-a-nobody.blogspot.de/2015/11/c-code-for-primepi-function.html
// algorithm http://am-just-a-nobody.blogspot.de/2015/11/algorithm-for-summing-all-primes-less.html
// I don't have a thorough understanding of the code ...
// but added these features:
// - allocate only as much memory as actually needed
// - fixed out-of-bounds errors
// - remove the MOD and return the actual result
// - made variables as local as possible
// - find quick result for small n (when primes[] contains enough numbers)
// if primes[] contains enough elements then run a fast binary search
if (!primes.empty() && primes.back() > n)
{
// find smallest number larger than n
auto i = std::upper_bound(primes.begin(), primes.end(), n);
return std::distance(primes.begin(), i);
}
auto v = (unsigned int)sqrt(n);
// about sqrt(n) * 12 bytes, for n = 10^12 => 12 MByte plus primes[]
std::vector<unsigned long long> higher(v+2, 0);
std::vector<unsigned int> lower (v+2, 0);
std::vector<bool> used (v+2, false);
// assume all numbers are prime numbers
unsigned long long result = n - 1;
// the remaining lines subtract composites until result contains the number of primes
// set up lower and upper bound
for (unsigned int p = 2; p <= v; p++)
{
lower [p] = p - 1;
higher[p] = n / p - 1;
}
for (unsigned int p = 2; p <= v; p++)
{
// composite ?
if (lower[p] == lower[p - 1])
continue;
// store prime numbers (if not already existing)
if (primes.empty() || p > primes.back())
primes.push_back(p);
auto temp = lower[p - 1];
// remove more composites
result -= higher[p] - temp;
auto pSquare = (unsigned long long)p * p;
auto end = std::min<unsigned long long>(v, n / pSquare);
// alternate between 1 and 2
auto j = 1 + (p & 1);
// adjust upper bound
for (auto i = p + j; i <= end + 1; i += j)
{
if (used[i])
continue;
auto d = i * p;
if (d <= v)
higher[i] -= higher[d] - temp;
else
higher[i] -= lower[n / d] - temp;
}
// adjust lower bound
for (auto i = v; i >= pSquare; i--)
lower[i] -= lower[i / p] - temp;
// cross off multiples
for (auto i = pSquare; i <= end; i += p*j)
used[i] = true;
}
return result;
}
// count all primes with exactly 8 divisors
unsigned long long fast(unsigned long long n)
{
// note: as a side effect, countPrimes generates all requires prime numbers
// which I need in the upcomig loops
// I don't actually care about the result of countPrimes(n) right now
countPrimes(n); // a dedicated sieve might be faster, though
// a * b * c where a < b < c
unsigned long long countABC = 0;
for (size_t indexA = 0; indexA < primes.size(); indexA++)
{
unsigned long long a = primes[indexA];
if (a * a * a > n)
break;
for (size_t indexB = indexA + 1; indexB < primes.size(); indexB++)
{
auto b = primes[indexB];
// min(c) = next prime after b
// max(c) = last prime before n / (a*b)
auto maxC = n / (a * b);
if (maxC <= b)
break;
// count all primes between min(c) and max(c)
auto high = countPrimes(maxC);
auto low = indexB + 1; // same as countPrimes(b);
countABC += high - low;
}
}
// a^3 * b
unsigned long long countA3B = 0;
for (auto a : primes)
{
auto maxB = n / ((unsigned long long)a * a * a);
if (maxB <= 1)
break;
// b can be any prime
auto numB = countPrimes(maxB);
// but b must differ from a
if (maxB >= a)
numB--;
countA3B += numB;
}
// a^7
unsigned long long countA7 = 0;
for (auto a : primes)
{
if ((unsigned long long) a * a * a * a * a * a * a > n)
break;
countA7++;
}
return countABC + countA3B + countA7;
}
int main()
{
auto limit = 1000000000000ULL;
std::cin >> limit;
// find all number with naive approach (= extremely slow)
//unsigned long long count = 0;
//for (unsigned long long i = 1; i <= limit; i++)
// if (hasEightDivisors(i))
// count++;
//std::cout << count << std::endl;
std::cout << fast(limit) << std::endl;
return 0;
}
This solution contains 31 empty lines, 47 comments and 4 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 56.3 seconds on an Intel® Core™ i7-2600K CPU @ 3.40GHz.
Peak memory usage was about 15 MByte.
(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
August 17, 2017 submitted solution
August 17, 2017 added comments
Difficulty
Project Euler ranks this problem at 40% (out of 100%).
Links
projecteuler.net/thread=501 - the best forum on the subject (note: you have to submit the correct solution first)
Code in various languages:
C++ github.com/evilmucedin/project-euler/blob/master/euler501/501.cpp (written by Den Raskovalov)
C++ github.com/ngthanhtrung23/ProjectEuler/blob/master/P5XX/501.cpp (written by Trung Nguyen)
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 |
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 500 - Problem 500!!! | Square on the Inside - problem 504 >> |