Cookie Monster is playing with his three cookies, which have distinct sizes $$$a$$$, $$$b$$$, and $$$c$$$, respectively. He wants to create a cookie tower by stacking them up. Cookie Monster wants the largest cookie to be on the bottom (for a strong foundation). However, he wants the top of the tower to be big as well, so he wants to place the second largest cookie on top and the smallest cookie in the middle.

Given the sizes of Cookie Monster's cookies $$$a$$$, $$$b$$$, and $$$c$$$, help him order the cookies in the order biggest, smallest, middle.

The first line will contain the three integers $$$a$$$, $$$b$$$, and $$$c$$$ ($$$1 \leq a, b, c \leq 100$$$). It is guaranteed that $$$a \neq b$$$, $$$b \neq c$$$, and $$$a \neq c$$$.

Print three integers corresponding to the cookies in the order Cookie Monster wants them in.

We can sort the array $$$[a, b, c]$$$ in decreasing order and then make a swap between the last two numbers and print our result.

blocks = list(map(int, input().split()))
blocks.sort(reverse=True)
temp = 0
temp = blocks[2]
blocks[2] = blocks[1]
blocks[1] = temp
print(" ".join(str(el) for el in blocks))

Cookie Monster has an empty board with $$$m$$$ rows and $$$n$$$ columns. He wants to place as many $$$2 \times 1$$$ dominoes as possible on the board without overlapping. Dominoes can be placed either vertically or horizontally, and each domino must cover exactly two adjacent cells.

Determine whether it is possible for Cookie Monster to tile the entire board with such dominoes so that no cell remains uncovered.

A single line containing two integers $$$m$$$ and $$$n$$$ ($$$1 \leq m, n \leq 100$$$).

Print YES if the board can be completely tiled with $$$2 \times 1$$$ dominoes. Otherwise, print NO.

If the total number of squares is odd, then it is impossible to tile the board as each domino contains an even number of squares. On the other hand, if either $$$m$$$ or $$$n$$$ is even, we can always fill the board with dominoes. Thus, we only have to check if $$$m$$$ or $$$n$$$ is even.

#include <bits/stdc++.h>
using namespace std;
 
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
 
    long long m, n;
    cin >> m >> n;
 
    if ((m * n) % 2 == 0) {
        cout << "YES\n";
    } else {
        cout << "NO\n";
    }
 
    return 0;
}

Look at the sample image given and try some cases on paper to see how the rotation operation works. We can think of the $$$180$$$ degree rotation as reversing the array of digits and rotating each individual digit by $$$180$$$ degrees. Notice that when rotated $$$180$$$, all digits remain the same except for $$$6$$$ and $$$9$$$ which swap. So to solve the problem, we need to reverse the order of the $$$n$$$ digits and swap any $$$6 \leftrightarrow 9$$$.

n = int(input())
nums = list(map(int, input().split()))
nums = nums[::-1]
nums = [6 if el == 9 else (9 if el == 6 else el) for el in nums]
print(" ".join(str(el) for el in nums))

We can find formulas for the index of the desk in the $$$i$$$-th row and $$$j$$$-th column ($$$0$$$-indexed) in both Cookie Monster's and Elmo's ordering.

• Cookie Monster: $$$i \cdot m + j + 1$$$
• Elmo: $$$i + j \cdot n + 1$$$

Setting them equal to each other, we get $$$i \cdot (m - 1) = j \cdot (n - 1)$$$ where $$$0 \le i \lt n$$$ and $$$0 \le j \lt m$$$. Let $$$g = \text{gcd}(m-1, n-1)$$$, $$$a=\frac{m-1}{g}$$$, and $$$b=\frac{n-1}{g}$$$. The equation simplifies to $$$i \cdot a = j \cdot b$$$, where $$$a$$$ and $$$b$$$ are coprime. This means that $$$i$$$ and $$$j$$$ must be multiples of $$$b$$$ and $$$a$$$, respectively. If we let $$$i=k \cdot b$$$, then we can see that $$$j=k \cdot a$$$. Within the constraints on $$$i$$$ and $$$j$$$, there are $$$g+1$$$ values of $$$k$$$ that satisfy the equation. So the answer is $$$\text{gcd}(m-1, n-1) + 1$$$.

#include <bits/stdc++.h>

using namespace std;

int main() {
    int n, m; cin >> n >> m;
    if (n == 1 && m == 1) cout << 1 << endl;
    else if (n == 1) cout << m << endl;
    else if (m == 1) cout << n << endl;
    else cout << 1 + __gcd(n-1, m-1) << endl;
}

Since the XOR operation is associative and commutative, if $$$A$$$ and $$$B$$$ have the same XOR, then the XOR of the whole array must be $$$0$$$. So if the XOR of the whole array is not $$$0$$$, then the answer is $$$0$$$. If the XOR is $$$0$$$, then any valid way to partition the array must result in two sets with equal XOR. The number of ways to partition the array as defined in the statement is $$$2^n - 2$$$.

#include <bits/stdc++.h>
 
using namespace std;
 
