<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

247. Difficult Choice
time limit per test: 0.25 sec.
memory limit per test: 65536 KB
input: standard
output: standard



In a bowling club "Odd Prime Ball" 2P colored bowling balls numbered from 1 to 2P are available for the game, where P is an odd prime number.
Antony noticed that he has a good fortune if he correctly chooses balls for the game. He thinks that his choice is correct if the sum of numbers of chosen balls is divisible by P. Moreover, the number of balls Antony choses must be equal to P, otherwise he will gain no fortune and lose the game. So he needs to choose P different balls for the game.
He wants to know how many ways he can choose the balls for the game to have a good fortune.
In the first example from balls with weights (1, 2, 3, 4, 5, 6) Antony can choose the following 8 combinations: (1, 2, 3); (1, 2, 6); (1, 3, 5); (1, 5, 6); (2, 3, 4); (2, 4, 6); (3, 4, 5); (4, 5, 6).

Input
The first line of input contains single integer N the number of tests (0<=N<=100). The following N lines contain tests descriptions (one line per test). Each test description consists of a single odd prime number P (1<P<1000).

Output
Output N lines (one line per test) with requested numbers.

Sample test(s)

Input
2
3
5

Output
8
52

<|response|>
## 1) Abridged problem statement (concise)

For each test case, given an odd prime \(P\) \((1<P<1000)\), consider numbers \(1,2,\dots,2P\). Count how many ways to choose **exactly \(P\) distinct numbers** such that their sum is divisible by \(P\). Output the count for each test.

---

## 2) Key observations needed to solve the problem

1. **Two obvious valid subsets exist**
   - \(A=\{1,2,\dots,P\}\) and \(B=\{P+1,P+2,\dots,2P\}\) both have size \(P\).
   - \(\sum_{i=1}^P i=\frac{P(P+1)}{2}\) is divisible by \(P\), so \(A\) is valid.
   - \(\sum_{i=P+1}^{2P} i\) is also divisible by \(P\), so \(B\) is valid.
   - So we already have **2** good subsets.

2. **All other \(P\)-subsets can be partitioned into groups of size \(P\)**
   - For any other subset \(S\) (not \(A\) or \(B\)), let
     \[
     k = |S \cap \{1,\dots,P\}|
     \]
     Then \(1 \le k \le P-1\).

3. **Cyclic shift trick**
   - “Shift” only the elements of \(S\) that lie in \(\{1,\dots,P\}\) by \(i\) modulo \(P\), for \(i=0,1,\dots,P-1\).
   - This creates exactly \(P\) distinct subsets \(S_0,S_1,\dots,S_{P-1}\) (an orbit).

4. **Exactly one subset per group has sum divisible by \(P\)**
   - Shifting increases the total sum by \(k\cdot i \pmod P\).
   - Since \(P\) is prime and \(1\le k\le P-1\), \(\gcd(k,P)=1\), so \(k i \bmod P\) runs through all residues exactly once.
   - Hence among the \(P\) sums, exactly **one** is \(0 \bmod P\).

---

## 3) Full solution approach

Let total number of \(P\)-subsets be:
\[
\binom{2P}{P}.
\]

- Two of these subsets are special: \(A\) and \(B\), both good.
- The remaining \(\binom{2P}{P}-2\) subsets are partitioned into disjoint groups (orbits) of size \(P\).
- Each such group contributes exactly 1 good subset.

Therefore, total good subsets:
\[
\boxed{\;2 + \frac{\binom{2P}{P}-2}{P}\;}
\]

### Computation
- For \(P\) up to 999, \(\binom{2P}{P}\) is enormous, so we need big integers.
- Compute \(\binom{2P}{P}\) using the multiplicative exact formula:
\[
\binom{2P}{P}=\prod_{i=0}^{P-1}\frac{2P-i}{i+1}
\]
Perform multiplication then exact division at each step.

Finally compute:
\[
\text{ans} = 2 + \frac{C-2}{P}.
\]

Time per test is \(O(P)\) big-integer operations, fine for \(P<1000\).

---

## 4) C++ implementation with detailed comments

