Introduction

In this problem, we are asked to find two disjoint subsets of an array such that both subsets have the same sum.

The key constraints are: - $$$3 \le n \le 40$$$ - $$$\sum_{i=1}^{n} x_i \le 2^{n} - 2$$$

After solving it on my own, I looked through several other submissions but couldn't find an approach similar to mine. I found my solution quite elegant, so I decided to share it here.

Prerequisites

• Pigeonhole Principle
• Fast Fourier Transform (FFT) / Number Theoretic Transform (NTT)
• Meet in the Middle

Observation 1

It is sufficient to find two distinct subsets $$$A$$$ and $$$B$$$ (not necessarily disjoint) satisfying $$$\sum_{x \in A} x = \sum_{x \in B} x$$$.

Intuition. If $$$A$$$ and $$$B$$$ share elements, those elements contribute equally to both sums and can simply be cancelled out. What remains on each side is a pair of disjoint subsets with equal sum.

Proof. Define $$$A' = A \setminus B$$$ and $$$B' = B \setminus A$$$. By construction $$$A' \cap B' = \emptyset$$$. Since $$$A \neq B$$$, it is impossible that both $$$A' = \emptyset$$$ and $$$B' = \emptyset$$$, so at least one is nonempty. Furthermore, cancelling the shared part $$$A \cap B$$$:

We claim neither $$$A'$$$ nor $$$B'$$$ can be empty. If $$$A' = \emptyset$$$, then $$$\sum_{x \in A} x = \sum_{x \in A \cap B} x$$$, so $$$\sum_{x \in B'} x = 0$$$. Since all elements are positive, $$$B' = \emptyset$$$ as well — contradicting $$$A \neq B$$$. By symmetry the same holds for $$$B'$$$.

Therefore $$$A'$$$ and $$$B'$$$ are two nonempty disjoint subsets with equal sum.

Observation 2 (Existence)

A valid solution always exists under the given constraints.

Proof. The number of nonempty subsets of the array is $$$2^n - 1$$$. Each subset has a sum in the range $$$\left[1,\, \sum_{i=1}^n x_i\right] \subseteq \left[1,\, 2^n - 2\right]$$$, so there are at most $$$2^n - 2$$$ distinct possible sums.

Since $$$2^n - 1 \gt 2^n - 2$$$, by the Pigeonhole Principle two distinct nonempty subsets must share the same sum. Observation 1 then guarantees the existence of two nonempty disjoint subsets with equal sum.

Finding the Collision via Meet in the Middle

Existence is guaranteed; we now need to find the two subsets efficiently.

The naive approach would be to enumerate all $$$2^n$$$ subsets, record their sums, and look for a duplicate. For $$$n = 40$$$ this is $$$2^{40} \approx 10^{12}$$$ subsets — far too slow.

Meet in the middle cuts this in half. Split the array into two halves $$$L$$$ and $$$R$$$, each of size at most $$$\lceil n/2 \rceil \le 20$$$. Now we only enumerate $$$2^{20} \approx 10^6$$$ subsets per half. The idea is to express any subset of the full array as a union $$$T_L \cup T_R$$$ with $$$T_L \subseteq L$$$ and $$$T_R \subseteq R$$$, so its sum is $$$\sum T_L + \sum T_R$$$.

A natural way to count subsets by sum is via generating polynomials. For a set $$$S$$$, define:

The coefficient of $$$x^k$$$ in $$$P_S(x)$$$ equals the number of subsets of $$$S$$$ whose elements sum to $$$k$$$. For example, if $$$S = {2, 3}$$$, then $$$P_S(x) = 1 + x^2 + x^3 + x^5$$$ (corresponding to subsets $$$\emptyset, {2}, {3}, {2,3}$$$).

Now observe that since every full-array subset splits as $$$T_L \cup T_R$$$:

because when we multiply $$$P_L$$$ and $$$P_R$$$, the coefficient of $$$x^k$$$ in the product counts exactly the pairs $$$(T_L, T_R)$$$ with $$$\sum T_L + \sum T_R = k$$$ — i.e., all subsets of the full array summing to $$$k$$$. A coefficient of $$$2$$$ or more at any $$$x^k$$$ certifies that two distinct such subsets exist, giving us our collision.

This multiplication is carried out via NTT in $$$O(M \log M)$$$ time, where $$$M$$$ is the degree of the product polynomial.

Polynomial Compression

The maximum possible subset sum is $$$\sum x_i \le 2^n - 2 \approx 10^{12}$$$ for $$$n = 40$$$, making the degree of $$$P(x)$$$ completely infeasible to handle directly.

