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

241. The United Fields of Chessboardia
time limit per test: 0.25 sec.
memory limit per test: 4096 KB
input: standard
output: standard



Little boy Vasya likes to play chess very much. Every Wednesday and Saturday at four o'clock after midday he goes to "An Check-Mate" (ACM) chess Club for five hours. Here Vasya solves many interesting chess problems and tasks.
Vasya has a great skill of that. He proves this skill everyday. Now, he can solve any simple chess problem faster than any his classmate. But he doesn't know how to solve difficult chess problems, and wants to find out a way to learn that. His trainer said, that the only way to learn how to solve difficult problem - is to know how others do that. And Vasya has thought up a problem for solving by some clever man. Please solve it and help Vasya to make an self-improvement.
Trying to create a really difficult problem Vasya decided not to use standard 8x8 chessboard.
Consider two chessboards NxN and MxM aligned to integer grid such that the left-bottom corner of the second chessboard located not to the bottom and not to the left from the left-bottom corner of the first chessboard. Let us call horizontal distance the number of columns located between two corners. Analogically, we will define the vertical distance (look at figure for clearance). Your task is, given sizes of chessboards, horizontal and vertical distances between their left-bottom corners and number K, to calculate a number of ways to place K rooks in peaceful position on union of these two boards. Position is called peaceful, if any two rooks don't attack each other. Rooks attack each other if located at same vertical or horizontal.


Input
On first line of input file, there are sizes of chessboards, N and M respectively (0<=N,M<=20), horizontal and vertical distances between left-bottom corners of chessboards, W and H (0<=W,H<=10^9) and K - number of rooks (0<=K<=10^9).

Output
First line of output file must contain only one number - quantity of peaceful positions.

Sample test(s)

Input
Test #1
8 2 6 8 1

Test #2
8 8 3 4 1

Output
Test #1
68

Test #2
108
Author:	Alexey Preobrajensky
Resource:	---
Date:	October, 2003

<|response|>
## 1. Abridged problem statement

We have two chessboards placed on the same infinite integer grid:

- first board: `N × N`, bottom-left corner at `(0, 0)`;
- second board: `M × M`, bottom-left corner at `(W, H)`.

Consider all cells belonging to the union of these two boards.

Count the number of ways to place exactly `K` rooks on these cells so that no two rooks share the same row or column.

Constraints:

```text
0 ≤ N, M ≤ 20
0 ≤ W, H ≤ 10^9
0 ≤ K ≤ 10^9
```

The answer can be very large, so it must be printed exactly.

---

## 2. Key observations needed to solve the problem

### Observation 1: Only occupied rows and columns matter

Although `W` and `H` can be huge, each board has size at most `20`.

The first board uses columns:

```text
[0, N)
```

The second board uses columns:

```text
[W, W + M)
```

So the total number of relevant columns is at most:

```text
N + M ≤ 40
```

The same is true for rows.

Large empty gaps between the boards are irrelevant.

---

### Observation 2: Split columns and rows into three types

Columns can be divided into:

1. columns belonging only to the first board;
2. columns shared by both boards;
3. columns belonging only to the second board.

Rows can be divided in the same way.

Let:

```text
c1 = number of first-only columns
cs = number of shared columns
c2 = number of second-only columns
```

and:

```text
r1 = number of first-only rows
rs = number of shared rows
r2 = number of second-only rows
```

The number of shared columns is the length of the intersection:

```text
[0, N) ∩ [W, W + M)
```

Similarly, the number of shared rows is:

```text
[0, N) ∩ [H, H + M)
```

---

### Observation 3: Allowed row types depend on column type

A first-only column contains cells only from the first board, so it can use:

- first-only rows;
- shared rows.

A second-only column contains cells only from the second board, so it can use:

- shared rows;
- second-only rows.

A shared column belongs to both boards, so it can use:

- first-only rows;
- shared rows;
- second-only rows.

So the allowed combinations are:

| Column type | first-only rows | shared rows | second-only rows |
|---|---:|---:|---:|
| first-only column | yes | yes | no |
| shared column | yes | yes | yes |
| second-only column | no | yes | yes |

---

### Observation 4: This is a rook placement / matching counting problem

A peaceful rook placement means:

- at most one rook in each chosen column;
- at most one rook in each chosen row.

We can process columns one by one.

For each column, either:

1. place no rook in it;
2. place one rook in an allowed unused row.

---

## 3. Full solution approach based on the observations

### Step 1: Compute overlaps

Define a helper function:

```text
overlap([0, A), [B, B + C))
```

It returns:

```text
max(0, min(A, B + C) - max(0, B))
```

For columns:

