# Project Euler 48: Self powers

The series, 1^{1} + 2^{2} + 3^{3} + ... + 10^{10} = 10405071317.

Find the last ten digits of the series, 1^{1} + 2^{2} + 3^{3} + ... + 1000^{1000}.

## Solution

With the knowledge of how a sum and a product behaves under a modular congruence relation, the problem can be solved quickly:

\[\sum\limits_{i=1}^n i^i \equiv \sum\limits_{i=1}^n (i^i \bmod m) \pmod{m}\]

Now we have isolated the exponential term \(i^i \bmod m\), which can be solved under a modular congruence like this:

\[i^i \equiv \prod\limits_{j=1}^i i \equiv \prod\limits_{j=1}^i (i\bmod m) \pmod{m}\]

That's enough already to solve this problem, but we can optimize the inner loop by using Exponentiation by squaring, by stating the products recursively so that

\[i^i \pmod{m} \equiv {\begin{cases}(i\bmod m)\,(i^2\bmod m)^{\frac {i-1}{2}},&{\mbox{if }}i{\mbox{ is odd}}\\(i^2\bmod m)^{\frac {i}{2}},&{\mbox{if }}i{\mbox{ is even}}.\end{cases}}\]

By checking the least significant bit of the exponent, we can decide which case must be used so that the algorithm can be implemented like this in Python:

def modpow(b, e, m): r = 1 while e > 0: if e & 1: r = (r * b) % m b = (b * b) % m e = e >> 1 return r

Using the modpow function, the formula stated before can be implemented easily with a proper chosen \(m\) to get the last 10 digits:

def powersum(n, m): s = 0 for i in range(1, n + 1): s+= modpow(i, i, m) s%= m return s print powersum(1000, 10000000000)

The `modpow` implementation isn't necessary, since Python's `pow` function has the same capabilities when passing three parameters. However, calling library functions is only half of the fun.