int main() {
    int n; cin >> n;
    int x = 0;
    for (int i = 0; i < n; i++) {
        int v; cin >> v;
        x ^= v;
    }
    if (x) cout << 0 << endl;
    else cout << (1LL << n) - 2 << endl;
}

Define a fixed point as an integer $$$i \in [1, n]$$$ such that $$$a_i = i$$$. Let $$$f(a)$$$ be the number of integers $$$j \in [1, n]$$$ which are not fixed points of the permutation $$$a_1, \ldots, a_n$$$.

Claim: Let $$$A^\prime$$$ refer to the permutation $$$A$$$ after applying the given operation once. We must have $$$f(A^\prime) \le \frac{f(A)}{2}$$$.

Proof: If $$$A$$$ has no fixed points, we are done. Otherwise, let $$$i$$$ be the index of the first $$$a_i$$$ which is not a fixed point. We will perform $$$\text{swap}(i, a_i)$$$, which converts $$$a_i$$$ into a fixed point. Additionally, $$$a_{a_i}$$$ is guaranteed to not be a fixed point and is moved to index $$$i$$$, where it may or may not be affected by future swaps. We can then continue to the next $$$a_j$$$ ($$$j \gt i$$$) which is not a fixed point. Each swap affects two numbers and converts at least one of them into a fixed point. Thus, the number of fixed points is at least halved during the operation.

A permutation is sorted when it has $$$n$$$ fixed points or when $$$f(A) = 0$$$. Initially, $$$f(A) \le n$$$, so after $$$O(\log{n})$$$ iterations we are guaranteed to have $$$f(A) = 0$$$. Thus, we can simulate operations until the permutation is sorted for a time complexity of $$$O(n \log n)$$$.

#include <bits/stdc++.h>
using namespace std;

void solve() {
	int n; cin >> n;

    int bound = 32 - __builtin_clz(n);

    vector<int> arr(n);
    for (int i = 0; i < n; i++) {
        cin >> arr[i];
        arr[i]--;
    }

    int cnt = 0;

    do {
        bool ok = true;
        for (int i = 0; i < n-1; i++) {
            if (arr[i] > arr[i+1]) {
                ok = false;
                break;
            }
        }
        if (ok) {
            cout << cnt << endl;
            break;
        }
        
        for (int i = 0; i < n; i++) {
            swap(arr[i], arr[arr[i]]);
        }

        cnt++;
        assert(cnt <= bound);
    } while (true);
}

int main() {
    int t = 1; 
    // cin >> t;
    while (t--) {
        solve();
    }
}

The problem is basically asking us to attach the trees that node $$$1$$$ is not a part of onto the tree that node $$$1$$$ is a part of. Starting with the tree that node $$$1$$$ is part of, it is optimal to keep attaching trees onto the deepest point of this tree. To compute the height of the final tree, we take the height of the initial tree plus the diameters of all the other trees, which is the farthest distance between any two nodes of the tree, as this tells us how deep we can make the tree go. To compute the diameter of each of these trees, we can perform a depth first search from an arbitrary node in the tree to find the farthest node from that node, then find the height of the tree when rooted at this new farthest node.

#include <bits/stdc++.h>
 
using namespace std;
 
using pi = pair<int, int>;
#define f first
#define s second
 
int main() {
	int n, m; cin >> n >> m;
 
    vector<int> vis(n, 0);
    vector<vector<int>> adj(n);
 
    for (int i = 0; i < m; i++) {
        int u, v; cin >> u >> v; 
        u--; v--;
        adj[u].push_back(v);
        adj[v].push_back(u);
    }
 
    auto dfs = [&](auto &&self, int x, int v) -> pi {
        vis[x]++;
        pi res{0, x};
        for (auto y : adj[x]) {
            if (vis[y] != v) continue;
            pi ret = self(self, y, v);
            if (ret.f >= res.f) {
                res = ret;
                res.f++;
            }
        }
        return res;
    };
 
    int ans = n-m-1;
    pi v = dfs(dfs, 0, 0);
    ans += v.f;
 
 
    for (int i = 0; i < n; i++) {
        if (vis[i] == 0) {
            pi z1 = dfs(dfs, i, 0);
            pi z2 = dfs(dfs, z1.s, 1);
            ans += z2.f;
        }
    }
 
    cout << ans << endl;
}

Cookie Monster and Elmo are playing a game on an $$$8 \times 8$$$ chessboard, but Cookie Monster is a prolific fraudster with a fake CS degree from MIT. The game is too hard for him, so Cookie Monster probably plans to cheat using a hidden vibrating device, but before this, he needs your help.

On each square of the board is a number $$$a_{ij}$$$. With one move, either player can increase the value of a single square $$$a_{ij}$$$ by $$$1$$$. First, Cookie Monster performs no more than $$$x$$$ moves on the board. Second, Elmo performs no more than $$$y$$$ moves on the board.

At the end of each of their turns, Cookie Monster's score is the number of $$$a_{ij}$$$ that are prime. Cookie Monster wants to maximize his score, while Elmo wants to minimize it. Determine the maximum score Cookie Monster can attain.

