Problem 250: $250250$

Project Euler Problem 250 Solution

EulerSolve provides an optimized solution for Project Euler Problem 250, $250250$, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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

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

Problem 250 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int limit = 250250;
    int divisor = 250;
    u64 modulo = 10000000000000000ULL;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit) ||
            parse_int_after_prefix(arg, "--divisor=", options.divisor) ||
            parse_u64_after_prefix(arg, "--mod=", options.modulo)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 1 && options.divisor >= 1 && options.modulo > 0;
}

int powmod_int(int base, int exp, const int mod) {
    int result = 1 % mod;
    int cur = base % mod;

    while (exp > 0) {
        if (exp & 1) {
            result = static_cast<int>((static_cast<long long>(result) * cur) % mod);
        }
        cur = static_cast<int>((static_cast<long long>(cur) * cur) % mod);
        exp >>= 1;
    }

    return result;
}

u64 solve(const int limit, const int divisor, const u64 modulo) {
    std::vector<u64> dp(static_cast<std::size_t>(divisor), 0);
    std::vector<u64> next(static_cast<std::size_t>(divisor), 0);
    dp[0] = 1;

    for (int i = 1; i <= limit; ++i) {
        const int residue = powmod_int(i, i, divisor);

        for (int r = 0; r < divisor; ++r) {
            next[static_cast<std::size_t>(r)] = dp[static_cast<std::size_t>(r)];
        }

        for (int r = 0; r < divisor; ++r) {
            const int to = (r + residue) % divisor;
            u64& cell = next[static_cast<std::size_t>(to)];
            cell += dp[static_cast<std::size_t>(r)];
            if (cell >= modulo) {
                cell -= modulo;
            }
        }

        dp.swap(next);
    }

    const u64 all_divisible = dp[0];
    return (all_divisible + modulo - 1) % modulo;  // Exclude empty subset.
}

bool run_checkpoints() {
    if (solve(10, 250, 10000000000000000ULL) != 5ULL) {
        std::cerr << "Checkpoint failed for first 10 values" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit, options.divisor, options.modulo) << '\n';
    return 0;
}

Python

def solve():
    limit = 250250
    divisor = 250
    modulo = 10**16

    dp = [0] * divisor
    dp[0] = 1

    for i in range(1, limit + 1):
        residue = pow(i, i, divisor)
        new_dp = list(dp)
        for r in range(divisor):
            to = (r + residue) % divisor
            new_dp[to] = (new_dp[to] + dp[r]) % modulo
        dp = new_dp

    return str((dp[0] + modulo - 1) % modulo)

if __name__ == '__main__':
    print(solve())

Java

public class Euler250 {
    static int powmodInt(int base, int exp, int mod) {
        int result = 1 % mod;
        int cur = base % mod;

        while (exp > 0) {
            if ((exp & 1) != 0) {
                result = (int) (((long) result * cur) % mod);
            }
            cur = (int) (((long) cur * cur) % mod);
            exp >>= 1;
        }
        return result;
    }

    public static String solve() {
        int limit = 250250;
        int divisor = 250;
        long modulo = 10000000000000000L;

        long[] dp = new long[divisor];
        long[] next = new long[divisor];
        dp[0] = 1;

        for (int i = 1; i <= limit; ++i) {
            int residue = powmodInt(i, i, divisor);

            for (int r = 0; r < divisor; ++r) {
                next[r] = dp[r];
            }

            for (int r = 0; r < divisor; ++r) {
                if (dp[r] == 0)
                    continue;
                int to = (r + residue) % divisor;
                next[to] += dp[r];
                if (next[to] >= modulo) {
                    next[to] -= modulo;
                }
            }

            long[] temp = dp;
            dp = next;
            next = temp;
        }

        long allDivisible = dp[0];
        long ans = (allDivisible + modulo - 1) % modulo;
        return String.valueOf(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}