LOJ3058 [HNOI2019]白兔之舞

过程比较繁琐，知识比较综合。

题意：在 $L + 1$ 行 $n$ 列的图形上从 $(0, x)$ 随机游走，要求只能向下走，从 $v_1$ 列走到 $v_2$ 列的方案数为 $W_{v_1, v_2}$（与行无关），只能在 $y$ 列停止。给定 $k$，求所有 $t\in[0, k - 1]$，求停留的行的编号满足 $\bmod k = t$ 的所有方案和，答案对 $p$ 取模。满足 $k\leq 2 ^ {16}$，$L\leq 10 ^ 8$，$10 ^ 8 < p < 2 ^ {30}$，$n\leq 3$，$k|p - 1$，$p$ 为质数。

仔细观察发现如果我们确定步数 $k$，行和列的答案是无关的，只需要最后乘起来即可。对于行来说，我们假设多给一步，强制走到 $L + 1$ 的位置，于是插板法容易得到方案数为 $\binom Lk$。对于列来说，这就是一个显然的矩阵乘法，为 $W ^ k[x][y]$。

于是可以得到：

$$ans_t = \sum_{i = 0} ^ L [i\bmod k = t] W ^ i[x][y]\binom Li$$

然后看到 $[i\bmod k = t]$，果断单位根反演，于是一波变换：

$$\begin{aligned} ans_t = & \sum_{i = 0} ^ L [i\bmod k = t] W ^ i[x][y]\binom Li\\ =& \dfrac 1k \sum_{i = 0} ^ L W ^ i [x][y]\binom Li \sum_{j = 0} ^ {k - 1} \omega_{k} ^ {ij - tj}\\ =& \dfrac 1k \sum_{j = 0} ^ {k - 1} \omega_k ^ {-tj}\sum_{i = 0} ^ L W ^ i[x][y] \omega_k ^ {ij}\binom Li\\ =& \dfrac 1k \sum_{j = 0} ^ {k - 1} \omega_k ^ {-tj} (\omega_k ^ j W + I) ^ L[x][y] \end{aligned}$$

最后一步我们在 单位根反演 的 T2 中已经见过矩阵合并二项式定理了，这里不再赘述。注意这里单位根需要我们先求出原根计算。

假设后面只关于 $j$ 的一坨记作 $f(j)$，其为 $\sum_{j = 0} ^ {k - 1} \omega_k ^ {-tj} (\omega_k ^ j W + I) ^ L[x][y]$。单个 $f(j)$ 容易在 $\log L$ 的时间完成。

于是我们得到：

$$k\times ans_t = \sum_{j = 0} ^ {k - 1} \omega_k ^ {-tj} f(j)$$

看样子这是一个卷积，但是指数上有乘积，不好拆开。

考虑如下两个变换：

$$\begin{aligned} (x + y) ^ 2 - x ^ 2 - y ^ 2 =& 2xy\\ \binom{x + y}2 - \binom x2 - \binom y2 =& xy \end{aligned}$$

两个式子展开一下即可证明。我们惊喜的发现，这个恰好把 $xy$ 化为了 $(x + y)$，这样看样子就可以卷积了！

套用第一个式子：

$$k\times ans_t = \omega_k ^ {\frac{t ^ 2}2}\sum_{j = 0} ^ {k - 1} f(j)\omega_k ^ {\frac{j ^ 2}{2}}\omega_k ^ {-\frac{(j + t) ^ 2}2}$$

但是除以 2 不好处理，因为我们可能没法开根。

考虑套用第二个式子：

$$\begin{aligned} k\times ans_t =& \sum_{j = 0} ^ {k - 1} \omega_k ^ {-tj} f(j)\\ =& \sum_{j = 0} ^ {k - 1} \omega_k ^ {\binom t2+ \binom j2 - \binom{t + j}2}f(j)\\ =& \omega_k ^ {\binom t2}\sum_{j = 0} ^ {k - 1} f(j)\omega_k ^ {\binom j2}\omega_k ^ {-\binom{j + t}{2}} \end{aligned}$$

显然这是一个差卷积，直接做即可，时间复杂度 $O(k\log k)$，可以通过。

另外，这后半部分拆开指数上的乘积其实就是一个叫 BlueStein 的算法。容易发现这个算法其实解决了任意长度的 FFT，因为 FFT 和这个算法都是计算 $n$ 个单位根位置的点值的，只不过 FFT 的 $n$ 是 2 的次幂。更详细的信息看 这篇 Blog。

int main()
{
	Matrix a;
	int k, L, x, y, G;
	std::cin >> n >> k >> L >> x >> y >> Mod, -- x, -- y;
	for (int i = 0; i < n; ++ i)
		for (int j = 0; j < n; ++ j) std::cin >> a[i][j];
	w = qpow(G = findrt(), (Mod - 1) / k);
	poly f(k);
	for (int l = 0; l < k; ++ l)
	{
		static Matrix tmp;
		int cur = qpow(w, l);
		for (int i = 0; i < n; ++ i)
			for (int j = 0; j < n; ++ j)
				tmp[i][j] = ((LL) a[i][j] * cur + (i == j)) % Mod;
		f[l] = qpow(tmp, L)[x][y];
	}
	for (int i = 0; i < k; ++ i)
		f[i] = (LL) f[i] * qpow(w, (i - 1LL) * i / 2 % (Mod - 1)) % Mod;
	std::reverse(f.begin(), f.end());
	poly g(2 * k);
	for (int i = 0; i < 2 * k; ++ i)
		g[i] = qpow(w, (Mod - 1 - i * (i - 1LL) / 2 % (Mod - 1)) % (Mod - 1));
	f = f * g, f = std::vector<int>(f.data() + k - 1, f.data() + 2 * k - 1);
	for (int i = 0; i < k; ++ i)
		f[i] = (LL) f[i] * qpow(w, (i - 1LL) * i / 2 % (Mod - 1)) % Mod * qpow(k) % Mod;
	for (int &x : f) printf("%d\n", x);
	return 0;
}