```text
cs = overlap([0, N), [W, W + M))
c1 = N - cs
c2 = M - cs
```

For rows:

```text
rs = overlap([0, N), [H, H + M))
r1 = N - rs
r2 = M - rs
```

---

### Step 2: Dynamic programming state

Use:

```text
dp[a][b][c]
```

where:

- `a` = number of unused first-only rows;
- `b` = number of unused shared rows;
- `c` = number of unused second-only rows.

Initially all rows are unused:

```text
dp[r1][rs][r2] = 1
```

---

### Step 3: Process columns

For every column, depending on its type, we know which row types it can use.

From state `dp[a][b][c]`, we can:

#### Skip this column

```text
next[a][b][c] += dp[a][b][c]
```

#### Place rook in a first-only row

Possible only if this column allows first-only rows and `a > 0`.

There are `a` choices for the actual row:

```text
next[a - 1][b][c] += dp[a][b][c] * a
```

#### Place rook in a shared row

Possible only if allowed and `b > 0`:

```text
next[a][b - 1][c] += dp[a][b][c] * b
```

#### Place rook in a second-only row

Possible only if allowed and `c > 0`:

```text
next[a][b][c - 1] += dp[a][b][c] * c
```

Process:

```text
c1 first-only columns
cs shared columns
c2 second-only columns
```

---

### Step 4: Extract the answer

At the end, the number of placed rooks is equal to the number of used rows:

```text
(r1 - a) + (rs - b) + (r2 - c)
```

Sum all states where this value equals `K`.

If `K` is too large, no state will match and the answer is `0`.

---

### Complexity

The number of states is at most:

```text
21 × 21 × 21
```

The number of processed columns is at most:

```text
40
```

So the algorithm is easily fast enough.

The answer can be huge, so we need arbitrary-precision integers.

---

## 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 n, m;
int64_t w, h, k;

void read() { cin >> n >> m >> w >> h >> k; }

void solve() {
    // The two boards are the squares [0, n) x [0, n) and [w, w+m) x [h, h+m) in
    // (column, row) coordinates. A rook on cell (c, r) belongs to the union iff
    // it lies in one of the two squares, and peaceful means no two rooks share
    // a column or a row of the infinite grid. So we are counting partial
    // matchings of size k in the bipartite graph "column -- row" restricted to
    // cells that exist in the union.
    //
    // Only columns that carry a cell ever hold a rook, so the relevant columns
    // are exactly those in board1's x-range [0,n) or board2's [w,w+m): at most
    // n+m <= 40 of them, and the enormous empty gap when w is large simply
    // contributes nothing. By x-position a column falls into one of three types
    // - board1 only, board2 only, or both - with counts c_a, c_ab, c_b. Rows
    // likewise split by y-position into board1-only (r1), shared (r_both), and
    // board2-only (r2). Since w and h enter only through the overlap lengths,
    // the huge ranges collapse to these six counts, each at most 20, and a k
    // larger than the rows available just leaves the target rook count
    // unreachable, i.e. zero.
    //
    // A column can host a rook on a row exactly when the two share a board:
    //
    // - a board1-only column reaches board1-only and shared rows;
    // - a board2-only column reaches board2-only and shared rows;
    // - a both-column reaches all three kinds of row.
    //
    // We sweep the relevant columns one at a time, keeping the number of rows
    // still free in each kind as the state (b, m, t). For each column we either
    // skip it or drop a rook on one of its reachable kinds; placing it on a
    // kind with f free rows multiplies by f, the choice of which concrete row
    // to use, and spends one row of that kind. Distinct columns are handled
    // separately, so no column or row is ever reused and each placement is
    // counted once. Every consumed row is one rook, so a state's rook count is
    // fixed by the rows it has spent, and the answer sums the states whose
    // spent rows equal k, kept exact with big integers.

    auto clamp_overlap = [](int64_t lo_b, int64_t len_b, int len_a) -> int {
        int64_t lo = max((int64_t)0, lo_b);
        int64_t hi = min((int64_t)len_a, lo_b + len_b);
        return (int)max((int64_t)0, hi - lo);
    };

    int c_ab = clamp_overlap(w, m, n);
    int c_a = n - c_ab;
    int c_b = m - c_ab;

    int r_both = clamp_overlap(h, m, n);
    int r1 = n - r_both;
    int r2 = m - r_both;

    vector<vector<vector<bigint>>> dp(
        r1 + 1, vector<vector<bigint>>(r_both + 1, vector<bigint>(r2 + 1))
    );
    dp[r1][r_both][r2] = bigint(1);

    auto place_column = [&](bool use_r1, bool use_both, bool use_r2) {
        vector<vector<vector<bigint>>> nd(
            r1 + 1, vector<vector<bigint>>(r_both + 1, vector<bigint>(r2 + 1))
        );
        for(int b = 0; b <= r1; b++) {
            for(int m = 0; m <= r_both; m++) {
                for(int t = 0; t <= r2; t++) {
                    const bigint& v = dp[b][m][t];
                    if(v.isZero()) {
                        continue;
                    }
                    nd[b][m][t] += v;
                    if(use_r1 && b > 0) {
                        nd[b - 1][m][t] += v * b;
                    }
                    if(use_both && m > 0) {
                        nd[b][m - 1][t] += v * m;
                    }
                    if(use_r2 && t > 0) {
                        nd[b][m][t - 1] += v * t;
                    }
                }
            }
        }

        dp = std::move(nd);
    };

    for(int i = 0; i < c_a; i++) {
        place_column(true, true, false);
    }
    for(int i = 0; i < c_ab; i++) {
        place_column(true, true, true);
    }
    for(int i = 0; i < c_b; i++) {
        place_column(false, true, true);
    }

    bigint answer;
    for(int b = 0; b <= r1; b++) {
        for(int m = 0; m <= r_both; m++) {
            for(int t = 0; t <= r2; t++) {
                int64_t rooks = (int64_t)(r1 - b) + (r_both - m) + (r2 - t);
                if(rooks == k) {
                    answer += dp[b][m][t];
                }
            }
        }
    }

    cout << answer << '\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();
        solve();
    }

    return 0;
}
```

---

## 5. Python implementation with detailed comments

```python
# Read input.
# N, M = sizes of the two boards.
# W, H = horizontal and vertical shift of the second board.
# K = number of rooks to place.
N, M, W, H, K = map(int, input().split())


