<< problem 401 - Sum of squares of divisors | Prime connection - problem 425 >> |

# Problem 407: Idempotents

(see projecteuler.net/problem=407)

If we calculate a^2 mod 6 for 0 <= a <= 5 we get: 0,1,4,3,4,1.

The largest value of a such that a^2 == a mod 6 is 4.

Let's call M(n) the largest value of a < n such that a^2 == a (mod n).

So M(6) = 4.

Find sum{M(n)} for 1 <= n <= 10^7.

# My Algorithm

The function `slow()`

just iterates over all a - that's nice for small n but obviously to slow for large n.

However, it reveals a few properties:

- a is often close to n

- if n is a prime number then a = 1

- if n has only one prime factor, that means is a perfect square/cube/... of a prime then a = 1 as well

Actually, those small solutions can be found online, too: oeis.org/A182665

Let's rewrite the equation a^2 == a mod n

(1) a^2 - a == 0 mod n

(2) a(a-1) == 0 mod n

The program has to find the largest a < n where a(a-1) is a multiple of n.

Each prime factor of n has to be either a prime factor of a or a - 1.

My program precomputes the largest prime factor of each n and stores it in `largestPrime`

(see `findLargestPrime`

).

Whenever `search(n)`

is started, it looks up the largest prime factor of `n`

.

If `largestPrime[n] == n`

then either n is prime or a power of a prime (e.g. 32 = 2^5 therefore `largestPrime[32] = 32`

).

Strictly speaking I don't need that special treatment for powers of a prime but it gives a little performance boost.

Then I analyze all multiples x of n's largest prime p, beginning with the largest multiple:

- if x = a is a multiple of p then check whether x^2 == x mod n is true → if yes, then x is the largest a

- if x = a - 1 → x + 1 = a is a multiple of p then check (x + 1)^2 == x mod n is true → if yes, then x + 1 is the largest a

## Note

I could reduce memory consumption by 50% if `largestPrime[n] = n / largestPrime`

and repeat the same division in `search()`

because then `unsigned short`

is sufficient.

Or compute `largestPrime`

in a segmented way (it's similar to a segmented prime sieve anyway).

Currently my solution of problem 407 is the "most expensive": it consumes lots of memory and CPU power - but both are within Project Euler's limit.

# 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 100 | ./407`

Output:

*Note:* the original problem's input `10000000`

__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>
// find largest prime factor of each number up to "limit"

std::vector<unsigned int> largestPrime;
// fill "largestPrime"

void findLargestPrimes(unsigned int limit)
{
largestPrime.resize(limit + 1, 1);
for (unsigned long long i = 2; i <= limit; i++)
{
// skip if i is not a prime
if (largestPrime[i] > 1)
continue;
// modify all multiples
for (auto j = i; j <= limit; j += i)
largestPrime[j] = i;
// and all exponents i^2, i^3, i^4, ...
auto j = i * i;
while (j <= limit)
{
largestPrime[j] = j;
j *= i;
}
}
}
// find maximum a that solves a^2 == a mod n

unsigned int search(unsigned int n)
{
// don't forget about that special case n = 1 ...
if (n == 1)
return 0;
// always return 1 if n has only one prime factor (is prime or perfect power of a prime, e.g. 32 = 2^5)
auto p = largestPrime[n];
if (p == n)
return 1;
// largest multiple of n's largest prime factor
auto start = n - (n % p);
// go backwards through all multiples of step
for (unsigned long long x = start; x >= p; x -= p)
{
// check a - 1 = x which is a = x + 1
auto next = x + 1;
if ((next * next) % n == next)
return next;
// check a
if ((x * x) % n == x)
return x;
}
return 1;
}
// simple algorithm (pretty slow)

unsigned int slow(unsigned int n)
{
unsigned long long a = n - 1;
while (a >= 1)
{
// found a solution ?
if ((a * a) % n == a)
break;
// no, keep going ...
a--;
}
return a;
}
int main()
{
unsigned int limit = 10000000;
std::cin >> limit;
// find largest prime factor of each number up to "limit"
findLargestPrimes(limit);
unsigned long long sum = 0;
for (unsigned int n = 1; n <= limit; n++)
{
//auto current = slow(n);
auto current = search(n);
sum += current;
}
// display result
std::cout << sum << std::endl;
return 0;
}

This solution contains 16 empty lines, 18 comments and 2 preprocessor commands.

# Benchmark

The correct solution to the original Project Euler problem was found in 47 seconds on a Intel® Core™ i7-2600K CPU @ 3.40GHz.

Peak memory usage was about 41 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

August 31, 2017 submitted solution

August 31, 2017 added comments

# Difficulty

Project Euler ranks this problem at **20%** (out of 100%).

# Links

projecteuler.net/thread=407 - **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 401 - Sum of squares of divisors | Prime connection - problem 425 >> |