<< problem 425 - Prime connection | Unfair wager - problem 436 >> |
Problem 429: Sum of squares of unitary divisors
(see projecteuler.net/problem=429)
A unitary divisor d of a number n is a divisor of n that has the property gcd(d, n/d) = 1.
The unitary divisors of 4! = 24 are 1, 3, 8 and 24.
The sum of their squares is 1^2 + 3^2 + 8^2 + 24^2 = 650.
Let S(n) represent the sum of the squares of the unitary divisors of n. Thus S(4!)=650.
Find S(10^8!) mod (10^9+9).
My Algorithm
The main insight regarding unitary divisors is that if I have a prime factorization of n! then any prime factor can only be found
in either d or n!/d but not both because otherwise they would share the prime factor p and therefore gcd(d, n!/d) >= p.
Mathworld has the formula including some explanations (see mathworld.wolfram.com/UnitaryDivisorFunction.html):
(1) sigma^{*}_2(p^{a_1}_1, p^{a_2}_2, ...) = (1 + p^{2a_1}_1) * (1 + p^{2a_2}_2) * ...
The resulting products will be huge and I have to mod 10^9 + 9 quite often.
I copied the powmod()
function form my toolbox instead of computing p^{2a_i}_i mod (10^9+9) in a loop.
Now all I need is a fast prime factorization of n! ... and it is pretty easy when using Legendre's formula (see en.wikipedia.org/wiki/Legendre's_formula):
(2) v_p(n!) = sum^{infinity}_{i=1} \lfloor dfrac{n}{p^i} \rfloor
In plain English: to figure how often a prime p can be found in n! I have to compute dfrac{n}{p} + dfrac{n}{p^2} + dfrac{n}{p^3} + ... where each division is an integer division (skip the remainder).
For example v_2(100!) = dfrac{100}{2} + dfrac{100}{4} + dfrac{100}{8} + dfrac{100}{16} + dfrac{100}{32} + dfrac{100}{64} = 50 + 25 + 12 + 6 + 3 + 1 = 97
That means 100! is a multiple of 2^97.
The inner loop of solve
starts with power = p
and in each iteration multiplies it by p
to get p*p
, then p*p*p
and so on.
A little trick to accelerate the process is to abort early when p
is a large prime 2p > n because then I know that dfrac{limit}{p} = 1 (about 30% faster due to avoiding slow division).
Note
About half of the time is spent in my prime sieve.
I could have merged the prime sieve and solve()
but this would mess up the code and I seriously doubt that a significant performance boost could be observed.
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 | ./429
Output:
Note: the original problem's input 100000000
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>
// ---------- 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;
}
}
}
// ---------- powmod is part of my toolbox, too ----------
// return (base^exponent) % modulo for 32-bit values, no need for mulmod
unsigned int powmod(unsigned int base, unsigned int exponent, unsigned int modulo)
{
unsigned int result = 1;
while (exponent > 0)
{
// fast exponentation:
// odd exponent ? a^b = a*a^(b-1)
if (exponent & 1)
result = (result * (unsigned long long)base) % modulo;
// even exponent ? a^b = (a*a)^(b/2)
base = (base * (unsigned long long)base) % modulo;
exponent >>= 1;
}
return result;
}
// ---------- problem-specific code ----------
unsigned long long solve(unsigned int limit, unsigned int modulo)
{
unsigned long long result = 1;
// find all prime numbers
fillSieve(limit);
for (unsigned int p = 2; p <= limit; p++)
{
// prime factors only
if (!isPrime(p))
continue;
// will be p^1, p^2, p^3, p^4, ...
unsigned long long power = p;
// count how often current prime appears in factorial(limit)
unsigned int count = 0;
while (power <= limit)
{
// fast check for large powers
if (2*power > limit)
{
count++;
break;
} // note: this if-block can be removed and you still get the correct result - albeit slower
// count multiples
auto multiples = limit / power;
count += multiples;
power *= p;
}
// (p*p)^x = p^2x
// => p*p can be 64 bits and powmod for 64 bits is more complex and much slower
result *= 1 + powmod(p, 2*count, modulo);
result %= modulo;
}
return result;
}
int main()
{
const unsigned int Modulo = 1000000009;
unsigned int limit = 100000000; // 10^8
std::cin >> limit;
std::cout << solve(limit, Modulo) << std::endl;
return 0;
}
This solution contains 20 empty lines, 26 comments and 2 preprocessor commands.
Benchmark
The correct solution to the original Project Euler problem was found in 0.8 seconds on an Intel® Core™ i7-2600K CPU @ 3.40GHz.
Peak memory usage was about 8 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 13, 2017 submitted solution
September 13, 2017 added comments
Difficulty
Project Euler ranks this problem at 20% (out of 100%).
Links
projecteuler.net/thread=429 - the best forum on the subject (note: you have to submit the correct solution first)
Code in various languages:
Python github.com/nayuki/Project-Euler-solutions/blob/master/python/p429.py (written by Nayuki)
Python github.com/smacke/project-euler/blob/master/python/429.py (written by Stephen Macke)
C++ github.com/evilmucedin/project-euler/blob/master/euler429/euler429.cpp (written by Den Raskovalov)
C++ github.com/Meng-Gen/ProjectEuler/blob/master/429.cc (written by Meng-Gen Tsai)
C++ github.com/roosephu/project-euler/blob/master/429.cpp (written by Yuping Luo)
Java github.com/nayuki/Project-Euler-solutions/blob/master/java/p429.java (written by Nayuki)
Java github.com/thrap/project-euler/blob/master/src/Java/Problem429.java (written by Magnus Solheim Thrap)
Mathematica github.com/nayuki/Project-Euler-solutions/blob/master/mathematica/p429.mathematica (written by Nayuki)
Mathematica github.com/steve98654/ProjectEuler/blob/master/429.nb
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 |
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 425 - Prime connection | Unfair wager - problem 436 >> |