The first line will contain two integers $$$x$$$ and $$$y$$$ ($$$0 \le x, y \le 64$$$).

The following $$$8$$$ lines will each contain $$$8$$$ integers $$$a_{ij}$$$ ($$$3 \lt a_{ij} \le 10^6$$$).

Output a single integer, the maximum possible attainable score for Cookie Monster if both players act optimally.

Since $$$a_{ij} \gt 1$$$, increasing a prime number by 1 will make it a non-prime, as it'll become divisible by $$$2$$$. Thus, Elmo's strategy is to increase as many prime numbers on the board as he can with exactly one operation, leaving all non-primes untouched.

Cookie Monster's strategy is based on trying to maximize the number of primes on the board. Let's first precompute all primes $$$\le 10^6$$$ with the Sieve of Eratosthenes. Then, we can find for each $$$a_{ij}$$$ the minimum number of operations to make it prime, which we will call $$$d_{ij}$$$. We can use a greedy algorithm to maximize the number of primes: keep taking the smallest $$$d_{ij}$$$ and converting it into prime until you run out of operations. We can do this by sorting $$$d_{ij}$$$, whcih gives a final time complexity of $$$O(n \log n + M \log \log M)$$$, where $$$n$$$ is the size of the board and $$$M$$$ is the maximum $$$a_{ij}$$$.

#include <bits/stdc++.h>
using namespace std;


int main() {
	ios_base::sync_with_stdio(false); cin.tie(nullptr);
	int x, y, sum = 0, tot, a;
    cin >> x >> y;
    vector<int> v;
    vector<bool> c(1'000'200);
    for (int i = 2; i * i < 1'000'200; i++) {
        if (!c[i]) {
            for (int j = i + i; j < 1'000'200; j += i) {
                c[j] = true;
            }
        }
    }
    for (int i = 0; i < 8; i++) {
        for (int j = 0; j < 8; j++) {
            cin >> a;
            tot = 0;
            for (int k = 0; k < 70; k++) {
                if (!c[a]) {
                    break;
                }
                a++;
                tot++;
            }
            v.push_back(tot);
        }
    }
    sort(v.begin(), v.end());
    for (int i = 0; i < 64; i++) {
        if (x >= v[i]) {
            x -= v[i];
            sum++;
        }
    }
    cout << max(sum - y, 0);
	return 0;
}

In this variation, a few critical details have been changed. We see that Elmo's can still make a number no longer "valid" in exactly one move using Fermat's Theorem on Sums of Two Squares.

Using the same logic from the easy variant, and combining our intuition from the aforementioned theorem, our main hurdle is determining quickly whether a very large number is prime. We use the deterministic variant of the Miller Rabin algorithm to overcome this. This algorithm is supported by a binary exponentiation algorithm, as well as utilizes 128-bit integers to avoid overflow with the scale of our grid's values. Besides this, we conduct a greedy and sorting approach as detailed previously.

#include <bits/stdc++.h>
using namespace std;

using u64 = uint64_t;
using u128 = __uint128_t;

u64 binpower(u64 base, u64 e, u64 mod) {
    u64 result = 1;
    base %= mod;
    while (e) {
        if (e & 1)
            result = (u128)result * base % mod;
        base = (u128)base * base % mod;
        e >>= 1;
    }
    return result;
}

bool check_composite(u64 n, u64 a, u64 d, int s) {
    u64 x = binpower(a, d, n);
    if (x == 1 || x == n &mdash; 1)
        return false;
    for (int r = 1; r < s; r++) {
        x = (u128)x * x % n;
        if (x == n &mdash; 1)
            return false;
    }
    return true;
};

bool MillerRabin(u64 n) {
    if (n < 2)
        return false;

    for (int a : {2, 3, 5, 13, 19, 73, 193, 407521, 299210837}) {
        if (n == a)
            return true;
    }

    int r = 0;
    u64 d = n - 1;
    while ((d & 1) == 0) {
        d >>= 1;
        r++;
    }

    for (int a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (n == a)
            return true;
        if (check_composite(n, a, d, r))
            return false;
    }
    return true;
}

int main() {
	ios_base::sync_with_stdio(false); cin.tie(nullptr);
	int x, y, sum = 0, tot;
    long long a;
    cin >> x >> y;
    vector<int> v;
    for (int i = 0; i < 8; i++) {
        for (int j = 0; j < 8; j++) {
            cin >> a;
            if (a % 2 == 1) {
                tot = 0;
                if (a % 4 == 3) {
                    tot += 1;
                    a += 2;
                }
                while (!MillerRabin(a)) {
                    tot += 2;
                    a += 4;
                }
                v.push_back(tot);
            } else if (a == 2) {
                v.push_back(0);
            } else {
                v.push_back(1000);
            }
        }
    }
    sort(v.begin(), v.end());
    for (int i = 0; i < 64; i++) {
        if (x >= v[i]) {
            x -= v[i];
            sum++;
        }
    }
    cout << max(sum &mdash; y, 0);
	return 0;
}

