Problem 553: Power Sets of Power Sets
View on Project EulerProject Euler Problem 553 Solution
EulerSolve provides an optimized solution for Project Euler Problem 553, Power Sets of Power Sets, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \([n]=\{1,\dots,n\}\). We choose a non-empty family \(\mathcal{F}\) of non-empty subsets of \([n]\). Its intersection graph has one vertex for each chosen subset, and two vertices are adjacent exactly when the corresponding subsets intersect. We must compute \(C(n,K)\), the number of such families whose intersection graph has exactly \(K\) connected components. The target input is \(n=10000\) and \(K=10\), so direct enumeration is completely infeasible. Mathematical Approach The implementations solve the problem by first counting all families on a fixed label set, then extracting the connected ones with exponential generating functions. Step 1: Count every non-empty family on \(r\) labels On a ground set of size \(r\), there are \(2^r-1\) non-empty subsets. Any non-empty selection of those subsets is a valid family, so the raw count is $$R_r=2^{2^r-1}-1.$$ Even before connectivity is considered, this quantity grows doubly exponentially, which explains why brute force is hopeless for the real parameters. Step 2: Force every label to appear Let \(U_n\) be the number of families whose union is exactly \([n]\). Inclusion-exclusion removes families that miss at least one label. If we restrict attention to a chosen set of \(r\) labels, then every selected subset must lie inside that \(r\)-set, so there are \(R_r\) possibilities....
Detailed mathematical approach
Problem Summary
Let \([n]=\{1,\dots,n\}\). We choose a non-empty family \(\mathcal{F}\) of non-empty subsets of \([n]\). Its intersection graph has one vertex for each chosen subset, and two vertices are adjacent exactly when the corresponding subsets intersect. We must compute \(C(n,K)\), the number of such families whose intersection graph has exactly \(K\) connected components. The target input is \(n=10000\) and \(K=10\), so direct enumeration is completely infeasible.
Mathematical Approach
The implementations solve the problem by first counting all families on a fixed label set, then extracting the connected ones with exponential generating functions.
Step 1: Count every non-empty family on \(r\) labels
On a ground set of size \(r\), there are \(2^r-1\) non-empty subsets. Any non-empty selection of those subsets is a valid family, so the raw count is
$$R_r=2^{2^r-1}-1.$$
Even before connectivity is considered, this quantity grows doubly exponentially, which explains why brute force is hopeless for the real parameters.
Step 2: Force every label to appear
Let \(U_n\) be the number of families whose union is exactly \([n]\). Inclusion-exclusion removes families that miss at least one label. If we restrict attention to a chosen set of \(r\) labels, then every selected subset must lie inside that \(r\)-set, so there are \(R_r\) possibilities. Therefore
$$U_n=\sum_{r=0}^{n}(-1)^{n-r}\binom{n}{r}R_r.$$
For the generating-function step we normalize by factorials:
$$a_n=\frac{U_n}{n!}\quad (n\ge 1),\qquad a_0=1.$$
The special value \(a_0=1\) represents the empty assembly of connected components; it is a standard bookkeeping convention for the exponential formula.
Step 3: Why the exponential formula applies
Consider a family whose union is all of \([n]\). If two graph components shared a label \(t\), then some chosen subset in the first component and some chosen subset in the second would both contain \(t\). Those two subsets would intersect, producing an edge between the components, which is impossible. Hence different connected components use disjoint blocks of labels.
Conversely, if we partition \([n]\) into disjoint non-empty blocks and place one connected family on each block, then the resulting intersection graph has exactly those blocks as its connected components. So families that use every label are exactly sets of connected families on disjoint labeled blocks.
Define
$$A(x)=\sum_{n\ge 0} a_n x^n,\qquad H(x)=\sum_{n\ge 1} h_n x^n,$$
where \(h_n\) is the number of connected families using all \(n\) labels, divided by \(n!\). The exponential formula gives
$$A(x)=\exp(H(x)).$$
Step 4: Recover the connected coefficients
Differentiating \(A(x)=\exp(H(x))\) gives
$$A'(x)=H'(x)A(x).$$
Matching the coefficient of \(x^{m-1}\) yields
$$m a_m=\sum_{i=1}^{m} i h_i a_{m-i}.$$
Because \(a_0=1\), this becomes the recurrence
$$h_m=a_m-\frac{1}{m}\sum_{i=1}^{m-1} i h_i a_{m-i}.$$
Thus the connected counts can be recovered one size at a time after the all-label counts are known.
Step 5: Allow unused labels and require exactly \(K\) components
The original problem does not force every label to appear in a chosen subset. An unused label contributes no graph vertex and can be chosen independently, which gives the labeled factor
$$e^x=\sum_{j\ge 0}\frac{x^j}{j!}.$$
A set of exactly \(K\) connected components contributes
$$\frac{H(x)^K}{K!}.$$
Therefore the required count is
$$\boxed{C(n,K)=n!\,\left[x^n\right]\left(e^x\frac{H(x)^K}{K!}\right).}$$
This is the coefficient formula evaluated by the implementations.
Worked Example: \(n=2\)
We have \(R_0=0\), \(R_1=1\), and \(R_2=2^3-1=7\). Hence
$$U_1=1,\qquad U_2=7-2\cdot 1=5,$$
so
$$a_1=1,\qquad a_2=\frac{5}{2}.$$
The recurrence gives
$$h_1=a_1=1,\qquad h_2=a_2-\frac{1}{2}h_1a_1=\frac{5}{2}-\frac{1}{2}=2.$$
Thus
$$H(x)=x+2x^2+O(x^3).$$
For one connected component,
$$C(2,1)=2!\,\left[x^2\right]\left(e^xH(x)\right)=2!\,(2+1)=6.$$
For two connected components,
$$C(2,2)=2!\,\left[x^2\right]\left(e^x\frac{H(x)^2}{2}\right)=2!\cdot\frac{1}{2}=1.$$
This matches direct inspection: among the seven non-empty families on \(\{1,2\}\), only the family containing \(\{1\}\) and \(\{2\}\) but not \(\{1,2\}\) has two connected components.
How the Code Works
The C++, Python, and Java implementations all work modulo \(10^9+7\). They precompute factorials, inverse factorials, and modular inverses so that binomial factors and exponential-generating-function coefficients can be handled with simple modular multiplications.
Next they evaluate the raw counts \(R_r=2^{2^r-1}-1\). Because the modulus is prime, the enormous exponent \(2^r-1\) is reduced modulo \(10^9+6\) before fast modular exponentiation is applied. This keeps the computation small while preserving the correct modular value.
The inclusion-exclusion step is implemented as one polynomial convolution of the signed sequence \(((-1)^r R_r/r!)_{r\ge 0}\) with \((1/r!)_{r\ge 0}\). After the final sign adjustment, this produces the normalized coefficients \(a_n\).
The connected coefficients \(h_n\) are then recovered sequentially from the recurrence above. To obtain exactly \(K\) components, the implementation raises the truncated series \(H(x)\) to the \(K\)-th power by binary exponentiation, multiplies by the truncated series for \(e^x\), extracts the coefficient of \(x^N\), and finally multiplies by \(N!/K!\).
The C++ implementation additionally parallelizes the expensive convolutions and includes optional checkpoint validations on small cases before tackling the large target. The Python and Java implementations use the same mathematics with straightforward quadratic convolutions.
Complexity Analysis
Let \(M=\max(n,K)\). Precomputing factorials, modular inverses, and the raw counts up to \(M\) costs \(O(M\log \text{MOD})\) time and \(O(M)\) memory, which is not the dominant part.
The main cost comes from truncated polynomial convolutions of degree \(M\). With the naive multiplication used here, one convolution costs \(O(M^2)\) time and \(O(M)\) additional memory. Recovering all connected coefficients from the recurrence is also \(O(M^2)\).
Exponentiating \(H(x)\) to the \(K\)-th power requires \(O(\log K)\) convolutions, so the total running time is \(O(M^2\log K)\) with an additional \(O(M^2)\) term from the recurrence. In practice the quadratic series operations dominate. Memory usage remains \(O(M)\) because only a small number of coefficient arrays are active at once.
Footnotes and References
- Problem page: https://projecteuler.net/problem=553
- Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
- Exponential formula in combinatorics: Wikipedia - Exponential formula
- Exponential generating function: Wikipedia - Generating function
- Connected component in graph theory: Wikipedia - Connected component
Problem 553 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <queue>
#include <stdexcept>
#include <thread>
#include <vector>
namespace {
using int64 = long long;
constexpr int kMod = 1'000'000'007;
constexpr int kModMinus1 = kMod - 1;
int64 mod_pow(int64 base, int64 exp, int64 mod) {
int64 res = 1 % mod;
int64 cur = base % mod;
while (exp > 0) {
if (exp & 1LL) res = static_cast<int64>((__int128)res * cur % mod);
cur = static_cast<int64>((__int128)cur * cur % mod);
exp >>= 1LL;
}
return res;
}
int mod_inv(int a) {
return static_cast<int>(mod_pow(a, kMod - 2, kMod));
}
int mulmod(int64 a, int64 b) {
return static_cast<int>((a * b) % kMod);
}
// Brute-force validator for n <= 4.
int64 brute_C(int n, int k_wanted) {
if (n > 4) return -1;
const int qsz = (1 << n) - 1;
const int total_families = (1 << qsz) - 1;
std::vector<int> subsets;
subsets.reserve(qsz);
for (int mask = 1; mask < (1 << n); ++mask) subsets.push_back(mask);
int64 count = 0;
for (int fam_mask = 1; fam_mask <= total_families; ++fam_mask) {
std::vector<int> verts;
for (int i = 0; i < qsz; ++i) {
if (fam_mask & (1 << i)) verts.push_back(subsets[i]);
}
const int vcount = static_cast<int>(verts.size());
std::vector<char> vis(vcount, 0);
int comps = 0;
for (int i = 0; i < vcount; ++i) {
if (vis[i]) continue;
++comps;
std::queue<int> q;
q.push(i);
vis[i] = 1;
while (!q.empty()) {
int u = q.front();
q.pop();
for (int v = 0; v < vcount; ++v) {
if (vis[v]) continue;
if ((verts[u] & verts[v]) != 0) {
vis[v] = 1;
q.push(v);
}
}
}
}
if (comps == k_wanted) ++count;
}
return count;
}
struct PolyMul {
int N;
int threads;
PolyMul(int n, int t) : N(n), threads(std::max(1, t)) {}
std::vector<int> mul(const std::vector<int>& a, const std::vector<int>& b) const {
std::vector<int> c(N + 1, 0);
if (threads <= 1 || N < 800) {
for (int i = 0; i <= N; ++i) {
__int128 sum = 0;
for (int j = 0; j <= i; ++j) sum += (__int128)a[j] * b[i - j];
c[i] = static_cast<int>(sum % kMod);
}
return c;
}
const long long total_ops = 1LL * (N + 1) * (N + 2) / 2;
std::vector<int> boundary(threads + 1);
boundary[0] = -1;
boundary[threads] = N;
int cur = -1;
auto prefix_ops = [](long long i) -> long long {
if (i < 0) return 0;
return (i + 1) * (i + 2) / 2;
};
for (int t = 1; t < threads; ++t) {
long long target = total_ops * t / threads;
while (cur < N && prefix_ops(cur) < target) ++cur;
boundary[t] = cur;
}
std::vector<std::thread> pool;
pool.reserve(threads);
for (int t = 0; t < threads; ++t) {
int L = boundary[t] + 1;
int R = boundary[t + 1] + 1;
pool.emplace_back([&, L, R]() {
for (int i = L; i < R; ++i) {
__int128 sum = 0;
for (int j = 0; j <= i; ++j) sum += (__int128)a[j] * b[i - j];
c[i] = static_cast<int>(sum % kMod);
}
});
}
for (auto& th : pool) th.join();
return c;
}
};
std::vector<int> poly_pow(std::vector<int> base, int exp, const PolyMul& pm) {
const int N = pm.N;
std::vector<int> res(N + 1, 0);
res[0] = 1;
auto is_identity = [&](const std::vector<int>& p) {
if (p[0] != 1) return false;
for (int i = 1; i <= N; ++i) {
if (p[i] != 0) return false;
}
return true;
};
while (exp > 0) {
if (exp & 1) {
if (is_identity(res)) res = base;
else res = pm.mul(res, base);
}
exp >>= 1;
if (exp == 0) break;
base = pm.mul(base, base);
}
return res;
}
int solve(int N, int K, int threads, bool validate) {
int max_n = std::max(N, K);
if (validate) max_n = std::max(max_n, 100);
std::vector<int> fact(max_n + 1), invfact(max_n + 1), invInt(max_n + 1);
fact[0] = 1;
for (int i = 1; i <= max_n; ++i) fact[i] = mulmod(fact[i - 1], i);
invfact[max_n] = mod_inv(fact[max_n]);
for (int i = max_n; i >= 1; --i) invfact[i - 1] = mulmod(invfact[i], i);
invInt[1] = 1;
for (int i = 2; i <= max_n; ++i) {
invInt[i] = kMod - static_cast<int>(1LL * (kMod / i) * invInt[kMod % i] % kMod);
}
std::vector<int> pow2_modM1(max_n + 1, 1);
for (int r = 1; r <= max_n; ++r) {
pow2_modM1[r] = static_cast<int>(2LL * pow2_modM1[r - 1] % kModMinus1);
}
std::vector<int> T(max_n + 1, 0);
for (int r = 0; r <= max_n; ++r) {
int e = pow2_modM1[r] - 1;
if (e < 0) e += kModMinus1;
int64 v = mod_pow(2, e, kMod);
int val = static_cast<int>(v) - 1;
if (val < 0) val += kMod;
T[r] = val;
}
if (validate && T[0] != 0) throw std::runtime_error("T(0) != 0");
PolyMul pm(max_n, threads);
std::vector<int> invfactSeries(max_n + 1), Bsign(max_n + 1);
for (int i = 0; i <= max_n; ++i) invfactSeries[i] = invfact[i];
for (int r = 0; r <= max_n; ++r) {
int64 b = static_cast<int64>(T[r]) * invfact[r] % kMod;
int bb = static_cast<int>(b);
if (r & 1) bb = (bb == 0 ? 0 : kMod - bb);
Bsign[r] = bb;
}
std::vector<int> convA = pm.mul(Bsign, invfactSeries);
std::vector<int> Acoef(max_n + 1, 0);
Acoef[0] = 1; // empty structure for the exponential formula.
for (int n = 1; n <= max_n; ++n) {
int v = convA[n];
if (n & 1) v = (v == 0 ? 0 : kMod - v);
Acoef[n] = v;
}
std::vector<int> Gcoef(max_n + 1, 0), iG(max_n + 1, 0);
for (int m = 1; m <= max_n; ++m) {
int64 sum = 0;
for (int i = 1; i < m; ++i) {
sum += static_cast<int64>(iG[i]) * Acoef[m - i] % kMod;
if (sum >= kMod) sum -= kMod;
}
int64 corr = sum * invInt[m] % kMod;
int gm = static_cast<int>(Acoef[m] - corr);
if (gm < 0) gm += kMod;
Gcoef[m] = gm;
iG[m] = static_cast<int>(1LL * m * Gcoef[m] % kMod);
}
if (validate && Gcoef[1] != 1) throw std::runtime_error("Gcoef[1] != 1");
std::vector<int> Gpow = poly_pow(Gcoef, K, pm);
std::vector<int> S = pm.mul(invfactSeries, Gpow);
int answer = static_cast<int>(1LL * fact[N] * S[N] % kMod * invfact[K] % kMod);
if (validate) {
if (brute_C(2, 1) != 6) throw std::runtime_error("brute C(2,1) mismatch");
if (brute_C(3, 1) != 111) throw std::runtime_error("brute C(3,1) mismatch");
if (brute_C(4, 2) != 486) throw std::runtime_error("brute C(4,2) mismatch");
if (N >= 100 && K == 10) {
int chk = static_cast<int>(1LL * fact[100] * S[100] % kMod * invfact[10] % kMod);
if (chk != 728209718) throw std::runtime_error("C(100,10) mismatch");
}
std::cerr << "Validation checkpoints passed.\n";
}
return answer;
}
} // namespace
int main(int argc, char** argv) {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
int N = 10'000;
int K = 10;
unsigned hw = std::thread::hardware_concurrency();
int threads = hw ? static_cast<int>(hw) : 1;
threads = std::max(1, std::min(threads, 8));
bool validate = true;
// Optional CLI: ./a.out [N] [K] [threads] [validate(0/1)]
if (argc >= 2) N = std::stoi(argv[1]);
if (argc >= 3) K = std::stoi(argv[2]);
if (argc >= 4) threads = std::max(1, std::stoi(argv[3]));
if (argc >= 5) validate = (std::stoi(argv[4]) != 0);
try {
int answer = solve(N, K, threads, validate);
std::cout << answer << "\n";
} catch (const std::exception& e) {
std::cerr << "Validation/runtime error: " << e.what() << "\n";
return 1;
}
return 0;
}
Python
def solve():
MOD = 1_000_000_007
N_VAL = 10_000
K_VAL = 10
def mod_pow(base, exp, mod):
result = 1 % mod
cur = base % mod
while exp > 0:
if exp & 1: result = result * cur % mod
cur = cur * cur % mod
exp >>= 1
return result
def mod_inv(a): return mod_pow(a, MOD - 2, MOD)
max_n = max(N_VAL, K_VAL, 100)
fact = [1] * (max_n + 1)
for i in range(1, max_n + 1): fact[i] = fact[i-1] * i % MOD
invfact = [1] * (max_n + 1)
invfact[max_n] = mod_inv(fact[max_n])
for i in range(max_n, 0, -1): invfact[i-1] = invfact[i] * i % MOD
invInt = [0] * (max_n + 1)
invInt[1] = 1
for i in range(2, max_n + 1):
invInt[i] = MOD - (MOD // i) * invInt[MOD % i] % MOD
kModM1 = MOD - 1
pow2_modM1 = [1] * (max_n + 1)
for r in range(1, max_n + 1):
pow2_modM1[r] = 2 * pow2_modM1[r-1] % kModM1
T = [0] * (max_n + 1)
for r in range(max_n + 1):
e = pow2_modM1[r] - 1
if e < 0: e += kModM1
v = mod_pow(2, e, MOD)
T[r] = (v - 1) % MOD
def poly_mul(a, b, n):
c = [0] * (n + 1)
for i in range(n + 1):
s = 0
for j in range(i + 1):
s += a[j] * b[i - j]
c[i] = s % MOD
return c
def poly_pow(base, exp, n):
res = [0] * (n + 1)
res[0] = 1
while exp > 0:
if exp & 1: res = poly_mul(res, base, n)
exp >>= 1
if exp > 0: base = poly_mul(base, base, n)
return res
invfactS = list(invfact[:max_n+1])
Bsign = [0] * (max_n + 1)
for r in range(max_n + 1):
b = T[r] * invfact[r] % MOD
Bsign[r] = (MOD - b) % MOD if r & 1 else b
convA = poly_mul(Bsign, invfactS, max_n)
Acoef = [0] * (max_n + 1)
Acoef[0] = 1
for n in range(1, max_n + 1):
v = convA[n]
if n & 1: v = (MOD - v) % MOD
Acoef[n] = v
Gcoef = [0] * (max_n + 1)
iG = [0] * (max_n + 1)
for m in range(1, max_n + 1):
s = 0
for i in range(1, m):
s = (s + iG[i] * Acoef[m - i]) % MOD
corr = s * invInt[m] % MOD
gm = (Acoef[m] - corr) % MOD
Gcoef[m] = gm
iG[m] = m * gm % MOD
Gpow = poly_pow(Gcoef, K_VAL, max_n)
S = poly_mul(invfactS, Gpow, max_n)
answer = fact[N_VAL] * S[N_VAL] % MOD * invfact[K_VAL] % MOD
return str(answer)
if __name__ == '__main__':
print(solve())
Java
public class Euler553 {
static final int kMod = 1000000007;
static final int kModMinus1 = kMod - 1;
static long modPow(long base, long exp, long mod) {
long res = 1 % mod;
long cur = base % mod;
while (exp > 0) {
if ((exp & 1) != 0) {
res = mulMod(res, cur, mod);
}
cur = mulMod(cur, cur, mod);
exp >>= 1;
}
return res;
}
static long mulMod(long a, long b, long mod) {
long res = (a * b - (long) ((double) a * b / mod) * mod) % mod;
if (res < 0)
res += mod;
return res;
}
static int modInv(int a) {
return (int) modPow(a, kMod - 2, kMod);
}
static int mulModInt(long a, long b) {
return (int) ((a * b) % kMod);
}
static class PolyMul {
int N;
PolyMul(int n) {
this.N = n;
}
int[] mul(int[] a, int[] b) {
int[] c = new int[N + 1];
for (int i = 0; i <= N; i++) {
long sum = 0;
for (int j = 0; j <= i; j++) {
sum += (long) a[j] * b[i - j];
if ((j & 7) == 7)
sum %= kMod;
}
c[i] = (int) (sum % kMod);
}
return c;
}
}
static int[] polyPow(int[] base, int exp, PolyMul pm) {
int N = pm.N;
int[] res = new int[N + 1];
res[0] = 1;
while (exp > 0) {
if ((exp & 1) != 0) {
boolean isIdentity = true;
if (res[0] != 1)
isIdentity = false;
for (int i = 1; i <= N; i++) {
if (res[i] != 0) {
isIdentity = false;
break;
}
}
if (isIdentity) {
System.arraycopy(base, 0, res, 0, N + 1);
} else {
res = pm.mul(res, base);
}
}
exp >>= 1;
if (exp == 0)
break;
base = pm.mul(base, base);
}
return res;
}
static int solveNK(int N, int K) {
int maxN = Math.max(N, K);
int[] fact = new int[maxN + 1];
int[] invfact = new int[maxN + 1];
int[] invInt = new int[maxN + 1];
fact[0] = 1;
for (int i = 1; i <= maxN; ++i)
fact[i] = mulModInt(fact[i - 1], i);
invfact[maxN] = modInv(fact[maxN]);
for (int i = maxN; i >= 1; --i)
invfact[i - 1] = mulModInt(invfact[i], i);
if (maxN >= 1)
invInt[1] = 1;
for (int i = 2; i <= maxN; ++i) {
invInt[i] = kMod - (int) (1L * (kMod / i) * invInt[kMod % i] % kMod);
}
int[] pow2ModM1 = new int[maxN + 1];
pow2ModM1[0] = 1;
for (int r = 1; r <= maxN; ++r) {
pow2ModM1[r] = (int) (2L * pow2ModM1[r - 1] % kModMinus1);
}
int[] T = new int[maxN + 1];
for (int r = 0; r <= maxN; ++r) {
int e = pow2ModM1[r] - 1;
if (e < 0)
e += kModMinus1;
long v = modPow(2, e, kMod);
int val = (int) v - 1;
if (val < 0)
val += kMod;
T[r] = val;
}
PolyMul pm = new PolyMul(maxN);
int[] invfactSeries = new int[maxN + 1];
int[] Bsign = new int[maxN + 1];
for (int i = 0; i <= maxN; ++i)
invfactSeries[i] = invfact[i];
for (int r = 0; r <= maxN; ++r) {
long b = 1L * T[r] * invfact[r] % kMod;
int bb = (int) b;
if ((r & 1) != 0)
bb = (bb == 0 ? 0 : kMod - bb);
Bsign[r] = bb;
}
int[] convA = pm.mul(Bsign, invfactSeries);
int[] Acoef = new int[maxN + 1];
Acoef[0] = 1;
for (int n = 1; n <= maxN; ++n) {
int v = convA[n];
if ((n & 1) != 0)
v = (v == 0 ? 0 : kMod - v);
Acoef[n] = v;
}
int[] Gcoef = new int[maxN + 1];
int[] iG = new int[maxN + 1];
for (int m = 1; m <= maxN; ++m) {
long sum = 0;
for (int i = 1; i < m; ++i) {
sum += 1L * iG[i] * Acoef[m - i] % kMod;
if (sum >= kMod)
sum -= kMod;
}
long corr = sum * invInt[m] % kMod;
int gm = (int) (Acoef[m] - corr);
if (gm < 0)
gm += kMod;
Gcoef[m] = gm;
iG[m] = (int) (1L * m * Gcoef[m] % kMod);
}
int[] Gpow = polyPow(Gcoef, K, pm);
int[] S = pm.mul(invfactSeries, Gpow);
int answer = (int) (1L * fact[N] * S[N] % kMod * invfact[K] % kMod);
return answer;
}
public static String solve() {
return Integer.toString(solveNK(10000, 10));
}
public static void main(String[] args) {
System.out.println(solve());
}
}