Problem 250: $250250$
View on Project EulerProject 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
Take the \(N=250250\) numbers \(1^1,2^2,3^3,\dots,N^N\). We must count the non-empty subsets of this set whose elements add up to a multiple of \(250\), and give the count modulo \(10^{16}\), that is, its rightmost \(16\) digits. [1] The numbers \(i^i\) grow strictly with \(i\), so the set has \(N\) different elements, and a subset is the same as a set of indices \(i\). The statement gives no checkpoint; the C++ program checks that the first ten numbers have exactly \(5\) such subsets (Section 7).
The whole solution in one sentence: whether a sum is divisible by \(250\) depends only on the remainders of its terms after division by \(250\), so we go through the numbers once and keep, for each of the \(250\) possible remainders, the number of subsets of the numbers seen so far whose sum has that remainder, all counts taken modulo \(10^{16}\).
Mathematical Approach
Brute force is hopeless: \(N^N\) alone has \(1\,350\,944\) digits, and the number \(2^N-1\) of non-empty subsets has \(75\,333\) digits. We start from a small version that can be done by hand.
1. A small version by hand
Replace \(250\) by \(5\) and take the first five numbers \(1\), \(4\), \(27\), \(256\) and \(3125\). Exactly five of their \(31\) non-empty subsets have a sum divisible by \(5\):
$$\begin{aligned}&\{1,4\},\ \{4,256\},\ \{3125\},\\&\{1,4,3125\},\ \{4,256,3125\},\end{aligned}$$
with the sums \(5\), \(260\), \(3125\), \(3130\) and \(3385\). The remainders of the five numbers after division by \(5\), namely \(1,4,2,1,0\), already decide everything: \(4+256\) is divisible by \(5\) because the remainders \(4\) and \(1\) add up to \(5\).
2. Only remainders matter
We write \(a\bmod m\) for the remainder of \(a\) after division by \(m\), and \(a\equiv b\pmod m\) when \(m\) divides \(a-b\). Congruences may be added: if \(a_1\equiv b_1\) and \(a_2\equiv b_2\pmod m\), then \(a_1+a_2\equiv b_1+b_2\pmod m\). [2] Applying this rule once per term, a sum of any number of terms is congruent to the sum of their remainders. So, with
$$r_i=i^i\bmod 250,$$
a set \(A\) of indices is good exactly when \(\sum_{i\in A}r_i\equiv0\pmod{250}\). The first ten remainders are \(1,4,27,6,125,156,43,216,239,0\); for example \(8^8=16777216\)\({}=250\cdot67108+216\).
3. Remainders of huge powers
To find \(r_i\) we never compute \(i^i\). Write the exponent in binary and read its digits from the lowest, starting with the result \(1\) and \(c=i\bmod250\). Squaring \(c\) after every digit keeps \(c\equiv i^{2^j}\) at digit \(j\), and a digit \(1\) multiplies \(c\) into the result; the selected powers \(2^j\) add up to the exponent, so the result becomes \(i^i\). Congruences may also be multiplied, so each product is reduced modulo \(250\) at once; [2] [3] a product of two reduced remainders is at most \(249^2=62001\). For \(7^7\), with \(7=111_2\) and congruences modulo \(250\):
- first digit: the result becomes \(7\), and \(c\) becomes \(7^2=49\);
- second digit: the result becomes \(7\cdot49=343\equiv93\), and \(c\) becomes \(49^2=2401\equiv151\);
- third digit: the result becomes \(93\cdot151=14043\equiv43\).
Indeed \(7^7=823543=250\cdot3294+43\). An exponent \(i\le N\lt2^{18}\) has at most \(18\) binary digits, so each \(r_i\) costs at most \(18\) squarings and \(18\) further multiplications.
4. The remainders repeat every 500 numbers
The remainders repeat: \(r_{i+500}=r_i\) for every \(i\ge1\). As \(250=2\cdot125\) and the factors \(2\) and \(125\) are coprime, a remainder modulo \(250\) is fixed by the remainders modulo \(2\) and modulo \(125\), by the Chinese remainder theorem. [5] We check both.
- Modulo \(2\): \(i^i\) is odd exactly when \(i\) is odd, and \(i+500\) has the same parity as \(i\).
- Modulo \(125\), if \(5\mid i\): \(5^i\) divides \(i^i\) and \(i\ge5\), so \(i^i\equiv0\), and likewise for \(i+500\).
- Modulo \(125\), if \(5\nmid i\): \(100\) of the numbers \(1,\dots,125\) are prime to \(5\), so Euler's theorem gives \(i^{100}\equiv1\pmod{125}\), and the exponent may be changed by multiples of \(100\). [4] Thus \(i^i\bmod125\) depends only on \(i\bmod125\) and \(i\bmod100\), and adding \(500\) changes neither.
No shorter shift works: for every \(1\le s\lt500\) some \(i\le500\) has \(r_{i+s}\ne r_i\). As \(N=500\cdot500+250\), the list \(r_1,\dots,r_N\) is \(500\) copies of one period followed by \(r_1,\dots,r_{250}\), with only \(152\) different values. The three cases also show that \(r_i=0\) exactly when \(i\) is even and a multiple of \(5\), that is, when \(10\mid i\), and \(r_i=125\) exactly when \(i\) is an odd multiple of \(5\); each happens for \(25025\) indices.
5. Counting subsets by the remainder of their sum
For \(0\le k\le N\) and every integer \(t\), let \(c_k(t)\) be the number of subsets \(A\) of \(\{1,\dots,k\}\), the empty set included, with \(\sum_{i\in A}r_i\equiv t\pmod{250}\). It depends only on \(t\bmod250\), so the row \(c_k\) consists of the \(250\) numbers \(c_k(0),\dots,c_k(249)\). Only the empty set uses no number, and its sum is \(0\), so \(c_0(t)\) is \(1\) for \(t\equiv0\) and \(0\) otherwise. For \(k\ge1\),
$$c_k(t)=c_{k-1}(t)+c_{k-1}(t-r_k).$$
Proof: the subsets of \(\{1,\dots,k\}\) with sum \(\equiv t\) that do not contain \(k\) are exactly the subsets of \(\{1,\dots,k-1\}\) with sum \(\equiv t\). Removing \(k\) from those that contain \(k\) leaves a subset of \(\{1,\dots,k-1\}\) with sum \(\equiv t-r_k\), and putting \(k\) back undoes this, so they are counted by \(c_{k-1}(t-r_k)\). This is dynamic programming, with the same include-or-leave-out choice as in the classical algorithm for the subset sum problem. [6] The empty set does not count, so the answer is \(c_N(0)-1\), taken modulo \(10^{16}\). Each row adds up to \(2^k\), since each of the \(2^k\) subsets of \(\{1,\dots,k\}\) is counted once; this gives a cheap check.
6. The small version as a table of counts
Back to the divisor \(5\) and the remainders \(1,4,2,1,0\) of Section 1. Writing each row as \((c_k(0),\dots,c_k(4))\), the rule of Section 5 with \(5\) in place of \(250\) gives:
- \(k=0\): \((1,0,0,0,0)\);
- \(k=1\), remainder \(1\): \((1,1,0,0,0)\);
- \(k=2\), remainder \(4\): \((2,1,0,0,1)\);
- \(k=3\), remainder \(2\): \((2,2,2,1,1)\);
- \(k=4\), remainder \(1\): \((3,4,4,3,2)\);
- \(k=5\), remainder \(0\): \((6,8,8,6,4)\).
Each new row is the old row plus a copy of it moved \(r_k\) places to the right, wrapping around at the end. For example \(c_3(2)=c_2(2)+c_2(0)=0+2\): the empty set and \(\{1,4\}\), both with sum \(\equiv0\), become \(\{27\}\) and \(\{1,4,27\}\) when \(27\) is added. The rows add up to \(1,2,4,\dots,32\). The last number has remainder \(0\) and doubles every count, and \(c_5(0)-1=5\) gives the five subsets found by hand.
7. The first ten numbers with the divisor 250
With the remainders of Section 2, the rule gives \(c_7(0)=1\), so no non-empty subset of the first seven numbers works, and \(c_9(0)=3\): besides the empty set, the index sets \(\{1,3,4,8\}\) and \(\{1,2,4,9\}\) qualify, with the remainder sums \(1+27+6+216=250\) and \(1+4+6+239=250\). With the full numbers,
$$\begin{aligned}&1+27+256+16777216\\&\quad=250\cdot67110,\\&1+4+256+387420489\\&\quad=250\cdot1549683.\end{aligned}$$
The tenth number \(10^{10}\) has remainder \(0\) and doubles the counts, so \(c_{10}(0)=6\): the five non-empty subsets are \(\{10^{10}\}\) and the two sets above, each with and without \(10^{10}\). This is the value \(5\) that the C++ program tests.
8. Why the update needs the old row
Without reduction, one array would suffice if the sums were updated from the largest down, because every entry read would still hold its old value. Modulo \(250\) no such order exists when \(r=r_k\ne0\): each new value at \(t+r\) needs the old value at \(t\), all around the cycle \(t,t+r,t+2r,\dots\), which returns to \(t\). Whichever entry of the cycle is overwritten first, its old value is still needed one step further on, and using the new value there would count subsets that contain \(k\) twice. Only \(r_k=0\), which doubles every entry, allows a plain update in place; otherwise one would have to walk each cycle carrying an old value along. The programs simply keep the old row intact: C++ and Java write the new row into a second array and swap the two, and Python builds a new list for every row.
9. Counting modulo \(10^{16}\)
The counts are gigantic: every subset of the \(25025\) numbers with remainder \(0\) is good, so \(c_N(0)\ge2^{25025}\), a number with \(7534\) digits, and by Section 10 the count has \(75331\) digits. Congruences modulo \(M=10^{16}\) may be added as well, [2] and the recurrence only adds, so every count may be replaced by its remainder modulo \(M\) after each addition. The empty set is removed as \((c+M-1)\bmod M\), where \(c\) is the remainder of \(c_N(0)\); this gives \(M-1\), not \(-1\), if \(c=0\). Stored counts are below \(10^{16}\), so a sum of two of them is below \(2\cdot10^{16}\), far below \(2^{63}-1\approx9.2\cdot10^{18}\): 64-bit integers suffice, and one subtraction of \(M\) restores the range.
10. Generating functions and the size of the count
Let \(F_k(x)=\prod_{i=1}^{k}(1+x^{r_i})\), computed with \(x^{250}=1\), so that exponents are taken modulo \(250\). Expanding the product, each subset \(A\) contributes \(x^{s}\) with \(s=\sum_{i\in A}r_i\), so the coefficient of \(x^t\) is \(c_k(t)\). Multiplying by \(1+x^{r_k}\) adds a copy shifted by \(r_k\), because multiplying generating functions convolves their coefficients; [7] this is the recurrence of Section 5. The roots of unity filter extracts the exponents divisible by \(250\): [8] with \(\omega=e^{2\pi\mathrm{i}/250}\), the average of \(1,\omega^s,\dots,\omega^{249s}\) is \(1\) if \(250\mid s\) and \(0\) otherwise, as the geometric sum shows, so
$$c_N(0)=\frac1{250}\sum_{j=0}^{249}F_N(\omega^j).$$
The term \(j=0\) is \(2^N/250\). For odd \(j\), the factor \(1+x^{125}\) from \(i=5\) becomes \(1+\omega^{125j}=0\). For each of the \(124\) even \(j\ne0\), a direct check finds an index \(i\le16\) with \(84\le jr_i\bmod250\le166\); then the angle of \(\omega^{jr_i}\) lies between \(2\pi/3\) and \(4\pi/3\), so \(|1+\omega^{jr_i}|\le1\). By Section 4 this factor occurs at least \(500\) times in \(F_N\), and no factor exceeds \(2\) in absolute value, so \(|F_N(\omega^j)|\le2^{N-500}\). Hence \(c_N(0)\) differs from \(2^N/250\) by at most \(\tfrac{124}{250}\,2^{N-500}\lt2^{N-501}\), a tiny fraction of it, and exact comparisons with powers of \(10\) show that it has \(75331\) digits. Its last \(16\) digits, however, are invisible to this estimate.
11. A shortcut through the period
By Section 4, \(F_N=P^{500}\,Q\), where \(P\) is the product of the factors \(1+x^{r_i}\) for \(1\le i\le500\) and \(Q\) the product for \(1\le i\le250\). Building \(P\) and \(Q\) takes \(750\) steps of the recurrence. Repeated squaring as in Section 3, without the unused square after the last digit, gives \(P^{500}\) with \(8\) squarings and \(6\) multiplications, since \(500=111110100_2\), and a final product with \(Q\) gives \(F_N\). Each of these \(15\) products needs \(250^2=62500\) multiplications and additions of coefficients, about a million of each in all, against \(62.5\) million additions in the row-by-row recurrence; but products of coefficients below \(10^{16}\) need 128-bit arithmetic or big integers. The programs do not use it; it serves as an independent check.
How the Code Works
The C++ program [9] keeps its settings in Options: limit \(=250250\), divisor \(=250\), modulo \(=10^{16}\) and the flag run_checkpoints. parse_arguments accepts --skip-checkpoints, --limit=, --divisor= and --mod=; any other argument, a limit or divisor below \(1\), or a modulus \(0\) ends the program with exit code \(1\). The parsers do not check for overflow, and a custom modulus above \(2^{63}\) can make an addition wrap around (Section 9). Unless the checkpoints are skipped, run_checkpoints requires solve to return \(5\) for the limit \(10\), the divisor \(250\) and the modulus \(10^{16}\) (Section 7), and otherwise ends the program with exit code \(2\). Then main prints the result of solve.
powmod_int(base, exp, mod) is the method of Section 3: for every binary digit of exp, from the lowest, it multiplies result by cur if the digit is \(1\) and squares cur, forming each product in long long. solve(limit, divisor, modulo) holds two vectors dp and next of divisor counts of the 64-bit unsigned type u64 and starts with dp[0] = 1, the row \(c_0\). For each i it computes residue \(=r_i\), copies dp into next for the subsets without \(i\), and adds every dp[r] to next[(r + residue) % divisor]: the subsets counted in dp[r] receive \(i\) and move to the remainder \(r+r_i\). This is the recurrence of Section 5 read forwards, with \(t=r+r_i\). A comparison and at most one subtraction keep each entry below modulo (Section 9), and dp.swap(next) exchanges the arrays without copying them (Section 8). Finally solve returns (dp[0] + modulo - 1) % modulo.
The Python solve has the same constants and loop. It gets \(r_i\) from the built-in pow(i, i, divisor), which never builds \(i^i\), [10] builds each new row as new_dp = list(dp) and reduces with % modulo; Python integers have no fixed size. The Java class has powmodInt, a copy of powmod_int, and solve with two long arrays that are exchanged by swapping references; a long reaches \(2^{63}-1\), [10] enough by Section 9. Its inner loop skips entries equal to \(0\); a replay of the computation shows that this happens only in the first \(12\) steps, while some remainders are still unreachable. Neither version has command-line options or the checkpoint.
Complexity and Verification
For \(N\) numbers and the divisor \(m=250\), the recurrence makes \(Nm=62\,562\,500\) additions of counts, of which Java skips \(1893\), each followed by a comparison in C++ and Java and by a reduction with % in Python. The power loops of C++ and Java make \(6\,466\,021\) multiplications modulo \(250\), at most \(35\) for one \(i\); Python leaves this work to pow. The time is \(O(Nm+N\log N)\) and the memory \(O(m)\): two arrays of \(250\) counts of \(8\) bytes, \(4000\) bytes of data, in C++ and Java, and two lists of references to integer objects in Python. On our machine, start-up included, the C++ and Java programs take about \(0.3\) seconds and the Python program about \(4\) seconds.
The C++ program checks the value \(5\) for the first ten numbers on every run unless --skip-checkpoints is given. For this article we also ran these checks, which are not part of the published programs:
- A brute force that adds up the actual numbers \(i^i\) of all \(2^N\) subsets agrees with the recurrence for every \(1\le N\le21\) with the divisor \(250\), and for every \(1\le N\le14\) with every divisor from \(1\) to \(60\), \(861\) cases in all. For \(N=20\) it finds \(4239\) subsets, as does the C++ program with
--limit=20. - The period \(500\) holds for every \(i\le N\), and so do the finite checks of Sections 4 and 10.
- At full size, the shortcut of Section 11 gives the same \(250\) counts modulo \(10^{16}\) as the recurrence, and they add up to \(2^N\) modulo \(10^{16}\).
- The three programs, with the C++ compiled by GCC and by Clang, print the same number as a separate implementation of the recurrence.
Footnotes and References
The references below supply the background facts; the period, the examples, the size estimate and the operation counts were worked out for this article.
- Project Euler 250 — 250250. The official statement asks for the number of non-empty subsets of {1¹, 2², 3³, …, 250250²⁵⁰²⁵⁰} whose sum of elements is divisible by 250, and for the rightmost 16 digits of that number.
- Wikipedia — Modular arithmetic. Congruence is compatible with addition and multiplication: if a₁ ≡ b₁ and a₂ ≡ b₂ (mod m), then a₁ + a₂ ≡ b₁ + b₂ and a₁a₂ ≡ b₁b₂ (mod m).
- Wikipedia — Modular exponentiation. The right-to-left binary method squares a running base once per bit, multiplies it into the result for each bit 1 and reduces every product modulo the modulus.
- Wikipedia — Euler's theorem. For coprime a and n, a^φ(n) ≡ 1 (mod n), with Euler's totient φ, so a^x ≡ a^y (mod n) whenever x ≡ y (mod φ(n)); here φ(125) = 100.
- Wikipedia — Chinese remainder theorem. The remainders of an integer after division by pairwise coprime integers determine its remainder after division by their product uniquely.
- Wikipedia — Subset sum problem. Its pseudo-polynomial dynamic programming algorithm passes from the first i numbers to the first i + 1 by leaving out or including the next one.
- Wikipedia — Generating function. Multiplying ordinary generating functions yields the discrete convolution, the Cauchy product, of their coefficient sequences.
- Wikipedia — Series multisection. The roots of unity filter: the terms of f(z) with exponents ≡ p (mod q) add up to (1/q) Σ ω^(−kp) f(ωᵏz) over k = 0, …, q − 1, where ω = e^(2πi/q).
- C++ — Euler250.cpp. The linked immutable C++ revision contains
parse_arguments,powmod_int,solveandrun_checkpoints. It is the implementation record for this article. - Python — built-in function pow: with three arguments it returns base to the power exp modulo mod, computed more efficiently than pow(base, exp) % mod. Java Language Specification, SE 17, §4.2.1: a long holds the integers from −2⁶³ to 2⁶³ − 1.
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());
}
}