Detailed mathematical approach

Problem Summary

Consider the sequence

$a_i=i^i\qquad(1\le i\le N),$

with \(N=250250\) in the main problem. We must count the non-empty subsets \(S\subseteq\{1,\dots,N\}\) such that

$\sum_{i\in S} a_i \equiv 0 \pmod{250},$

and report the result modulo

$M=10^{16}.$

The direct search space has size \(2^N-1\), so brute force is impossible. The key observation is that divisibility by \(250\) depends only on residue classes modulo \(250\), not on the full size of the numbers \(i^i\).

Mathematical Approach

Let

$r_i=a_i\bmod 250=i^i\bmod 250.$

Then for any subset \(S\),

$\sum_{i\in S} a_i \equiv \sum_{i\in S} r_i \pmod{250}.$

So the original problem is equivalent to counting subsets whose residue sum is \(0\) modulo \(250\).

Why Reduction Modulo \(250\) Is Sufficient

The congruence rule

$x\equiv y\pmod{250}\quad\Longrightarrow\quad x+z\equiv y+z\pmod{250}$

shows that divisibility of a sum by \(250\) depends only on each summand's class in \(\mathbb{Z}/250\mathbb{Z}\). Therefore the gigantic integers \(i^i\) never need to be stored explicitly. Only the 250 possible residues matter.

Subset-Sum DP on the Cyclic Group \(\mathbb{Z}/250\mathbb{Z}\)

For \(k=0,1,\dots,N\), define

$dp_k(t)=\#\left\{S\subseteq\{1,\dots,k\}:\sum_{i\in S} r_i\equiv t\pmod{250}\right\}.$

The desired count is \(dp_N(0)-1\), because \(dp_N(0)\) includes the empty subset.

Initial Condition

Before processing any terms, there is exactly one subset of \(\{1,\dots,0\}\): the empty subset. Its sum is \(0\). Hence

$dp_0(0)=1,\qquad dp_0(t)=0\quad(1\le t\le 249).$

Transition Formula

When we process residue \(r_k\), every old subset has two possibilities:

1. Exclude \(k\): its residue class stays the same.

2. Include \(k\): its residue class shifts by \(r_k\) modulo \(250\).

Therefore

$dp_k(t)=dp_{k-1}(t)+dp_{k-1}(t-r_k)\pmod{M},$

where \(t-r_k\) is interpreted modulo \(250\).

The code implements the same recurrence in forward form with a temporary array next:

$next[t]=dp[t],$

$next[(t+r_k)\bmod 250]\leftarrow next[(t+r_k)\bmod 250]+dp[t]\pmod{M}.$

Using a second array is important: it guarantees 0/1 behavior. If we updated dp in place, a newly created subset could be reused in the same iteration, which would incorrectly allow the same element to be chosen multiple times.

Generating-Function View

The same recurrence can be written as a product of polynomials modulo \(x^{250}-1\):

$F_N(x)=\prod_{i=1}^{N}\left(1+x^{r_i}\right)\pmod{x^{250}-1}.$

The coefficient of \(x^t\) in \(F_N(x)\) counts subsets whose residue sum is \(t\) modulo \(250\). The dynamic program is simply a coefficient-update algorithm for this product.

Computing \(i^i\bmod 250\) Efficiently

The code uses repeated squaring to compute

$r_i=i^i\bmod 250$

in \(O(\log i)\) time. This is the standard modular exponentiation recurrence:

$b^{2m}\equiv (b^2)^m\pmod{250},\qquad b^{2m+1}\equiv b\cdot (b^2)^m\pmod{250}.$

Because multiplication is reduced modulo \(250\) at every step, intermediate numbers stay small.

Worked Example: First 10 Terms

For \(N=10\), the residues are

$1^1,2^2,\dots,10^{10}\equiv 1,4,27,6,125,156,43,216,239,0 \pmod{250}.$

The checkpoint in the code states that the answer is \(5\). Indeed, the non-empty divisible subsets are:

$\{1,3,4,8\},\qquad \{1,2,4,9\},\qquad \{10\},$

$\{1,3,4,8,10\},\qquad \{1,2,4,9,10\}.$

For example,

$1^1+3^3+4^4+8^8\equiv 1+27+6+216=250\equiv 0\pmod{250},$

and since \(10^{10}\equiv 0\pmod{250}\), appending the 10th term preserves divisibility. Hence there are exactly \(5\) valid non-empty subsets.

Why the Final Subtraction Is Necessary

The DP starts from the empty subset, so \(dp_N(0)\) counts all subsets with residue \(0\), including

$S=\varnothing.$

Since the problem asks only for non-empty subsets, the final answer is

$\boxed{(dp_N(0)-1)\bmod M.}$

Complexity Analysis

There are \(N\) iterations, and each iteration updates all \(250\) residue classes. Hence the dynamic-programming work is

$O(N\cdot 250)=O(ND).$

The memory usage is \(O(250)\), since only the current and next DP rows are stored. Modular exponentiation adds \(O(\log i)\) work per term, so the full cost is \(O(ND+N\log N)\); for fixed \(D=250\), this is still easily practical.

How the Code Works

The C++ solution initializes dp[0]=1 and all other states to zero. For each \(i\) from \(1\) to limit, it computes residue = powmod_int(i, i, divisor), copies dp into next, applies the include-transition to every residue class, and then swaps the arrays. After the loop, dp[0] contains the number of divisible subsets including the empty one, so the code returns (dp[0] + modulo - 1) % modulo. The Python and Java versions implement the same recurrence.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=250
2. Dynamic programming: Wikipedia — Dynamic programming
3. Modular exponentiation: Wikipedia — Modular exponentiation
4. Cyclic groups and modular arithmetic: Wikipedia — Modular arithmetic
5. Generating functions in combinatorics: Wikipedia — Generating function