Problem 48 — Self Powers

Problem Analysis

The task is to find the last ten decimal digits of the self-power sum 1¹ + 2² + ... + 1000¹⁰⁰⁰. The pivotal observation is that "the last ten digits" of an integer is precisely that integer reduced modulo 10¹⁰. Modular arithmetic is the arithmetic of remainders: rather than tracking a full value, I track only its remainder after division by a fixed modulus m. The crucial property is that the modulus distributes over both addition and multiplication, so (a + b) mod m = ((a mod m) + (b mod m)) mod m, and the same holds for products. This means I never need the true sum at all - the largest term, 1000¹⁰⁰⁰, has roughly 3000 decimal digits - because I can keep every intermediate value below 10¹⁰.

The problem then reduces to computing iⁱ mod 10¹⁰ for each i from 1 to 1000. Multiplying i by itself i times would cost O(i) multiplications per term. The standard speedup is modular exponentiation via exponentiation by squaring, which decomposes the exponent into its binary digits and computes the power in O(log i) multiplications by repeatedly squaring the base. Across all 1000 terms the total work is O(N log N) - effectively instantaneous - and every reduced value fits comfortably in a machine integer.

Solution Approach

Modular Exponentiation with a Reduced Running Sum

Both the Python and C files implement the same algorithm. I iterate i from 1 to the limit inclusive, compute the term iⁱ mod 10¹⁰ by fast modular exponentiation, add it to a running total, and reduce that total modulo 10¹⁰ after every addition. The final running total is the answer.

In Python I rely on the built-in three-argument pow(i, i, modulo). This is not the same as pow(i, i) % modulo; the three-argument form performs the modular reduction inside the squaring loop, so the full multi-thousand-digit power is never materialised. Python integers are arbitrary precision, so there is no overflow to consider, and the routine runs in optimised C under the hood.

In C I wrote mod_pow by hand as a classic binary exponentiation by squaring loop: while the exponent is positive, I test its lowest bit with exp & 1 to decide whether to fold the current base into the result, then square the base and right-shift the exponent. The subtlety is overflow. With a modulus of 10¹⁰, two reduced operands can each approach 10¹⁰, so their product can reach about 10²⁰ - well beyond the roughly 9.2x10¹⁸ ceiling of a signed 64-bit long long. I avoid this by casting one operand to __int128 before each multiply, computing the product at 128-bit precision, and immediately reducing modulo mod back into 64-bit range. The 128-bit value is purely transient and never stored across iterations.

Both versions accept the limit as a runtime argument rather than hard-coding 1000, so the same binary or script handles the small development case (limit 10, answer 405071317) and the full case (limit 1000, answer 9110846700) through the identical code path.