Let \(f_{k, i, j}\) be the number of combinations of \(k\) operations after which the rook is at \((i, j)\). The answer is \(f_{K, x_2, y_2}\).

If \(i = x_1\) and \(j = y_1\), then \(f_{0, i, j} = 1, \), and otherwise \(f_{0, i, j} = 0\). For \(k > 0\), the following relation holds.

• \(f_{k, i, j} = \sum_{x \neq i} f_{k - 1, x, j} + \sum_{y \neq j} f_{k - 1, i, y}\)

Computing this naively requires a total of \(\Theta (N^3 K)\) computation time, so we have to optimize it using the property of \(f\).

Actually calculating \(f_{k, i, j}\), you will observe that most of the values are the same, and for each \(k\) there are at most \(4\) distinct values of \(f_{k, i, j}\). So let us partition the set of squares on the grid into the following four:

• \(S_{0, 0} = \{ (x_1, y_1) \}\)
• \(S_{0, 1} = \{ (x_1, y) \, \vert \, y \neq y_1 \}\)
• \(S_{1, 0} = \{ (x, y_1) \, \vert \, x \neq x_1 \}\)
• \(S_{1, 1} = \{ (x, y) \, \vert \, x \neq x_1 \land y \neq y_1 \}\) It can be proved by induction that, for a fixed \(k\), all values \(f_{k, i, j}\) of the squares in each set \(S_{a, b}\) is all equal.

Let \(g_{k, a, b} = \sum_{(i, j) \in S_{a, b}} f(k, i, j)\). Then, for \(k > 0\) the following relations hold.

• \(g_{k, 0, 0} = g_{k - 1, 0, 1} + g_{k - 1, 1, 0}\)
• \(g_{k, 0, 1} = (Y - 1) g_{k - 1, 0, 0} + (Y - 2) g_{k - 1, 0, 1} + g_{k - 1, 1, 1}\)
• \(g_{k, 1, 0} = (X - 1) g_{k - 1, 0, 0} + (X - 2) g_{k - 1, 1, 0} + g_{k - 1, 1, 1}\)
• \(g_{k, 1, 1} = (X - 1)g_{k - 1, 0, 1} + (Y - 1)g_{k - 1, 1, 0} + (X - 2 + Y - 2) g_{k - 1, 1, 1}\)

If \((x_2, y_2)\) belongs to some set \(S_{a, b}\), then the answer is \(\displaystyle \frac{g_{K, a, b}}{ |S_{a, b}|}\). The time complexity is \(O(K)\).

Sample code (C++) :

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

using Fp = modint998244353;

int main() {
    int h, w, k;
    cin >> h >> w >> k;
    array<Fp, 2> row = {1, h - 1}, col = {1, w - 1};
    array<array<Fp, 2>, 2> dp = {};
    dp[0][0] = 1;
    while (k--) {
        array<array<Fp, 2>, 2> next = {};
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                for (int a = 0; a < 2; ++a) {
                    next[i][a] += dp[i][j] * (col[a] - (j == a));
                    next[a][j] += dp[i][j] * (row[a] - (i == a));
                }
            }
        }
        dp = move(next);
    }
    int x1, y1, x2, y2;
    cin >> x1 >> y1 >> x2 >> y2;
    int i = (x1 != x2), j = (y1 != y2);
    cout << (dp[i][j] / row[i] / col[j]).val() << '\n';
    return 0;
}

Bonus: one can use matrix exponentiation to solve it in a total of \(O(\log K)\) time.