```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

// base and base_digits must be consistent
const int base = 1000000000;
const int base_digits = 9;

struct bigint {
    vector<int> z;
    int sign;

    bigint() : sign(1) {}

    bigint(long long v) { *this = v; }

    bigint(const string& s) { read(s); }

    void operator=(const bigint& v) {
        sign = v.sign;
        z = v.z;
    }

    void operator=(long long v) {
        sign = 1;
        if(v < 0) {
            sign = -1, v = -v;
        }
        z.clear();
        for(; v > 0; v = v / base) {
            z.push_back(v % base);
        }
    }

    bigint operator+(const bigint& v) const {
        if(sign == v.sign) {
            bigint res = v;

            for(int i = 0, carry = 0;
                i < (int)max(z.size(), v.z.size()) || carry; ++i) {
                if(i == (int)res.z.size()) {
                    res.z.push_back(0);
                }
                res.z[i] += carry + (i < (int)z.size() ? z[i] : 0);
                carry = res.z[i] >= base;
                if(carry) {
                    res.z[i] -= base;
                }
            }
            return res;
        }
        return *this - (-v);
    }

    bigint operator-(const bigint& v) const {
        if(sign == v.sign) {
            if(abs() >= v.abs()) {
                bigint res = *this;
                for(int i = 0, carry = 0; i < (int)v.z.size() || carry; ++i) {
                    res.z[i] -= carry + (i < (int)v.z.size() ? v.z[i] : 0);
                    carry = res.z[i] < 0;
                    if(carry) {
                        res.z[i] += base;
                    }
                }
                res.trim();
                return res;
            }
            return -(v - *this);
        }
        return *this + (-v);
    }

    void operator*=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = 0, carry = 0; i < (int)z.size() || carry; ++i) {
            if(i == (int)z.size()) {
                z.push_back(0);
            }
            long long cur = z[i] * (long long)v + carry;
            carry = (int)(cur / base);
            z[i] = (int)(cur % base);
            // asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur),
            // "c"(base));
        }
        trim();
    }

    bigint operator*(int v) const {
        bigint res = *this;
        res *= v;
        return res;
    }

    friend pair<bigint, bigint> divmod(const bigint& a1, const bigint& b1) {
        int norm = base / (b1.z.back() + 1);
        bigint a = a1.abs() * norm;
        bigint b = b1.abs() * norm;
        bigint q, r;
        q.z.resize(a.z.size());

        for(int i = a.z.size() - 1; i >= 0; i--) {
            r *= base;
            r += a.z[i];
            int s1 = b.z.size() < r.z.size() ? r.z[b.z.size()] : 0;
            int s2 = b.z.size() - 1 < r.z.size() ? r.z[b.z.size() - 1] : 0;
            int d = ((long long)s1 * base + s2) / b.z.back();
            r -= b * d;
            while(r < 0) {
                r += b, --d;
            }
            q.z[i] = d;
        }

        q.sign = a1.sign * b1.sign;
        r.sign = a1.sign;
        q.trim();
        r.trim();
        return make_pair(q, r / norm);
    }

    friend bigint sqrt(const bigint& a1) {
        bigint a = a1;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        int n = a.z.size();

        int firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int norm = base / (firstDigit + 1);
        a *= norm;
        a *= norm;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        bigint r = (long long)a.z[n - 1] * base + a.z[n - 2];
        firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int q = firstDigit;
        bigint res;

        for(int j = n / 2 - 1; j >= 0; j--) {
            for(;; --q) {
                bigint r1 =
                    (r - (res * 2 * base + q) * q) * base * base +
                    (j > 0 ? (long long)a.z[2 * j - 1] * base + a.z[2 * j - 2]
                           : 0);
                if(r1 >= 0) {
                    r = r1;
                    break;
                }
            }
            res *= base;
            res += q;

            if(j > 0) {
                int d1 =
                    res.z.size() + 2 < r.z.size() ? r.z[res.z.size() + 2] : 0;
                int d2 =
                    res.z.size() + 1 < r.z.size() ? r.z[res.z.size() + 1] : 0;
                int d3 = res.z.size() < r.z.size() ? r.z[res.z.size()] : 0;
                q = ((long long)d1 * base * base + (long long)d2 * base + d3) /
                    (firstDigit * 2);
            }
        }

        res.trim();
        return res / norm;
    }

    bigint operator/(const bigint& v) const { return divmod(*this, v).first; }

    bigint operator%(const bigint& v) const { return divmod(*this, v).second; }

    void operator/=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = (int)z.size() - 1, rem = 0; i >= 0; --i) {
            long long cur = z[i] + rem * (long long)base;
            z[i] = (int)(cur / v);
            rem = (int)(cur % v);
        }
        trim();
    }

    bigint operator/(int v) const {
        bigint res = *this;
        res /= v;
        return res;
    }

    int operator%(int v) const {
        if(v < 0) {
            v = -v;
        }
        int m = 0;
        for(int i = z.size() - 1; i >= 0; --i) {
            m = (z[i] + m * (long long)base) % v;
        }
        return m * sign;
    }

    void operator+=(const bigint& v) { *this = *this + v; }
    void operator-=(const bigint& v) { *this = *this - v; }
    void operator*=(const bigint& v) { *this = *this * v; }
    void operator/=(const bigint& v) { *this = *this / v; }

    bool operator<(const bigint& v) const {
        if(sign != v.sign) {
            return sign < v.sign;
        }
        if(z.size() != v.z.size()) {
            return z.size() * sign < v.z.size() * v.sign;
        }
        for(int i = z.size() - 1; i >= 0; i--) {
            if(z[i] != v.z[i]) {
                return z[i] * sign < v.z[i] * sign;
            }
        }
        return false;
    }

    bool operator>(const bigint& v) const { return v < *this; }
    bool operator<=(const bigint& v) const { return !(v < *this); }
    bool operator>=(const bigint& v) const { return !(*this < v); }
    bool operator==(const bigint& v) const {
        return !(*this < v) && !(v < *this);
    }
    bool operator!=(const bigint& v) const { return *this < v || v < *this; }

    void trim() {
        while(!z.empty() && z.back() == 0) {
            z.pop_back();
        }
        if(z.empty()) {
            sign = 1;
        }
    }

    bool isZero() const { return z.empty() || (z.size() == 1 && !z[0]); }

    bigint operator-() const {
        bigint res = *this;
        res.sign = -sign;
        return res;
    }

    bigint abs() const {
        bigint res = *this;
        res.sign *= res.sign;
        return res;
    }

    long long longValue() const {
        long long res = 0;
        for(int i = z.size() - 1; i >= 0; i--) {
            res = res * base + z[i];
        }
        return res * sign;
    }

    friend bigint gcd(const bigint& a, const bigint& b) {
        return b.isZero() ? a : gcd(b, a % b);
    }
    friend bigint lcm(const bigint& a, const bigint& b) {
        return a / gcd(a, b) * b;
    }

    void read(const string& s) {
        sign = 1;
        z.clear();
        int pos = 0;
        while(pos < (int)s.size() && (s[pos] == '-' || s[pos] == '+')) {
            if(s[pos] == '-') {
                sign = -sign;
            }
            ++pos;
        }
        for(int i = s.size() - 1; i >= pos; i -= base_digits) {
            int x = 0;
            for(int j = max(pos, i - base_digits + 1); j <= i; j++) {
                x = x * 10 + s[j] - '0';
            }
            z.push_back(x);
        }
        trim();
    }

    friend istream& operator>>(istream& stream, bigint& v) {
        string s;
        stream >> s;
        v.read(s);
        return stream;
    }

    friend ostream& operator<<(ostream& stream, const bigint& v) {
        if(v.sign == -1) {
            stream << '-';
        }
        stream << (v.z.empty() ? 0 : v.z.back());
        for(int i = (int)v.z.size() - 2; i >= 0; --i) {
            stream << setw(base_digits) << setfill('0') << v.z[i];
        }
        return stream;
    }

    static vector<int> convert_base(
        const vector<int>& a, int old_digits, int new_digits
    ) {
        vector<long long> p(max(old_digits, new_digits) + 1);
        p[0] = 1;
        for(int i = 1; i < (int)p.size(); i++) {
            p[i] = p[i - 1] * 10;
        }
        vector<int> res;
        long long cur = 0;
        int cur_digits = 0;
        for(int i = 0; i < (int)a.size(); i++) {
            cur += a[i] * p[cur_digits];
            cur_digits += old_digits;
            while(cur_digits >= new_digits) {
                res.push_back(int(cur % p[new_digits]));
                cur /= p[new_digits];
                cur_digits -= new_digits;
            }
        }
        res.push_back((int)cur);
        while(!res.empty() && res.back() == 0) {
            res.pop_back();
        }
        return res;
    }

    typedef vector<long long> vll;

    static vll karatsubaMultiply(const vll& a, const vll& b) {
        int n = a.size();
        vll res(n + n);
        if(n <= 32) {
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n; j++) {
                    res[i + j] += a[i] * b[j];
                }
            }
            return res;
        }

        int k = n >> 1;
        vll a1(a.begin(), a.begin() + k);
        vll a2(a.begin() + k, a.end());
        vll b1(b.begin(), b.begin() + k);
        vll b2(b.begin() + k, b.end());

        vll a1b1 = karatsubaMultiply(a1, b1);
        vll a2b2 = karatsubaMultiply(a2, b2);

        for(int i = 0; i < k; i++) {
            a2[i] += a1[i];
        }
        for(int i = 0; i < k; i++) {
            b2[i] += b1[i];
        }

        vll r = karatsubaMultiply(a2, b2);
        for(int i = 0; i < (int)a1b1.size(); i++) {
            r[i] -= a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            r[i] -= a2b2[i];
        }

        for(int i = 0; i < (int)r.size(); i++) {
            res[i + k] += r[i];
        }
        for(int i = 0; i < (int)a1b1.size(); i++) {
            res[i] += a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            res[i + n] += a2b2[i];
        }
        return res;
    }

    bigint operator*(const bigint& v) const {
        vector<int> a6 = convert_base(this->z, base_digits, 6);
        vector<int> b6 = convert_base(v.z, base_digits, 6);
        vll a(a6.begin(), a6.end());
        vll b(b6.begin(), b6.end());
        while(a.size() < b.size()) {
            a.push_back(0);
        }
        while(b.size() < a.size()) {
            b.push_back(0);
        }
        while(a.size() & (a.size() - 1)) {
            a.push_back(0), b.push_back(0);
        }
        vll c = karatsubaMultiply(a, b);
        bigint res;
        res.sign = sign * v.sign;
        for(int i = 0, carry = 0; i < (int)c.size(); i++) {
            long long cur = c[i] + carry;
            res.z.push_back((int)(cur % 1000000));
            carry = (int)(cur / 1000000);
        }
        res.z = convert_base(res.z, 6, base_digits);
        res.trim();
        return res;
    }
};

int p;

void read() { cin >> p; }

void solve() {
    // There is a classic DP solution but that would be too slow, especially as
    // we need big integers. It's possible to precompute all answers which might
    // be doable as there are 1000 unique inputs, but that solution isn't
    // interesting.
    //
    // Instead, because P is an odd prime we have nice properties. Let's first
    // consider the two sets {1, ..., P} and {P+1, ..., 2P} which already have
    // sums divisible by P. Now we can consider all other subsets (there are
    // C(2*P,P) - 2 of them). We can make groups of exactly P subsets each:
    //
    //   Take any such subset. Let k = how many of its numbers
    //   are from 1 to P (so 1 <= k <= P-1).
    //
    //   Now create P versions by adding the same number i (0 to P-1)
    //   to each of those k small numbers, then wrap around mod P
    //   so they stay between 1 and P. The big numbers stay the same.
    //
    //   This gives P different P-sized subsets. Their sums change
    //   by k*i each time. Since k and P are coprime, k*i mod P
    //   hits every possible remainder exactly once as i goes 0..P-1.
    //
    //   Therefore exactly one subset in each team of P has sum
    //   divisible by P.
    //
    // So the ordinary subsets contribute (C(2*P,P) - 2) / P good ones.
    // Adding the two special ones gives 2 + (C(2*P,P) - 2) / P.
    //
    // This is explanation is heavily based on the below and it turns out this
    // problem also appeared at IMO 1995:
    //
    //     https://math.stackexchange.com/questions/314788/let-p-be-an-odd-prime-
    //             number-how-many-p-element-subsets-of-1-2-3-4-ldo

    bigint res = 1;
    for(int i = 0; i < p; ++i) {
        res *= (2 * p - i);
        res /= (i + 1);
    }
    res -= 2;
    res /= p;
    res += 2;
    cout << res << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

---

## 5) Python implementation with detailed comments

```python
import sys

def solve():
    data = sys.stdin.read().strip().split()
    if not data:
        return
    t = int(data[0])
    idx = 1
    out = []

    for _ in range(t):
        p = int(data[idx]); idx += 1

        # Compute C(2p, p) exactly using multiplicative formula.
        # Python integers are arbitrary precision.
        c = 1
        for i in range(p):
            c *= (2 * p - i)
            c //= (i + 1)   # exact division

        # Apply the proven counting formula:
        # ans = 2 + (c - 2) / p
        ans = 2 + (c - 2) // p
        out.append(str(ans))

    sys.stdout.write("\n".join(out))

if __name__ == "__main__":
    solve()
```

These implementations directly encode the math result:
\[
\boxed{\text{answer} = 2 + \frac{\binom{2P}{P}-2}{P}}
\]
and handle the large values via big integers.