<< problem 613 - Pythagorean Ant |
Problem 615: The millionth number with at least one million prime factors
(see projecteuler.net/problem=615)
Consider the natural numbers having at least 5 prime factors, which don't have to be distinct.
Sorting these numbers by size gives a list which starts with:
32 = 2 * 2 * 2 * 2 * 2
48 = 2 * 2 * 2 * 2 * 3
64 = 2 * 2 * 2 * 2 * 2 * 2
72 = 2 * 2 * 2 * 3 * 3
80 = 2 * 2 * 2 * 2 * 5
96 = 2 * 2 * 2 * 2 * 2 * 3
...
So, for example, the fifth number with at least 5 prime factors is 80.
Find the millionth number with at least one million prime factors.
Give your answer modulo 123454321.
My Algorithm
Based on the example I figured out pretty soon that most of the "one million prime factors" will be 2s.
The one millionth prime number is 15485863, therefore the smallest number with at least one million prime factors cannot exceed 2^999999 * 15485863 (see tooLarge
).
I assumed that all but the last 32 prime factors will be 2s (in fact I later found empirically that all but the last 27 must be 2s).
A minimum-priority-queue, that means smallest numbers first, is initialized with 2^32 (later changed to 2^27).
Then I pick and remove the smallest number current
one million times and:
- "append a prime": multiply
current
by p where p will be any of the first million primes - but this time I take care that p is no smaller than the largest prime factor ofcurrent
- "replace a prime": replace one of the 2s of
current
by p where p will be any of the first million primes - but this time I take care that p is no smaller than the largest prime factor ofcurrent
(same condition as above) - if the new numbers exceed the limit mentioned above, then I continue with the next number
top()
element will be the smallest number with at least 32 (later reduced to 27) prime factors.Obviously I can't multiply that number by 2^{1000000-32} with C++ native data types.
However, only the residue modulo 123454321 is asked for:
(a * b) mod m = ((a mod m) * (b mod m)) mod m
Multiplying that number 999968 times by 2 while repeatedly applying modulo 123454321 seems to be pretty slow
but GCC does an amazing job at optimizing this loop and it finished in a few milliseconds.
Note
The code was slightly modified after I submitted the correct result: initially, I had an std::set<Number> candidates
which spent most of the time allocating (and freeing) memory.
The std::priority_queue
is about twice as fast when its internal container is an std::vector
.
Most of the one million primes will never be used so I looked for the smallest prime that still produces the correct result:
Starting with 173207, which is the 15770th prime, instead of 15485863 brings down the execution time to about one third of a second (see #define FAST
).
I found that number solely by trial'n'error.
There is probably a smarter way to avoid the duplicates in candidates
but I wasn't willing to spend much time on it.
Interactive test
This feature is not available for the current problem.
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>
#include <queue>
#include <deque>
#include <functional>
#include <set>
// ---------- a few routines from my toolbox ----------
// odd prime numbers are marked as "true" in a bitvector
std::vector<bool> sieve;
// return true, if x is a prime number
bool isPrime(unsigned int x)
{
// handle even numbers
if ((x & 1) == 0)
return x == 2;
// lookup for odd numbers
return sieve[x >> 1];
}
// find all prime numbers from 2 to size
void fillSieve(unsigned int size)
{
// store only odd numbers
const unsigned int half = size >> 1;
// allocate memory
sieve.resize(half, true);
// 1 is not a prime number
sieve[0] = false;
// process all relevant prime factors
for (unsigned int i = 1; 2*i*i < half; i++)
// do we have a prime factor ?
if (sieve[i])
{
// mark all its multiples as false
unsigned int current = 3*i+1;
while (current < half)
{
sieve[current] = false;
current += 2*i+1;
}
}
}
// ---------- problem-specific code ----------
// a tiny wrapper for a number and its largest prime factor
struct Number
{
unsigned long long value; // a number
unsigned int largestFactor; // and its largest prime factor
// for std::greater / min-queue
bool operator>(const Number& other) const
{
return value > other.value;
}
};
int main()
{
unsigned int numCandidates = 1000000;
std::cin >> numCandidates;
const unsigned int Modulo = 123454321;
// generate the first million prime numbers
unsigned int maxPrime = 15485863; // http://www.wolframalpha.com/input/?i=1000000th+prime
#define FAST
#ifdef FAST
if (numCandidates <= 1000000)
maxPrime = 173207; // empirically determined limit
#endif
// compute enough prime numbers
fillSieve(maxPrime + 1);
// copy prime numbers to a consecutive array (just for faster access)
std::vector<unsigned int> primes = { 2 };
primes.reserve(numCandidates);
for (unsigned int i = 3; i <= maxPrime; i += 2)
if (isPrime(i))
primes.push_back(i);
// empirically determined
unsigned int numVariableFactors = 27;
if (numCandidates < 27)
numVariableFactors = numCandidates; // for live-test
unsigned long long seed = 1 << numVariableFactors;
// create one million numbers with exactly 27 prime factors of the form prime*2^26
std::priority_queue<Number, std::vector<Number>, std::greater<Number>> candidates;
candidates.push({ seed, 2 });
// currently 1000000th number (there are a few smaller numbers that need to be discovered now)
auto tooLarge = seed / 2 * primes.back(); // largest number in the candidates container
// detect duplicates
unsigned long long previous = 0;
// look at the candidates in ascending order
for (unsigned int numProcessed = 0; numProcessed < numCandidates; numProcessed++)
{
// pick next number
auto i = candidates.top();
auto current = i.value;
candidates.pop();
// skip duplicates
while (current == previous)
{
i = candidates.top();
current = i.value;
candidates.pop();
}
auto maxFactor = i.largestFactor;
previous = current;
// "append" a prime to the current candidate
for (auto p : primes)
{
auto next = current * p;
if (next >= tooLarge)
break;
if (p >= maxFactor)
candidates.push({ next, p });
}
// replace a 2 by a prime (pretty much the same loop, only "next" is computed slightly different)
for (auto p : primes)
{
auto next = (current / 2) * p;
if (next >= tooLarge)
break;
if (p >= maxFactor && p > 2) // no use in replacing 2 by 2
candidates.push({ next, p });
}
}
// pick the 1000000th number
unsigned int result = candidates.top().value % Modulo;
// multiply by 2^1000000 (except for the 27x 2's "used")
for (unsigned int i = numVariableFactors; i < numCandidates; i++)
result = (result * 2) % Modulo; // I thought that simple loop would be slow (because of "expensive" modulo),
// but GCC converts it to impressively fast code !
std::cout << result << std::endl;
return 0;
}
This solution contains 24 empty lines, 30 comments and 9 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.
Peak memory usage was about 37 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
December 4, 2017 submitted solution
December 4, 2017 added comments
Links
projecteuler.net/thread=615 - the best forum on the subject (note: you have to submit the correct solution first)
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 |
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 |
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 613 - Pythagorean Ant |