Consider a slow method of answering a single query. Binary search on the first time the range is less than $$$x$$$. When trying to check a specific time, let's say $$$T$$$, first build a segment/Fenwick tree with the initial values. Then apply all updates until $$$T$$$, then finally query for the range in question and compare that with $$$x$$$.

The issue with the previous solution is that we are repeatedly doing the process of building the segment tree and then updating from time $$$1 \ldots t$$$. We can expand on this solution by using a sweepline, performing the binary searches in parallel. Repeatedly perform the process of building the segment tree and updating times $$$1 \ldots t$$$. Whenever you reach a time $$$T$$$ that needs to be checked for a query's binary search, do a range sum on the tree and update that query's binary search state, and then mark the new $$$T$$$ that needs to be checked. If the binary search has terminated, set that value as the answer for the query.

There is a more detailed explanation here.

There are other alternative solutions, such as square root decomposition and powerful data structures (offline 2D BIT).

#include <bits/stdc++.h>
using namespace std;
 
using ll = long long;
using pl = pair<ll, ll>;
#define f first
#define s second
 
struct SegTree {
    int len;
    vector<ll> arr;
    SegTree(int _len) : len(_len), arr(_len * 2, 0) {}
    void set(int ind, ll val) {
        //cerr<<"set: " << ind << " " << val << endl;
        ind += len; arr[ind] = val;
        for (; ind > 1; ind /= 2) arr[ind / 2] = arr[ind] + arr[ind ^ 1];
    }
    ll query(int start, int end) {
        // cerr<<"query: "<< start << " " << end;
        ll res = 0;
        for (start += len, end += len; start < end; start /= 2, end /= 2) {
            if (start % 2 == 1) { res += arr[start++]; }
            if (end % 2 == 1) { res += arr[--end]; }
        }
        //cerr << " " << res << endl;
        return res;
    }
};
 
int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
 
    int n, q, t; cin >> n >> t >> q;
    vector<ll> arr(n);
    for (auto &x : arr) cin >> x;
 
    vector<pl> ops(t);
    for (auto &x : ops) cin >> x.f >> x.s;
    for (auto &x : ops) x.f--;
 
    vector<array<ll, 3>> queries(q);
    for (auto &x : queries) cin >> x[0] >> x[1] >> x[2];
    for (auto &x : queries) {x[0]--; x[1]--;}
 
    vector<pl> bound(q, {-1, t+1});
 
    int rounds = 64 - __builtin_clz(q) + 1;
    while (rounds--) {
        vector<int> check[t+1];
        for (int i = 0; i < q; i++) {
            int m = (bound[i].f + bound[i].s) / 2;
            if (m >= 0 && m <= t) check[m].push_back(i);
        }
 
        vector<ll> car(arr);
        SegTree st(n);
        for (int i = 0; i < n; i++) st.set(i, car[i]);
 
        for (auto j : check[0]) {
            ll res = st.query(queries[j][0], queries[j][1]+1);
            if (res <= queries[j][2]) bound[j].s = 0;
            else bound[j].f = 0;
        }
 
        for (int i = 0; i <= t; i++) {
            for (auto j : check[i]) {
                ll res = st.query(queries[j][0], queries[j][1]+1);
                if (res <= queries[j][2]) bound[j].s = i;
                else bound[j].f = i;
            }
            if (i == t) break;
            car[ops[i].f] -= ops[i].s;
            st.set(ops[i].f, car[ops[i].f]);
        }
    }
 
    for (auto &lr : bound) cout << (lr.s == t+1 ? -1 : lr.s) << endl;
}

We will use intuition from linear algebra to approach this problem. We will use the notation that indexes are always taken $$$\bmod n$$$ and $$$0$$$-indexed. In other words, $$$x_i$$$ and other variables are defined as $$$x_{(i \bmod n)}$$$ for integers $$$i$$$ outside the range $$$[0, n-1]$$$.

Let $$$c_i$$$ be the total amount that you increase the $$$i$$$-th lock by. For each integer $$$i \in [0, n-1]$$$, we have an equation $$$\sum\limits_{j=i}^{i+m} c_{j} = a_{i+m}$$$, forming a system of $$$n$$$ linear equations. If we subtract all pairs of equations with adjacent $$$i$$$, then we get $$$n$$$ equations of the form $$$c_{i+m+1} - c_i = a_{i+m+1} - a_{i+m}$$$. Our new system of equations contains $$$2n$$$ variables instead of $$$mn$$$, making it significantly simpler to solve. However, our new system contains less information as we cannot represent the original system as a linear combination of the equations in our new system. To fix this, we can add one equation from our original system (i.e. $$$c_0 + \ldots + c_m = a_m$$$) to our new system at the end, which we will refer to as the "final" equation.

