Project Euler Solution #48: Self powers

Problem 48

The series, 1¹ + 2² + 3³ + ... + 10¹⁰ = 10405071317.

Find the last ten digits of the series, 1¹ + 2² + 3³ + ... + 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

\[ b^e \pmod{m} \equiv \begin{cases} {(b\bmod m)} , &\text{if }e=1\\ {(b^2\bmod m)}^{\frac {e-1}{2}}, &\text{if }e\text{ is odd}\\ {(b^2\bmod m)}^{\frac{e}{2}} , &{\text{if }}e{\text{ 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.