<< problem 329 - Prime Frog | Maximix Arrangements - problem 336 >> |
Problem 333: Special partitions
(see projecteuler.net/problem=333)
All positive integers can be partitioned in such a way that each and every term of the partition can be expressed as 2^i * 3^j, where i,j >= 0.
Let's consider only those such partitions where none of the terms can divide any of the other terms.
For example, the partition of 17 = 2 + 6 + 9 = 2^1 * 3^0 + 2^1 * 3^1 + 2^0 * 3^2 would not be valid since 2 can divide 6.
Neither would the partition 17 = 16 + 1 = 2^4 * 3^0 + 2^0 * 3^0 since 1 can divide 16.
The only valid partition of 17 would be 8 + 9 = 2^3 * 3^0 + 2^0 * 3^2.
Many integers have more than one valid partition, the first being 11 having the following two partitions.
11 = 2 + 9 = 2^1 * 3^0 + 2^0 * 3^2
11 = 8 + 3 = 2^3 * 3^0 + 2^0 * 3^1
Let's define P(n) as the number of valid partitions of n. For example, P(11) = 2.
Let's consider only the prime integers q which would have a single valid partition such as P(17).
The sum of the primes q < 100 such that P(q)=1 equals 233.
Find the sum of the primes q < 1000000 such that P(q)=1.
My Algorithm
My bruteForce()
attempt can verify P(q < 100) = 1: it creates a set of all powers 2^i 3^j < 100 and then analyzes each subset.
Unfortuately there are 143 such powers for q < 10^6 so that bruteForce()
would have to check 2^143 subsets - that's impossible.
So I "restarted" my thinking process: if I have a term 2^a 3^b and then each other term 2^c 3^d must be neither a multiple nor a divisor.
That means if a > c then b < d because otherwise 2^a 3^b would be a multiple of 2^c 3^d.
And if a < c then b > d because otherwise 2^a 3^b would be a divisor of 2^c 3^d.
I am free to "sort" all elements of a partition and decided to apply the following rule:
the first element has the smallest 2^x and the largest 3^y of all terms of the partition.
That's exactly how the example is sorted (coincidence ?):
17 = 2^3 * 3^0 + 2^0 * 3^2
The first step is iterate over each possible 2^a * 3^b so I can identify all numbers that can be reached by partitions consisting of only one term.
Then I go through all numbers in ascending order:
only if there was a way to reach a number with (at least) one term then I could add one more term - that must obey the rule if a < c then b > d
My actual implementation may look a bit confusing because I tried hard to trim memory consumption:
- a unique ID is assigned to each of 143 terms 2^a * 3^b < 10^6 (see container
ids
) - create 143 bit vectors
reachable
with 10^6 elements: an element is true if there is at least one partition and the "last" term carries that ID - create 143 bit vectors
multiple
with 10^6 elements: an element is true if there is more than one partition and the "last" term carries that ID
reachable
and multiple
are filled step-by-step as I iterate over all 10^6 numbers.- ignore each ID of the current number
i
wherereachable[id][i] == false
- else add one more term: iterate over each possible successor where a < c and b > d
- set
reachable[i + newTerm]
accordingly, updatemultiple[i + newTerm]
as well - if the current number
i
has at most onereachable[id][i]
that is true (and even this must bemultiple[id][i] == false
) then I found a unique partition
Alternative Approaches
It seems you could write a Dynamic Programming algorithm - I tried but failed.
Note
My final program was heavily affected by the desire to keep memory consumption low (it still needs about 38 MByte).
The bit vectors reachable
and multiple
could be merged into a single vector of unsigned char
that just counts all combinations.
But this would require at least one byte per ID and per number, thus 143 MByte.
The whole concept of using IDs could be replaced by temporarily allowing terms 2^a * 3^b > 10^6 and filtering them in the last step.
However, this further increases memory consumption.
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 | ./333
Output:
Note: the original problem's input 1000000
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 <set>
#include <algorithm>
// ---------- standard prime sieve 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) + 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 ----------
// slow algorithm
unsigned int bruteForce(unsigned int limit)
{
// create all terms 2^i * 3^j
std::vector<unsigned int> terms;
for (unsigned int two = 1; two <= limit; two *= 2)
for (unsigned int three = 1; two * three <= limit; three *= 3)
terms.push_back(two * three);
// each term is a multiple of 1, so I can safely remove it
terms.erase(terms.begin());
// in ascending order
std::sort(terms.begin(), terms.end());
// mainly for debugging: keep track of all terms of all solutions
typedef std::vector<unsigned int> Solution;
std::vector<std::vector<Solution>> solutions(limit + 1);
// analyze all subsets
unsigned int maxMask = 1 << terms.size();
for (unsigned int mask = 1; mask < maxMask; mask++)
{
unsigned int sum = 0;
Solution currentTerms;
// collect all terms
for (unsigned int pos = 0; pos < terms.size(); pos++)
{
unsigned int posMask = 1 << pos;
if ((posMask & mask) == 0)
continue;
// yes, use that term
sum += terms[pos];
currentTerms.push_back(terms[pos]);
}
// too large ?
if (sum > limit)
continue;
// no term must a multiple of any other term
bool isGood = true;
for (unsigned int i = 0; i < currentTerms.size() && isGood; i++)
for (unsigned int j = i + 1; j < currentTerms.size() && isGood; j++)
isGood = (currentTerms[j] % currentTerms[i] > 0); // note: currentTerms was sorted ascendingly
// all terms are coprime
if (isGood)
solutions[sum].push_back(currentTerms);
}
// add all primes with exactly one solution
unsigned int result = 0;
for (unsigned int i = 0; i <= limit; i++)
if (solutions[i].size() == 1 && isPrime(i))
{
// display terms
auto current = solutions[i].front();
std::cout << i << "=";
for (auto x : current)
std::cout << x << " ";
std::cout << std::endl;
result += i;
}
return result;
}
// return 2^two * 3^three
unsigned int power(unsigned int two, unsigned int three)
{
// fast exponentiation for 3^three
unsigned int result = 1;
unsigned int base = 3;
while (three > 0)
{
if (three & 1)
result *= base;
base *= base;
three >>= 1;
}
// multiplying by 2^two is a bitshift
return result << two;
}
int main()
{
unsigned int limit = 1000000;
std::cin >> limit;
// solve small inputs
//std::cout << bruteForce(limit) << std::endl;
// count how many powers of 2 and 3 are below the limit
unsigned int numPowerOfTwo = 0;
unsigned int numPowerOfThree = 0;
while (power(numPowerOfTwo, 0) <= limit)
numPowerOfTwo++;
while (power(0, numPowerOfThree) <= limit)
numPowerOfThree++;
// assign a unique ID to each 2^i * 3^j <= limit
const unsigned char TooLarge = 0;
unsigned char nextId = 1;
unsigned char ids[30][20]; // more than big enough
std::vector<std::pair<unsigned int, unsigned int>> reverseId(1); // 0 is dummy ID (TooLarge)
for (unsigned int expTwo = 0; expTwo < numPowerOfTwo; expTwo++)
for (unsigned int expThree = 0; expThree < numPowerOfThree; expThree++)
{
auto current = power(expTwo, expThree);
if (current <= limit)
{
reverseId.push_back({ expTwo, expThree });
ids[expTwo][expThree] = nextId++;
}
else
ids[expTwo][expThree] = TooLarge;
}
// => only 143 out of 260 pairs
// allocate memory, will store the number of ways to reach a value
std::vector<std::vector<bool>> reachable(nextId);
for (auto& x : reachable)
x.resize(limit + 1, 0);
std::vector<std::vector<bool>> multiple (nextId);
for (auto& x : multiple)
x.resize(limit + 1, 0);
// seed values
for (unsigned int id = 1; id < reverseId.size(); id++)
{
auto current = reverseId[id];
auto value = power(current.first, current.second);
reachable[id][value] = true;
}
// add all primes with unique partitions
fillSieve(limit);
unsigned int sum = 0;
for (unsigned int i = 1; i <= limit; i++)
{
// is there at least one combination to reach that number ?
bool possible = false;
// is there at most one combination to reach that number ?
bool unique = true;
// => look for numbers that are possible and unique
for (unsigned int id = 1; id < reverseId.size(); id++)
{
auto expTwo = reverseId[id].first;
auto expThree = reverseId[id].second;
// can I reach this number if one term has the current ID ?
if (!reachable[id][i])
continue;
// store new sum
if (possible || multiple[id][i]) // multiple ways to reach current position ?
unique = false;
possible = true;
// try to add every possible coprime number with 2^c * 3^d (where c > a and b < d)
for (unsigned int nextExpTwo = expTwo + 1; nextExpTwo < numPowerOfTwo; nextExpTwo++)
for (unsigned int nextExpThree = 0; nextExpThree < expThree; nextExpThree++)
{
auto nextId = ids[nextExpTwo][nextExpThree];
if (nextId == TooLarge)
break;
// add one term
auto next = i + power(nextExpTwo, nextExpThree);
if (next > limit)
break;
// if the next number was already reachable then there are now multiple ways
if (reachable[nextId][next])
multiple [nextId][next] = true;
else
reachable[nextId][next] = true;
}
}
// unique sum of a prime ?
if (possible && unique && isPrime(i))
sum += i;
}
// display result
std::cout << sum << std::endl;
return 0;
}
This solution contains 35 empty lines, 47 comments and 4 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 0.4 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
September 20, 2017 submitted solution
September 20, 2017 added comments
Difficulty
Project Euler ranks this problem at 35% (out of 100%).
Links
projecteuler.net/thread=333 - 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/euler313/euler333.py (written by Den Raskovalov)
Python github.com/Meng-Gen/ProjectEuler/blob/master/333.py (written by Meng-Gen Tsai)
Python github.com/smacke/project-euler/blob/master/python/333.py (written by Stephen Macke)
C++ github.com/roosephu/project-euler/blob/master/333.cpp (written by Yuping Luo)
C++ github.com/smacke/project-euler/blob/master/cpp/333.cpp (written by Stephen Macke)
C github.com/LaurentMazare/ProjectEuler/blob/master/e333.c (written by Laurent Mazare)
Java github.com/thrap/project-euler/blob/master/src/Java/Problem333.java (written by Magnus Solheim Thrap)
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 |
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 329 - Prime Frog | Maximix Arrangements - problem 336 >> |