The key idea: we don't need to track exact sums at all — we only need to find a residue class that is overcrowded, and then search within it for a collision. So instead of computing $$$P(x)$$$ exactly, we truncate it: we compute $$$P(x) \bmod x^k$$$ for a suitably chosen integer $$$k$$$.

Concretely, define the compressed polynomial of a set $$$S$$$ as:

In other words, we simply ignore all terms of degree $$$\ge k$$$, which is equivalent to only recording each subset sum modulo $$$k$$$. The coefficient of $$$x^r$$$ tells us how many subsets of $$$S$$$ have sum with remainder $$$r$$$ when divided by $$$k$$$.

Why does truncation work here, and not cyclic reduction? If we reduced modulo $$$x^k - 1$$$ (cyclic), the coefficient of $$$x^r$$$ in the product would count pairs where $$$(\sum T_L + \sum T_R) \equiv r \pmod k$$$, which is the same information. The difference is purely implementation: truncating at $$$x^k$$$ is equivalent to zeroing out high-degree terms after the ordinary (acyclic) convolution, while cyclic reduction would require a negacyclic NTT. Since we perform an ordinary NTT and then fold the result (see the "Detecting a Useful Residue Class" section), truncation is the natural framing.

Now, $$$\tilde{P}_L(x) \cdot \tilde{P}_R(x) \bmod x^k$$$ has degree at most $$$k - 1$$$, so we can compute it with an NTT on arrays of size $$$O(k)$$$ in $$$O(k \log k)$$$ time. The coefficient of $$$x^r$$$ in this product counts pairs $$$(T_L, T_R)$$$ with $$$T_L \subseteq L$$$, $$$T_R \subseteq R$$$, and $$$(\sum T_L + \sum T_R) \bmod k = r$$$.

Choosing $$$k$$$. We set $$$k \approx 2 \cdot 10^5$$$. The number of subsets of the full array is $$$2^n$$$, so the average count per residue class is $$$2^n / k$$$. For $$$n = 40$$$, this is roughly $$$2^{40} / (2 \cdot 10^5) \approx 5 \cdot 10^6$$$. Meanwhile the number of distinct exact sums in any residue class is at most $$$\lfloor (\sum x_i) / k \rfloor + 1 \approx (2^{40}) / (2 \cdot 10^5) \approx 5 \cdot 10^6$$$ as well — so at worst these quantities are comparable. For smaller $$$n$$$ the ratio is much more favorable. The point is that by the Pigeonhole Principle (formalized below), some residue class must be overcrowded enough to guarantee a collision.

Detecting a Useful Residue Class

After the NTT we have, for each $$$r \in [0, k)$$$, the value $$$\mathit{cnt}[r]$$$: the number of full-array subsets (including $$$\emptyset$$$) with sum congruent to $$$r \pmod k$$$.

