<< problem 425 - Prime connection | Double pandigital number divisible by 11 - problem 491 >> |

# 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 to`echo 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 a 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=c++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)

# 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 |

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 |

I scored 13,386 points (out of 15600 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 | Double pandigital number divisible by 11 - problem 491 >> |