Let's define $$$d_{i} = a_{i+m+1} - a_{i+m}$$$ for convenience. We now have our system of equations $$$c_{i+m+1} - c_i = d_{i}$$$ for each integer $$$i \in [0, n-1]$$$. Let's consider this as a functional graph with an edge from every node $$$i$$$ to node $$$(i+m+1) \bmod n$$$ with weight $$$d_i$$$. This graph has several useful properties. First, it consists of $$$g=\text{gcd}(m+1, n)$$$ disjoint cycles of equal length. Additionally, if there is a path from node $$$i$$$ and to node $$$j$$$ with sum of the edge weights equal to $$$x$$$, then $$$c_j - c_i \equiv x \pmod {998\,244\,353}$$$. Thus, if we traverse one full cycle, the edge weights along our path must add to $$$0 \pmod {998\,244\,353}$$$. If any cycle does not satisfy this condition, then no solution exists and we can terminate our algorithm. Otherwise, we can always find values $$$c_i$$$ that satisfy the constraints. If we fix the value $$$c_i$$$ for one node in a cycle, it determines the values of $$$c_i$$$ for the remaining nodes in the cycle. In linear algebra terms, each cycle, if satisfiable, consists of $$$1$$$ free variable and $$$\frac{n}{g} - 1$$$ basic variables.

After processing each cycle, we have a total of $$$g$$$ free variables and all other variables $$$c_i$$$ are determined from one of these free variables plus a constant offset. Now, we need to consider the "final" equation. This equation essentially reduces the number of our free variables by $$$1$$$ by allowing us to write one of our free variables as a linear combination of the others plus some constant. Since $$$g \ge 1$$$, we will always have at least one free variable to accommodate this equation. After this step, we have fully solved our system of equations.

To output a construction as the problem asks, we can set each the free variables to any value and solve for the rest of the variables. For instance, the model solution sets all the free variables to $$$0$$$. The final time complexity is $$$O(n)$$$.

#include <bits/stdc++.h>

using namespace std;

using ll = long long;
const int M = 998'244'353;

ll modpow(ll x, ll p) {
    ll a = 1;
    while (p) {
        if (p & 1) a = (a * x) % M;
        x = (x * x) % M;
        p >>= 1;
    }
    return a;
}

ll modinv(ll x) {
    return modpow(x, M-2);
}



int main() {
	int n, m; cin >> n >> m;
    m++;
    int g = __gcd(m, n);

    // input, adjacent differences
    ll a[n], b[n];
    for (int i = 0; i < n; i++) cin >> a[i];
    for (int i = 0; i < n; i++) b[i] = (a[(i+1)%n] - a[i] + M) % M;

    // free variable, difference
    // c[i] = c[fv[i]] + d[i]
    int fv[n];
    ll d[n];
    for (int i = 0; i < g; i++) {
        fv[i] = i; d[i] = 0LL;

        int x = i, y = (i + m) % n;
        while (y != i) {
            fv[y] = i;
            d[y] = (d[x] + b[x]) % M;
            x = y;
            y = (y + m) % n;
        }

        if ((d[x] + b[x]) % M != d[i]) {
            cout << "NO" << endl;
            return 0;
        }
    }

    cout << "YES" << endl;

    // construction
    ll x = a[0];
    for (int i = 0; i < m; i++) x += M - d[i];
    x %= M;

    int fac = m / g;
    ll sm = modinv((ll)m/g) * x % M; // c1 + ... + cg = sm (mod M)

    ll c[n];      
    c[0] = sm;
    for (int i = 1; i < g; i++) c[i] = 0;
    for (int i = g; i < n; i++) c[i] = (c[fv[i]] + d[i]) % M;

    for (int i = 0; i < n; i++) cout << c[(i + m - 1) % n] << " \n"[i==n-1];

    // verify
    ll arr[2*n];
    for (int i = 0; i < 2*n; i++) arr[i] = 0;
    for (int i = 0; i < n; i++) {
        arr[i] += c[(i + m - 1) % n];
        arr[i + m] -= c[(i + m - 1) % n];
    }
    for (int i = 1; i < 2*n; i++) arr[i] += arr[i-1];
    for (int i = 0; i < n; i++) arr[i] += arr[i+n];

    for (int i = 0; i < n; i++) {
        arr[i] %= M;
        assert(arr[i] == a[i]);
        // cerr << arr[i] << " \n"[i==n-1];
    }
}

We can notice that the available seats for Person $$$i$$$ depends only on where everyone with a higher number sat. This motivates us to have people sit from highest to lowest order.

We can approach this problem with dynamic programming. Let $$$\text{dp}[i][j]$$$ be the number of ways to add the $$$i$$$ highest people so that there are $$$j$$$ valid places for the next highest person and at least $$$1$$$ invalid place. The lattermost condition ensures that there are exactly $$$2$$$ ways to increase $$$j$$$ or keep it the same: one on either side of the current valid range.

In other words, there is a contribution of $$$2 \text{dp}[i][j]$$$ for each $$$\text{dp}[i+1][j+\ell]$$$ for $$$0 \le \ell \lt k$$$. The number of ways to reduce $$$j$$$ by $$$1$$$ is $$$j - 2k$$$, so there is a contribution of $$$(j - 2k) \cdot \text{dp}[i][j]$$$ to $$$\text{dp}[i+1][j-1]$$$. We might try to transition to a state $$$\text{dp}[i+1][j+\ell]$$$ where $$$i + 1 + j + \ell \ge n$$$, which means that now all remaining places are valid for future people. This contributes $$$(n-i-1)! \cdot \text{dp}[i+1][j+\ell]$$$ to our final answer, and we can ignore it in our subsequent transitions.

