学了学 EI 的做法……
总算是看懂了。
设 \(n\) 对情侣都不和睦的方案数为 \(f_n\)，则恰有 \(k\) 对情侣和睦的方案数为 \[{\binom nk}^2 f_{n-k} k!2^k\]
两个组合数分别是从 \(n\) 对情侣里选 \(k\) 对，从 \(n\) 排里选 \(k\) 排；\(k!\) 是这 \(k\) 对情侣的排列方案；\(2^k\) 是情侣两个人互相换位置。
那么枚举和睦情侣的对数有 \[ \sum\limits_{k=0}^n {\binom nk}^2 f_{n-k} k!2^k = (2n)! \]
众所周知 EGF 的卷积会附带一个组合数，而上式中除了组合数以外的都显然是卷积的形式。
于是可以如法炮制 EGF，设计一种新的生成函数： \[
F(x) = \sum\limits_{n=0}^{\infty} f_n \frac{x^n}{(n!)^2}
\]
这样的话卷积就可以有两个组合数了。
那么我们考虑上式中除了 \(F\) 以外的两个数列 \[ \begin{aligned} \sum\limits_{n=0}^{\infty} n!2^n\frac{x^n}{(n!)^2} &= \sum\limits_{n=0}^{\infty} \frac{2^nx^n}{n!} = \mathrm e^{2x} \\ \sum\limits_{n=0}^{\infty} (2n)!\frac{x^n}{(n!)^2} &= \sum\limits_{n=0}^{\infty} \binom{-\frac12}n (-4)^n x^n = (1-4x)^{-1/2} \end{aligned} \]
即 \[ \begin{aligned} \mathrm e^{2x} F(x) &= (1-4x)^{-1/2} \\ F(x) &= \mathrm e^{-2x}(1-4x)^{-1/2} \end{aligned} \]
求个导看看 \[ \begin{aligned} F'(x) &= 8x \mathrm e^{-2x}(1-4x)^{-3/2} \\ &= \frac{8x}{1-4x} F(x) \\ (1-4x)F'(x) &= 8xF(x) \\ F'(x) &= 4xF'(x) + 8xF(x) \\ \end{aligned} \]
又因为 \[ \begin{aligned} \left(\sum\limits_{n=0}^{\infty} f_n\frac{x^n}{(n!)^2}\right)' &= \sum\limits_{n=0}^{\infty} \left(f_n\frac{x^n}{(n!)^2}\right)' \\ &= \sum\limits_{n=1}^{\infty} f_n\frac{nx^{n-1}}{(n!)^2} \\ &= \sum\limits_{n=0}^{\infty} (n+1)^{-1} f_{n+1}\frac{x^n}{(n!)^2} \\ xF'(x) &= \sum\limits_{n=0}^{\infty} (n+1)^{-1}f_{n+1}\frac{x^{n+1}}{(n!)^2} \\ &= \sum\limits_{n=1}^{\infty} nf_n\frac{x^n}{(n!)^2} \\ xF(x) &= \sum\limits_{n=0}^{\infty} f_n\frac{x^{n+1}}{(n!)^2} \\ &= \sum\limits_{n=1}^{\infty} n^2 f_{n-1}\frac{x^n}{(n!)^2} \end{aligned} \]
于是 \[ \begin{aligned} (n+1)^{-1} f_{n+1} &= 4nf_n + 8n^2f_{n-1} \\ f_{n+1} &= 4n(n+1)f_n + 8n^2(n+1)f_{n-1} \qquad(n \ge 1) \end{aligned} \]
初值为 \(f_0=1,f_1=0\)。
代码： 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
using namespace std;
const int N = 5e6;
const int mod = 998244353;
int T,n,k;
int fac[N + 5],ifac[N + 5];
int f[N + 5];
inline int fpow(int a,int b)
{
int ret = 1;
for(;b;b >>= 1)
(b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
return ret;
}
int main()
{
f[0] = fac[0] = 1;
for(register int i = 2;i <= N;++i)
f[i] = (4LL * (i - 1) % mod * i % mod * f[i - 1] % mod + 8LL * (i - 1) % mod * (i - 1) % mod * i % mod * f[i - 2] % mod) % mod;
for(register int i = 1;i <= N;++i)
fac[i] = (long long)fac[i - 1] * i % mod;
ifac[N] = fpow(fac[N],mod - 2);
for(register int i = N;i;--i)
ifac[i - 1] = (long long)ifac[i] * i % mod;
scanf("%d",&T);
for(;T;--T)
scanf("%d%d",&n,&k),printf("%lld\n",(long long)fac[n] * fac[n] % mod * ifac[k] % mod * ifac[n - k] % mod * ifac[n - k] % mod * fpow(2,k) % mod * f[n - k] % mod);
}