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

221. Big Bishops
time limit per test: 0.25 sec.
memory limit per test: 65536 KB
input: standard
output: standard




A bishop is a piece used in the game of chess which is played on a board of square grids. A bishop can only move diagonally from its current position and two bishops attack each other if one is on the path of the other.

Given two numbers n and k, your job is to determine the number of ways one can put k bishops on an n × n chessboard so that no two of them are in attacking positions.

Input

The input file contains two integers n (1 ≤ n ≤ 50) and k (0 ≤ k ≤ n2).

Output

Print a line containing the total number of ways one can put the given number of bishops on a chessboard of the given size so that no two of them are in attacking positions.

Sample test(s)

Input
Test #1

4 4

Test #2

8 6

Output
Test #1

260

Test #2

5599888
Author:	Andrew Stankevich
Resource:	Little Chess Pieces Series, SPb IFMO 2003-2004 Authumn Training Sessions
Date:	2003-09-27

<|response|>
1. Abridged Problem Statement  
Given an n×n chessboard and an integer k, count the number of ways to place k bishops so that no two attack each other. Bishops attack along diagonals. Output the exact count (it may be large).

2. Key Observations  
- A bishop on a white square can never attack a bishop on a black square. The board splits into two independent “parities.”  
- On each color, if you list all diagonals in one direction (say “\”-diagonals) of that color, no two bishops can occupy the same such diagonal, and choosing one square on each selected diagonal already guarantees they don’t share the other “/”-diagonal.  
- Therefore, for each color we can reduce the problem to: given a list of diagonal lengths v₁, v₂, …, vₘ, how many ways to choose j diagonals and put one bishop on each so that no two attack?  
- That subproblem admits a simple DP: dp[t][j] = ways to place j bishops among the first t diagonals; for diagonal t of length vₜ you either skip it or place one bishop, avoiding squares attacked by the already chosen j–1 bishops on other diagonals.

3. Full Solution Approach  
a. Build two lists per_diag_even and per_diag_odd of diagonal lengths (in “\” direction) for the two square colors:  
   - Diagonals are indexed d = 0 … 2n–2, length = (d<n ? d+1 : 2n–1–d).  
   - Even d go to per_diag_even, odd d go to per_diag_odd.  
b. For each list per_diag of length m, build a DP table dp of size (m+1)×(k+1):  
   - dp[0][0] = 1; dp[0][j>0] = 0.  
   - For t from 1 to m, let v = per_diag[t–1]. Then for j=0..k:  
       dp[t][j] = dp[t–1][j]                           (skip diagonal t)  
                + dp[t-1][j-1] * max(0, v – (j–1))    (place one bishop, avoid j–1 attacked squares)  
c. After computing dp_even and dp_odd, the answer is  
   sum_{i=0..k} dp_even[m_even][i] * dp_odd[m_odd][k–i].  

Time complexity O(n·k). We use big integers because counts grow large.

4. C++ Solution
```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, k;

void read() { cin >> n >> k; }

vector<bigint> solve_one_parity(int p) {
    vector<int> per_diag;
    for(int i = p; i < n; i += 2) {
        per_diag.push_back(i + 1);
        if(i != n - 1) {
            per_diag.push_back(i + 1);
        }
    }

    vector<bigint> prev(k + 1, bigint(0)), cur(k + 1, bigint(0));
    prev[0] = bigint(1);
    for(int v: per_diag) {
        cur.assign(k + 1, bigint(0));
        cur[0] = bigint(1);
        for(int j = 1; j <= min(k, v); j++) {
            cur[j] = prev[j] + prev[j - 1] * (v - j + 1);
        }

        swap(prev, cur);
    }

    return prev;
}

void solve() {
    // Bishops attack along the two diagonal directions. Rotating the board 45
    // degrees turns the diagonals into independent rows so the board splits by
    // cell colour ((r+c) parity) into two sub-boards whose diagonals never
    // interact. For one colour the available cells form anti-diagonals whose
    // lengths are 1, 2, 2, 3, 3, ... (per_diag): processing these diagonals in
    // increasing length order, dp[j] = number of ways to place j non-attacking
    // bishops, with the transition cur[j] = prev[j] + prev[j-1] * (v - j + 1),
    // since a new diagonal of length v offers (v - j + 1) free squares given
    // j-1 already placed on shorter diagonals. Counts grow factorially so we
    // accumulate them with the vendored bigint. The final answer for k bishops
    // is the convolution sum over i of black[i] * white[k - i].

    if(k > 2 * n - 1) {
        cout << 0 << '\n';
        return;
    }

    vector<bigint> data_0 = solve_one_parity(0);
    vector<bigint> data_1 = solve_one_parity(1);

    bigint ans(0);
    for(int i = 0; i <= k; i++) {
        ans += data_0[i] * data_1[k - i];
    }

    cout << ans << '\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
sys.setrecursionlimit(10000)

def build_dp(n, k, start):
    """
    Build dp table for one parity of diagonals.
    Returns dp[t][j] = ways to place j bishops among first t diagonals.
    start = 0 for even diagonals, 1 for odd.
    """
    # Collect lengths of diagonals in "\" direction for this parity
    per_diag = []
    for d in range(start, 2*n - 1, 2):
        if d < n:
            length = d + 1
        else:
            length = 2*n - 1 - d
        per_diag.append(length)

    m = len(per_diag)
    # Initialize dp with zeros
    dp = [ [0]*(k+1) for _ in range(m+1) ]
    dp[0][0] = 1  # one way to place zero bishops

    # Fill DP
    for t in range(1, m+1):
        v = per_diag[t-1]     # length of diagonal t
        dp[t][0] = 1          # zero bishops => always one way
        # For j>=1, either skip this diagonal or place one bishop
        for j in range(1, k+1):
            # Skip diagonal t
            dp[t][j] = dp[t-1][j]
            # Place one bishop: v squares minus (j-1) attacked squares
            if j <= v:
                dp[t][j] += dp[t-1][j-1] * (v - (j - 1))
    return dp

def solve():
    n, k = map(int, sys.stdin.readline().split())
    # Too many bishops for one color => zero ways
    if k > 2*n - 1:
        print(0)
        return

    # Build dp for even and odd diagonals
    dp_even = build_dp(n, k, 0)
    dp_odd  = build_dp(n, k, 1)
    me = len(dp_even) - 1
    mo = len(dp_odd)  - 1

    # Combine counts: sum over distributing k bishops among the two parities
    total = 0
    for i in range(k+1):
        total += dp_even[me][i] * dp_odd[mo][k - i]
    print(total)

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

Explanation of the DP transition:  
- dp[t-1][j] counts ways if we skip diagonal t.  
- dp[t-1][j-1] counts ways if we place one bishop on diagonal t; there are v squares on that diagonal, but the j–1 bishops already chosen occupy one square on each of j–1 other “/”-diagonals, each attacking exactly one square here, so only v–(j–1) squares remain safe. Multiplying by that gives the new placements.