Consider a DP where \(dp[i][j] =\) the number of ways to determine the first \(i\) terms of \(P\) and \(Q\) such that there are \(j\) integers that have appeared in \(P\) but not in \(Q\).

Suppose \(dp[N-M][0],dp[N-M][1],\dots,dp[N-M][M]\) have been computed. The number of ways to determine the remaining \(M\) terms given \(dp[N-M][k]\) equals the number of pairs of permutations of \((1,2,\dots,M)\) such that \(1,2,\dots,M-k\) do not appear in the same position in both; by inclusion-exclusion, this is \(\sum_{i=0}^{M-k} \binom{M-k}{i} (-1)^i M! (M-i)!\) ways. Let this value be \(A_k\). This can be computed via convolution in \(\mathrm{O}(M \log M)\).

Thus, it suffices to compute \(dp[N-M][0],dp[N-M][1],\dots,dp[N-M][M]\). Define \(f_j(x) = \sum_i dp[i][j] x^i\). Actually, the transition differs before and after \(dp[N-M]\) depending on whether the number of usable values increases by \(1\) or not, but we will use the transition where the number of usable values keeps increasing by \(1\) even after \(dp[N-M]\).

From the DP transition, we have \(f_0(x) = 1 + xf_1(x)\), and for \(1 \le i \le M\), \(f_i(x) = x((M-i+1)(M-i+2)f_{i-1}(x) + 2i(M-i+1)f_i(x) + (i+1)^2f_{i+1}(x))\).

Since \(f_i(x)\) is divisible by \(x^i\), substitute \(g_i(x) = \frac{f_i(x)}{x^i}\). Then, we get the equation \((M-i+1)(M-i+2)g_{i-1}(x)+ (2i(M-i+1)x-1)g_i(x) + (i+1)^2x^2g_{i+1}(x) = 0\).

Consider solving this recurrence from the end.

Let:

\[h_i = \begin{pmatrix} \frac{1-2i(M-i+1)x}{(M-i+1)(M-i+2)} & \frac{-(i+1)^2x^2}{(M-i+1)(M-i+2)} \\ 1 & 0 \end{pmatrix}\]

Then \(\begin{pmatrix} g_{i-1}(x) \\ g_i(x) \end{pmatrix} = h_i \begin{pmatrix} g_{i}(x) \\ g_{i+1}(x) \end{pmatrix}\) holds. That is, \(\begin{pmatrix} g_{0}(x) \\ g_1(x) \end{pmatrix} = h_1 h_2 ... h_M \begin{pmatrix} g_{M}(x) \\ g_{M+1}(x) \end{pmatrix}\). Here, let \(H = h_1 h_2 ... h_M\).

Noting \(g_{M+1}(x) = 0\), we get \(g_0(x) = H_{0,0} g_{M}(x),g_1(x) = H_{1,0} g_M(x)\). Also, using \(g_0(x) = 1 + x^2 g_1(x)\), we get \((H_{0,0} - x^2H_{1,0})g_M(x) = 1\). Thus, \(g_M(x) = \frac{1}{H_{0,0} - x^2H_{1,0}}\). The constant term of \(g_M(x)\) is \((M+1)!M!\), which under the given constraints is nonzero modulo \(998244353\), so \(g_M(x)\) can be computed using the inverse of a formal power series.

The answer is \(\sum_{i=0}^{M} A_i dp[N-M][i] = [x^{N-M}]\sum_{i=0}^{M} A_i f_i(x) = [x^{N-M}] \sum_{i=0}^{M} A_i x^i g_i(x)\). To find this polynomial, let:

\[h'_i = \begin{pmatrix} \frac{1-2i(M-i+1)x}{(M-i+1)(M-i+2)} & \frac{-(i+1)^2x^2}{(M-i+1)(M-i+2)} & 0 \\ 1 & 0 & 0 \\ A_i & 0 & x \\ \end{pmatrix}\]

and let \(H' = h'_1 h'_2 ... h'_M\); then it suffices to compute \([x^{N-M}] xH'_{2,0} g_M(x) + A_0 g_0(x)\). Once \(H'\) is known, this value can be computed in \(\mathrm{O}(M \log N \log M)\) using the Bostan–Mori algorithm.

Thus, by computing the matrix product part via divide and conquer, the answer can be found in \(\mathrm{O}(M \log M (\log M + \log N))\). Also, since the top-left \(2 \times 2\) block of \(H'\) coincides with \(H\), it suffices to actually compute only \(H'\). (This does not change the order, but improves the constant factor.)