<< problem 267 - Billionaire | Sum of Squares - problem 273 >> |
Problem 268: Counting numbers with at least four distinct prime factors less than 100
(see projecteuler.net/problem=268)
It can be verified that there are 23 positive integers less than 1000 that are divisible by at least four distinct primes less than 100.
Find how many positive integers less than 10^16 are divisible by at least four distinct primes less than 100.
My Algorithm
There are 25 primes less than 100.
I start a loop running from 0 to 2^25. I multiply all the n-th primes if the n-th bit of the counter (named mask
) is set.
If mask
has at least four bits set then the product of those primes needs to be further analyzed:
product
can exceed 10^16 and cause strange results, reject all those large products- there are 10^16 / product numbers divisible by
product
for example 2310 is divisible by 2,3,5,7 and 11 (in fact 2*3*5*7*11 = 2310). But 2310 is also divisible by (2,3,5,7), (2,3,5,11), (2,5,7,11) and (3,5,7,11).
That means (2,3,5,7,11) "overlaps" all combinations of its primes where I take one prime away.
The tetrahedral numbers can compute this "overlap" (see en.wikipedia.org/wiki/Tetrahedral_number) together with
the inclusion-exclusion principle (see en.wikipedia.org/wiki/Inclusionâ€“exclusion_principle).
In plain English: I add the count of numbers divisible by 4 primes, subtract the count divisible by 5 primes, add 6s, subtract 7s, etc.
→ if
numBits
is even, then add else subtract, the fast way is just to look at the lowest bit of the numBits
counter
Interactive test
You can submit your own input to my program and it will be instantly processed at my server:
This is equivalent toecho "4 25 1000" | ./268
Output:
Note: the original problem's input 4 25 10000000000000000
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 <vector>
// all 25 primes less than 100
std::vector<unsigned int> primes = { 2, 3, 5, 7, 11,
13, 17, 19, 23, 29,
31, 37, 41, 43, 47,
53, 59, 61, 67, 71,
73, 79, 83, 89, 97 };
// count numbers with at least "minPrimes" prime factors
// slow version - helping me while debugging the fast code below
unsigned int bruteForce(unsigned int limit, unsigned int minPrimes)
{
unsigned int result = 0;
for (unsigned int i = 1; i <= limit; i++)
{
// count distinct prime factors
unsigned int numPrimeFactors = 0;
for (auto p : primes)
if (i % p == 0)
numPrimeFactors++;
// at least four prime factors ?
if (numPrimeFactors >= minPrimes)
result++;
}
return result;
}
// return the n-th tetrahedral number
unsigned int tetrahedral(unsigned int n)
{
// see https://en.wikipedia.org/wiki/Tetrahedral_number
return n * (n + 1) * (n + 2) / 6;
}
int main()
{
// at least 4 primes, search up to 10^16
unsigned int minPrimes = 4;
unsigned int numPrimes = primes.size();
unsigned long long limit = 10000000000000000ULL;
std::cin >> minPrimes >> numPrimes >> limit;
// find result for small limits
//std::cout << bruteForce(limit, minPrimes) << std::endl;
// precompute tetrahedral numbers
std::vector<unsigned long long> count(numPrimes + 1, 0);
for (unsigned int i = minPrimes; i < count.size(); i++)
count[i] = tetrahedral(i - minPrimes + 1); // same as choose(i - 1, minPrimes - 1)
unsigned long long sum = 0;
unsigned int maxMask = (1 << numPrimes) - 1;
for (unsigned int mask = 0; mask < maxMask; mask++)
{
// multiply all primes
unsigned long long product = 1;
// but don't exceed 10^16
bool tooLarge = false;
// at least four bits must be set (=> pick at least four different primes)
unsigned int numBits = 0;
// look at each prime
for (unsigned int current = 0; current < numPrimes; current++)
{
// is this prime used ?
unsigned int bitpos = 1 << current;
if ((mask & bitpos) == 0)
continue;
// include this prime factor
numBits++;
product *= primes[current];
// avoid too large products (beware of 64 bit overflows ...)
if (product > limit)
{
tooLarge = true;
// speed optimization:
// if the next mask has the same number of bits set (or more), then it will exceed the limit as well
// because the primes are getting larger
// => get closer to the next number with the less bits set
// add the lowest set bit to the current number
// => worst case: same number of bits if the last bit is a single bit
// => best case: several bits next to the lowest bit are set to and all flip to zeros (and a zero becomes a one)
auto withLowestBit = mask & (mask - 1);
auto lowestBit = mask - withLowestBit;
mask += lowestBit;
// hope for the best case ...
while (mask & (lowestBit << 1))
{
lowestBit <<= 1;
mask += lowestBit;
}
mask--; // because the for-loop will add 1 in each iteration
break;
}
}
// skip if too few primes or their product became too large
if (tooLarge || numBits < minPrimes)
continue;
// count numbers divisible by the current primes
auto divisible = limit / product;
// they were part of products with less primes, too, hence compute the "overlap"
divisible *= count[numBits];
// add if even number of primes, else subtract (inclusion-exclusion principle)
if (numBits & 1)
sum -= divisible;
else
sum += divisible;
}
// that's it !
std::cout << sum << std::endl;
return 0;
}
This solution contains 17 empty lines, 31 comments and 2 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 0.9 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 27, 2017 submitted solution
September 27, 2017 added comments
Difficulty
Project Euler ranks this problem at 70% (out of 100%).
Links
projecteuler.net/thread=268 - the best forum on the subject (note: you have to submit the correct solution first)
Code in various languages:
Python github.com/evilmucedin/project-euler/blob/master/euler268/268.py (written by Den Raskovalov)
Python github.com/smacke/project-euler/blob/master/python/268.py (written by Stephen Macke)
C++ github.com/evilmucedin/project-euler/blob/master/euler268/268.cpp (written by Den Raskovalov)
C++ github.com/roosephu/project-euler/blob/master/268.cpp (written by Yuping Luo)
C github.com/LaurentMazare/ProjectEuler/blob/master/e268.c (written by Laurent Mazare)
Java github.com/thrap/project-euler/blob/master/src/Java/Problem268.java (written by Magnus Solheim Thrap)
Perl github.com/shlomif/project-euler/blob/master/project-euler/268/euler-268-v2.pl (written by Shlomi Fish)
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 |
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 267 - Billionaire | Sum of Squares - problem 273 >> |