#include <bits/stdc++.h>
using namespace std;
 
const int mod = 998'244'353;
 
long long ex(long long a, long long b) {
    if (b == 0) {
        return 1;
    }
    long long r = ex(a, b / 2);
    r *= r;
    r %= mod;
    if (b % 2 == 1) {
        r *= a;
        r %= mod;
    }
    return r;
}
 
int main() {
	ios_base::sync_with_stdio(false); cin.tie(nullptr);
    int n, k;
    cin >> n >> k;
    // long long dp[n+k][n+1], p[n+k][n+1];
    // memset(dp, 0, sizeof(dp));
    // memset(p, 0, sizeof(p));
    vector<vector<long long>> dp(n + k, vector<long long>(n + 1)), p(n + k, vector<long long>(n + 1));
    if (2 * k + 1 > n) {
        cout << 1;
    } else {
        for (int j = 2 * k; j < n; j++) {
            for (int i = n + k - 1; i >= 2 * k; i--) {
                if (i != n + k - 1) {
                    p[i][j] += p[i + 1][j];
                    p[i][j] %= mod;
                }
                if (i >= j) {
                    dp[i][j] = 1;
                } else {
                    long long den = ex(j, mod - 2);
                    dp[i][j] = (((((i - 2 * k) * dp[i - 1][j - 1]) % mod) + p[i][j]) * den) % mod;
                    dp[i][j] %= mod;
                }
                p[i][j + 1] += 2 * dp[i][j];
                p[i][j + 1] %= mod;
                p[i - k][j + 1] -= 2 * dp[i][j];
                p[i - k][j + 1] += 2 * mod;
                p[i - k][j + 1] %= mod;
            }
        }
        cout << dp[2 * k][n - 1];
    }
    return 0;
}

We can notice that the availible seats for Person $$$i$$$ depends only on where everyone with a higher number sat. This motivates us to have people sit from highest to lowest order.

Define $$$p_i$$$ as the probability that $$$i$$$ people sitting at a table with $$$i$$$ seats satisfy the condition. Denote $$$a_i$$$ as the number of valid ways people numbered $$$1$$$ through $$$i-1$$$ can be arranged in a table with $$$i$$$ seats after fixing person $$$i$$$'s seat. We can see that $$$p_i = \frac{a_i}{(i-1)!}$$$. Our goal is to compute $$$p_n$$$.

Let's see how we can compute $$$a_i$$$. Our trivial base cases are $$$p_i = 1$$$ and $$$a_i = (i-1)!$$$ for $$$1 \le i \le 2k+1$$$, because then any configuration satisfies the condition.

Claim 1: Suppose Person $$$x$$$ and Person $$$y$$$ are seated and have $$$m$$$ empty seats between them. There are $$$m$$$ people, each with a number smaller than $$$\text{min}(x,y)$$$, who want to fill in these seats. The number of ways that they can do this is $$$a_{m+1}$$$.

Proof: Every seat other than the $$$m$$$ empty ones and the $$$2$$$ occupied ones are irrelevant. If we remove them from the problem and merge the two seated people, we now have a $$$m+1$$$ seats with the highest-numbered person already seated. By definition, the number of ways to arrange the remaining $$$m$$$ people is $$$a_{m+1}$$$.

Claim 2: Suppose Person $$$i$$$ and Person $$$i-1$$$ have sat a distance $$$d \le k$$$ from each other on a table with $$$i \ge 3$$$ seats. The number of ways to place the remaining $$$i-2$$$ people is given by $$$a_{i-d} \cdot {i-2 \choose d-1} \cdot (d-1)! = a_{i-d} \cdot \frac{(i-2)!}{(i-d-1)!}$$$

Proof: We choose any $$$d-1$$$ people to put between Person $$$n-1$$$ and Person $$$n$$$, and there are $$$(d-1)!$$$ ways to order them. Since $$$d \le k$$$, these $$$d-1$$$ people will always be covered) For the remaining $$$n-d-1$$$ seats, we can apply Claim 1 and determine there are $$$a_{n-d}$$$ ways to do so.

Let's use Claim 2 to find a recurrence formula for $$$a_i$$$ for $$$i \gt 2k+1$$$. If we have a table with $$$i$$$ seats and we fix person $$$i$$$, for each $$$d$$$ from $$$1$$$ to $$$k$$$ inclusive there are $$$2$$$ ways to place person $$$i-1$$$ a distance $$$d$$$ away. Thus, our recurrence is

This formula simplifies further if we use $$$p_i$$$ instead of $$$a_i$$$. Since $$$p_i = \frac{a_i}{(i-1)!}$$$, we have

.