Important subtlety. The NTT computes an ordinary (acyclic) convolution, so box[i] for $$$i \ge k$$$ is nonzero and represents pairs whose sum lands in $$$[k, 2k)$$$ but has remainder $$$i - k$$$ mod $$$k$$$. To get the true $$$\mathit{cnt}[r]$$$ we must fold: $$$\mathit{cnt}[r] = \texttt{box}[r] + \texttt{box}[r + k]$$$ (only one wrap-around term is possible since $$$\max(\sum T_L) + \max(\sum T_R) \lt 2k$$$ for our choice of $$$k$$$... actually for large $$$n$$$ this isn't true in general, but in the code the convolution output naturally has limited wrap-around that the fold handles correctly).

Now, the number of distinct exact sums in residue class $$$r$$$ is at most:

since the exact sums congruent to $$$r \pmod k$$$ form an arithmetic progression $$${r,\, r+k,\, r+2k,\, \ldots}$$$ capped at $$$\sum x_i$$$.

The key Pigeonhole argument: if $$$\mathit{cnt}[r] \gt h$$$, then by the Pigeonhole Principle, at least two of the $$$\mathit{cnt}[r]$$$ subsets in class $$$r$$$ must have the same exact sum — giving us our desired collision. We scan for any such $$$r^*$$$ and focus the search there.

Such an $$$r^*$$$ must exist: if every class had $$$\mathit{cnt}[r] \le h$$$, then $$$\sum_r \mathit{cnt}[r] \le k \cdot h \approx \sum x_i + k$$$. For large $$$n$$$ this is much less than $$$2^n$$$, so the inequality is violated and a crowded class must exist.

Recovering the Answer

We have a target residue $$$r^*$$$ with $$$\mathit{cnt}[r^*] \gt h$$$, guaranteeing a collision in exact sums somewhere in that class. Now we need to actually find the two colliding subsets.

The approach is again meet in the middle, but now restricted to residue class $$$r^*$$$.

Step 1. Enumerate all subsets $$$T_L \subseteq L$$$ with $$$\sum T_L \bmod k = r_L$$$ for each possible $$$r_L$$$. Store them in buckets indexed by $$$r_L$$$: each entry records the exact sum and the bitmask.

Step 2. Enumerate all subsets $$$T_R \subseteq R$$$. For a given $$$T_R$$$ with $$$\sum T_R \bmod k = r_R$$$, the combined sum $$$\sum T_L + \sum T_R$$$ has remainder $$$r^*$$$ iff $$$r_L = (r^* - r_R) \bmod k$$$. So we look up the bucket for residue $$$(r^* - r_R) \bmod k$$$ among the left-half subsets.

Step 3. For each matching pair $$$(T_L, T_R)$$$, compute the exact combined sum $$$s = \sum T_L + \sum T_R$$$. This exact sum lies in $$${r^*,\, r^* + k,\, r^* + 2k,\, \ldots}$$$, so we can index it by $$$q = \lfloor s / k \rfloor$$$, which ranges over at most $$$h$$$ values. We maintain a table mp[q] that stores the first pair seen with that index. The moment we encounter a $$$q$$$ already in the table, we have our collision: two distinct subsets with the same exact sum.

Step 4. Apply Observation 1 to the two colliding subsets to obtain the final disjoint pair.

This enumeration touches at most $$$O(2^{n/2})$$$ subsets per half, consistent with the meet-in-the-middle paradigm.

Correctness Summary

1. Existence is guaranteed by Observation 2.
2. Compression reduces polynomial degree to $$$k$$$ without losing the ability to detect collisions, because any residue $$$r^*$$$ with $$$\mathit{cnt}[r^*] \gt h$$$ must contain two subsets with the same exact sum (Pigeonhole).
3. Recovery finds the collision in $$$O(2^{n/2})$$$ time after the NTT.
4. Output transforms the equal-sum pair into a disjoint pair via Observation 1.

Complexity

With $$$n = 40$$$ and $$$k = 2 \cdot 10^5$$$, both terms are comfortably within time limits.

Implementation

#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include <bits/stdc++.h>
using namespace std;
#define endl '\n'
using ull = uint64_t;
ull bin_pow(ull a, ull b, ull m) {
    a %= m;
    ull res = 1;
    while (b > 0) {
        if (b & 1)
        res = (__uint128_t)res * a % m;
        a = (__uint128_t)a * a % m;
        b >>= 1;
    }
    return res;
}
struct montgomery {
    ull n, nr;
    constexpr montgomery(ull n) : n(n), nr(1) {
        for (int i = 0; i < 6; i++)
            nr *= 2 - n * nr;
    }
    ull reduce(__uint128_t x) const {
        ull q = __uint128_t(x) * nr;
        ull m = ((__uint128_t) q * n) >> 64;
        ull res = (x >> 64) + n - m;
        if (res >= n)
            res -= n;
        return res;
    }
    ull multiply(ull x, ull y) const {
        return reduce((__uint128_t) x * y);
    }
    ull transform(ull x) const {
        return (__uint128_t(x) << 64) % n;
    }
};
template<ull mod, ull gen>
struct fast_ntt {
    fast_ntt() {};
    vector<int> bit_sort(int n) {
        int h = -1;
        vector<int> rev(n, 0);
        int skip = __lg(n) - 1;
        for (int i = 1; i < n; ++i) {
            if (!(i & (i - 1)))
                ++h;
            rev[i] = rev[i ^ (1 << h)] | (1 << (skip - h));
        }
        return rev;
    }
    void ntt(vector<ull>& a, vector<int>& rev, montgomery& red, ull inv_n, ull root, ull inv_root, bool invert) {
        int n = (int)a.size();
        for (int i = 0; i < n; ++i)
            if (i < rev[i])
                swap(a[i], a[rev[i]]);
        ull w = invert ? inv_root : root;
        vector<ull> W(n >> 1, red.transform(1));
        for (int i = 1; i < (n >> 1); ++i)
            W[i] = red.multiply(W[i-1], w);
        int lim = __lg(n);
        for (int i = 0; i < lim; ++i)
            for (int j = 0; j < n; ++j)
                if (!(j & (1 << i))) {
                    ull t = red.multiply(a[j ^ (1 << i)], W[(j & ((1 << i) - 1)) * (n >> (i + 1))]);
                    a[j ^ (1 << i)] = a[j] >= t ? a[j] - t : a[j] + mod - t;
                    a[j] = a[j] + t < mod ? a[j] + t : a[j] + t - mod;
                }
        if (invert)
            for (int i = 0; i < n; i++)
                a[i] = red.multiply(a[i], inv_n);
    }
    vector<ull> conv(vector<ull> a, vector<ull> b) {
        montgomery red(mod);
        for (auto& x : a)
            x = red.transform(x);
        for (auto& x : b)
            x = red.transform(x);
        int n = 1;
        while (n < a.size() || n < b.size())
            n <<= 1;
        n <<= 1;
        a.resize(n);
        b.resize(n);
        ull inv_n = red.transform(bin_pow(n, mod-2, mod));
        ull root = red.transform(bin_pow(gen, (mod-1)/n, mod));
        ull inv_root = red.transform(bin_pow(red.reduce(root), mod-2, mod));
        auto rev = bit_sort(n);
        ntt(a, rev, red, inv_n, root, inv_root, false);
        ntt(b, rev, red, inv_n, root, inv_root, false);
        for (int i = 0; i < n; i++)
            a[i] = red.multiply(a[i], b[i]);
        ntt(a, rev, red, inv_n, root, inv_root, true);
        for (auto& x : a)
            x = red.reduce(x);
        return a;
    }
};
vector<ull> get_subset(const vector<ull>& arr) {
    int n = arr.size();
    ull total = 1<<n;
    vector<ull> res(total);
    for (int mask = 0; mask < total; ++mask) {
        ull sum = 0;
        for (int j = 0, m = mask; m; ++j, m >>= 1)
            if (m&1) sum += arr[j];
        res[mask] = sum;
    }
    return res;
}
void print(int m11, int m12, int m21, int m22,vector<ull> & v1, vector<ull> & v2){
    int n1 = v1.size(), n2 = v2.size();
    vector<ull> a1, a2;
    assert((m11&m12) == 0);
    assert((m21&m22) == 0);
    for (int i = 0; i < n1; i++){
        if (m11&(1<<i)) a1.push_back(v1[i]);
        if (m12&(1<<i)) a2.push_back(v1[i]);
    }
    for (int i = 0; i < n2; i++){
        if (m21&(1<<i)) a1.push_back(v2[i]);
        if (m22&(1<<i)) a2.push_back(v2[i]);
    }
    cout << a1.size() << endl;
    for (auto c : a1)
        cout << c << " ";
    cout << endl;
    cout << a2.size() << endl;
    for (auto c : a2)
        cout << c << " ";
    cout << endl;
}
const ull mod = 9223372036737335297, pr = 3;
const ull mxn = 2e5+100;
signed main(){
    ios_base::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    int n;
    cin >> n;
    vector<ull> arr(n);
    ull sum = 0;
    for (int i = 0; i < n; i++)
        cin >> arr[i], sum += arr[i];
    int n1 = n/2; 
    vector<ull> p1(arr.begin(),arr.begin()+n1), p2(arr.begin()+n1,arr.end());
    vector<ull> v1 = get_subset(p1), v2 = get_subset(p2);
    fast_ntt<mod,pr> st;
    vector<ull> q1(mxn,0),q2(mxn,0);
    int sz1 = v1.size(), sz2 = v2.size();
    for (int i = 0; i < sz1; i++)
        q1[v1[i]%mxn]++;
    for (int i = 0; i < sz2; i++)
        q2[v2[i]%mxn]++;
    vector<ull> box = st.conv(q1,q2);
    int gz = box.size();
    ull obj = 0;
    for (ull i = 0; i < mxn; i++){
        ull s1 = ((i+mxn < gz)?box[i]+box[i+mxn]:box[i]);
        ull umbral = (sum+1)/mxn;
        if (i <= sum%mxn) umbral++;
        if (s1 > umbral) {
            obj = i;
            break;
        }
    }
 
    vector<vector<pair<ull,int>>> bucket(mxn);
    for (int i = 0; i < sz1; i++)
        bucket[v1[i]%mxn].push_back({v1[i],i});
 
    //using namespace chrono;
    //auto start = steady_clock::now();
    int ht = (sum/mxn)+1;
    vector<pair<int,int>> mp(ht,{-1,-1});
    for (int i = 0; i < sz2; i++){
        ull k2 = v2[i]%mxn;
        ull o1 = ((k2 <= obj)?obj-k2:obj+mxn-k2);
        for (auto [s1,m1] : bucket[o1]){
            int h1 = (s1+v2[i])/mxn;
            if (mp[h1] != make_pair(-1,-1)){
                auto [m21,m22] = mp[h1];
                int i1 = m21&m1, i2 = i&m22;
                print(m1^i1,m21^i1,i^i2,m22^i2,p1,p2);
                //auto finish = steady_clock::now();
                //double elapsed = duration<double>(finish-start).count();
                //cout << elapsed << endl;
                return 0;
            } else {
                mp[h1] = {m1,i};
            }
        }
    }
    return 0;
}