def overlap(start_b, len_b, len_a):
    """
    Return the length of intersection of intervals:

        [0, len_a)
        [start_b, start_b + len_b)

    The first board starts at 0.
    The second board starts at start_b.
    """

    left = max(0, start_b)
    right = min(len_a, start_b + len_b)

    return max(0, right - left)


# Count column types.
shared_cols = overlap(W, M, N)
first_only_cols = N - shared_cols
second_only_cols = M - shared_cols

# Count row types.
shared_rows = overlap(H, M, N)
first_only_rows = N - shared_rows
second_only_rows = M - shared_rows

"""
dp[a][b][c] = number of ways after processing some columns, where:

a = number of unused first-only rows
b = number of unused shared rows
c = number of unused second-only rows
"""
dp = [
    [
        [0 for _ in range(second_only_rows + 1)]
        for _ in range(shared_rows + 1)
    ]
    for _ in range(first_only_rows + 1)
]

# Initially, every row is unused.
dp[first_only_rows][shared_rows][second_only_rows] = 1


def process_column(can_first, can_shared, can_second):
    """
    Process one column.

    can_first:
        whether this column can use first-only rows

    can_shared:
        whether this column can use shared rows

    can_second:
        whether this column can use second-only rows
    """

    global dp

    next_dp = [
        [
            [0 for _ in range(second_only_rows + 1)]
            for _ in range(shared_rows + 1)
        ]
        for _ in range(first_only_rows + 1)
    ]

    for a in range(first_only_rows + 1):
        for b in range(shared_rows + 1):
            for c in range(second_only_rows + 1):
                cur = dp[a][b][c]

                if cur == 0:
                    continue

                # Option 1: place no rook in this column.
                next_dp[a][b][c] += cur

                # Option 2: place a rook in a first-only row.
                if can_first and a > 0:
                    # There are a choices for the actual row.
                    next_dp[a - 1][b][c] += cur * a

                # Option 3: place a rook in a shared row.
                if can_shared and b > 0:
                    # There are b choices for the actual row.
                    next_dp[a][b - 1][c] += cur * b

                # Option 4: place a rook in a second-only row.
                if can_second and c > 0:
                    # There are c choices for the actual row.
                    next_dp[a][b][c - 1] += cur * c

    dp = next_dp


# First-only columns can use first-only and shared rows.
for _ in range(first_only_cols):
    process_column(True, True, False)

# Shared columns can use all row types.
for _ in range(shared_cols):
    process_column(True, True, True)

# Second-only columns can use shared and second-only rows.
for _ in range(second_only_cols):
    process_column(False, True, True)


answer = 0

# Sum all final states where exactly K rooks were placed.
for a in range(first_only_rows + 1):
    for b in range(shared_rows + 1):
        for c in range(second_only_rows + 1):
            placed = (
                (first_only_rows - a)
                + (shared_rows - b)
                + (second_only_rows - c)
            )

            if placed == K:
                answer += dp[a][b][c]

print(answer)
```