Thus, we have proved that $$$p_i$$$ is equal to $$$\frac{2}{i-1}$$$ times the sum of the previous $$$k$$$ values of $$$p_i$$$. The previous $$$k$$$ values can be maintained with a sliding window technique, allowing us to compute $$$p_n$$$ in $$$O(n)$$$.

#include <bits/stdc++.h>
using namespace std;

const int mod = 998'244'353;

long long ex(long long a, long long b) {
    if (b == 0) {
        return 1;
    }
    long long r = ex(a, b / 2);
    r *= r;
    r %= mod;
    if (b % 2 == 1) {
        r *= a;
        r %= mod;
    }
    return r;
}

int main() {
	ios_base::sync_with_stdio(false); cin.tie(nullptr);
	int n, k;
    long long sum;
    cin >> n >> k;
    vector<long long> dp(n);
    if (2 * k > n) {
        cout << 1;
    } else {
        for (int i = 0; i < 2 * k; i++) {
            dp[i] = 1;
        }
        sum = 2 * k;
        for (int i = 2 * k; i < n; i++) {
            dp[i] = (sum * ex(i, mod - 2)) % mod;
            dp[i] %= mod;
            sum += 2 * dp[i];
            sum -= 2 * dp[i - k];
            sum += 2 * mod;
            sum %= mod;
        }
        cout << dp[n - 1];
    }
	return 0;
}

Given an array of length $$$n$$$ and a number $$$m$$$, there is a straightforward way to find the maximum bitwise AND of a length $$$k$$$ subsequence. Iterate from the most significant bit to the lowest. If there are $$$k$$$ or more numbers in the array with the current bit set to $$$1$$$, remove any numbers that have the bit set to $$$0$$$. Otherwise, do nothing. The answer is the bitwise AND of the remaining numbers.

To solve the original problem, we will define $$$\text{dp}[i][j]$$$ ($$$1 \le i \le m$$$, $$$k \le j \le n$$$) to be the number of ways to fill in the $$$i$$$ highest bits in an array of length $$$n$$$ so that the algorithm above ends with $$$j$$$ numbers remaining whose bitwise AND is equal to the $$$i$$$ highest bits of $$$x$$$. Transitions from $$$i$$$ to $$$i+1$$$ depend on the value of the $$$i+1$$$-th bit in $$$x$$$.

If the $$$i+1$$$-th bit of $$$x$$$ is $$$0$$$, we need to make sure that there are less than $$$k$$$ numbers in the array with $$$i+1$$$-th bit set to $$$1$$$. The remaining $$$n-j$$$ numbers that we discarded earlier can have their bits set to either $$$0$$$ or $$$1$$$.

The factor that each $$$\text{dp}[i][j]$$$ is multiplied by depends only on $$$j$$$ and can be precomputed beforehand.

If the $$$i+1$$$-th bit of $$$x$$$ is $$$1$$$, we can reduce the size of our group. The discarded $$$n-j$$$ numbers again can be either $$$0$$$ or $$$1$$$.

We can optimize this transition with FFT/NTT. Letting $$$d = \ell - j$$$, we can rewrite the transition as:

Since our transitions are either $$$O(n)$$$ or $$$O(n \log n)$$$ and we have $$$m$$$ layers of our $$$\text{dp}$$$, the final time complexity is $$$O(nm \log n)$$$.

#include <bits/stdc++.h>

using namespace std;

template <const long long mod>
struct Mod_Int {
  int x;
  Mod_Int() : x(0) {}
	Mod_Int(int y) : x(y >= 0 ? y % mod : (mod - (-y) % mod) % mod) {}
  Mod_Int(long long y) : x((int)y >= 0 ? (int)y % mod : (mod - (-(int)y) % mod) % mod) {}
  static int get_mod() { return mod; }
  Mod_Int &operator+=(const Mod_Int &p) {
    if ((x += p.x) >= mod) x -= mod;
    return *this;
  }
  Mod_Int &operator-=(const Mod_Int &p) {
    if ((x += mod - p.x) >= mod) x -= mod;
    return *this;
  }
  Mod_Int &operator*=(const Mod_Int &p) {
    x = (int)(1LL * x * p.x % mod);
    return *this;
  }
  Mod_Int &operator/=(const Mod_Int &p) {
    *this *= p.inv();
    return *this;
  }
  Mod_Int operator-() const { return Mod_Int(-x); }
  Mod_Int operator+(const Mod_Int &p) const { return Mod_Int(*this) += p; }
  Mod_Int operator-(const Mod_Int &p) const { return Mod_Int(*this) -= p; }
  Mod_Int operator*(const Mod_Int &p) const { return Mod_Int(*this) *= p; }
  Mod_Int operator/(const Mod_Int &p) const { return Mod_Int(*this) /= p; }
  bool operator==(const Mod_Int &p) const { return x == p.x; }
  bool operator!=(const Mod_Int &p) const { return x != p.x; }
  Mod_Int inv() const {
    assert(*this != Mod_Int(0));
    return pow(mod - 2);
  }
  Mod_Int pow(long long k) const {
    Mod_Int now = *this, ret = 1;
    for (; k > 0; k >>= 1, now *= now) {
      if (k & 1) ret *= now;
    }
    return ret;
  }
  friend ostream &operator<<(ostream &os, const Mod_Int &p) {
    return os << p.x;
  }
  friend istream &operator>>(istream &is, Mod_Int &p) {
    long long a;
    is >> a;
    p = Mod_Int<mod>(a);
    return is;
  }
};
 
template <typename T> struct NumberTheoreticTransform {
  static int max_base;
  static T root;
  static vector<T> r, ir;
  NumberTheoreticTransform() {}
  static void init() {
    if (!r.empty()) return;
    int mod = T::get_mod();
    int tmp = mod - 1;
    root = 2;
    while (root.pow(tmp >> 1) == 1) root += 1;
    max_base = 0;
    while (~tmp & 1) tmp >>= 1, max_base += 1;
    r.resize(max_base), ir.resize(max_base);
    for (int i = 0; i < max_base; i++) {
      r[i] = - root.pow((mod - 1) >> (i + 2));
      ir[i] = r[i].inv();
    }
  }
  static void ntt(vector<T> &a) {
    init();
    int n = (int) a.size();
    assert((n & (n - 1)) == 0);
    assert(n <= (1 << max_base));
    for (int k = n; k >>= 1;) {
      T w = 1;
      for (int s = 0, t = 0; s < n; s += 2 * k) {
        for (int i = s, j = s + k; i < s + k; i++, j++) {
          T x = a[i], y = w * a[j];
          a[i] = x + y, a[j] = x - y;
        }
        w *= r[__builtin_ctz(++ t)];
      }
    }
  }
  static void intt(vector<T> &a) {
    init();
    int n = (int) a.size();
    assert((n & (n - 1)) == 0);
    assert(n <= (1 << max_base));
    for (int k = 1; k < n; k <<= 1) {
      T w = 1;
      for (int s = 0, t = 0; s < n; s += 2 * k) {
        for (int i = s, j = s + k; i < s + k; i++, j++) {
          T x = a[i], y = a[j];
          a[i] = x + y, a[j] = (x - y) * w;
        }
        w *= ir[__builtin_ctz(++ t)];
      }
    }
    T iv = T(n).inv();
    for (auto &x : a) x *= iv;
  }
  static vector<T> convolution(vector<T> a, vector<T> b) {
    int k = (int) a.size() + (int) b.size() - 1, n = 1;
    while (n < k) n <<= 1;
    a.resize(n), b.resize(n);
    ntt(a), ntt(b);
    for (int i = 0; i < n; i++) a[i] *= b[i];
    intt(a), a.resize(k);
    return a;
  }
};
template <typename T> int NumberTheoreticTransform<T>::max_base = 0;
template <typename T> T NumberTheoreticTransform<T>::root = T();
template <typename T> vector<T> NumberTheoreticTransform<T>::r  = vector<T>();
template <typename T> vector<T> NumberTheoreticTransform<T>::ir = vector<T>();

using mint = Mod_Int<998244353>;
using NTT = NumberTheoreticTransform<mint>;

const int mxn = 2e5+5;

mint fc[mxn], iv[mxn], pow2[mxn];
void init() {
	fc[0] = mint(1);	pow2[0] = mint(1);
	for (int i = 1; i < mxn; i++) fc[i] = fc[i-1] * mint(i);
	for (int i = 0; i < mxn; i++) iv[i] = fc[i].inv();
	for (int i = 1; i < mxn; i++) pow2[i] = pow2[i-1] * mint(2);
}
mint choose(int n, int k) {
  assert(n >= 0 && k >= 0 && k <= n);
	return fc[n] * iv[k] * iv[n-k];
}

int main() {
	init();

	int n, m, k, x; 
	cin >> n >> m >> k >> x;

	vector<mint> coefs(n+1);
	mint cur = pow2[k] - mint(1);
	for (int i = k; i <= n; i++) {
		coefs[i] = cur * pow2[n-i];
		cur *= 2; 
		cur -= choose(i, k-1);
	}

	vector<mint> dp(n+1);
	dp[n] = mint(1);

	for (int i = m-1; i >= 0; i--) {
		if (x & (1 << i)) {
			for (int j = 0; j <= n; j++) dp[j] *= pow2[n-j] * fc[j];

			vector<mint> invs(n+1);
			for (int j = 0; j <= n; j++) invs[n-j] = iv[j];

			vector<mint> res = NTT::convolution(dp, invs);
			for (int j = n; j <= 2*n; j++) dp[j-n] = res[j];
			for (int j = 0; j <= n; j++) {
				if (j < k) dp[j] = mint(0); 
				else dp[j] *= iv[j];
			}
		} else {
			for (int j = 0; j <= n; j++)[user:flight][user:treewave] dp[j] *= coefs[j];
		}
	}

	mint ans = 0;
	for (int i = k; i <= n; i++) ans += dp[i];
	cout << ans << endl;
}

To be added by flight

To be added by flight

To be added by flight

To be added by flight