构造 DFT ：某少女附中的体育课 题解

构造 DFT ：某少女附中的体育课 题解

题意

神仙题。一题能当十题写。

**一句话题意：**给出一个下标系统 $A$，求序列 $P$ 在该下标系统下的 $N$ 维卷积的 $K+1$ 次幂。

先来解释一下什么意思。

通常来讲，卷积是一个以下的形式 $$ c_r = \sum_{p \circ q = r} a_pb_q $$ 当 $\circ$ 为 $+$ 的时候，就是咱们熟悉的多项式乘法。

那若是是位运算呢？就是 $\text{FWT}$ 干的事情了。

那若是是其余的东西呢？这个题就是给了你一个 $m \times m$ 的矩阵 $A$ ，表示一个运算表。

而后你有一个长度为 $m^n$ 的序列 $P$ ，每一个下标是一个 $n$ 位 $m$ 进制数，下标之间的运算就是按位进行矩阵 $A$ 下的运算。你须要在这个意义下算出 $P$ 的 $K+1$ 次幂。

为了方便，如下设 $N = m^n$ 。

算法一

直接暴力计算卷积便可。

预计得分： $7$。

算法二

咱们注意到下标运算系统其实就是同或。直接上同或 $\text{FWT}$ 就行。

什么，你不会同或 $\text{FWT}$？没事，我也不会。你只须要会异或 $\text{FWT}$ 就好了。

咱们注意到两个数的同或其实就是异或值的取反，因此作完异或 $\text{FWT}$ 以后 std::reverse一下就好了。

可是这是作一遍卷积，两遍卷积就至关于又 std::reverse 回来了。因此 $K+1$ 为奇数的时候才须要 std::reverse。

而后题目保证了 $A$ 是关于主对角线对称的（交换律），因而 $m=2$ 的状况其实只有四种 $A$ ，分别对应与，或，同或，异或，直接上四种 $\text{FWT}$ 就好了。

预计得分：$25$ 。

前二十五分并不须要什么冷门知识。

算法三

剩下的就有点意思了。

Subtask4中的 $A$ 表示按位取 max ，类比 $\text{FWT}$ ，咱们只须要将 $0$ 的位置的贡献加到 $1$ 上，而后再将 $1$ 的位置的贡献加到 $2$ 上，写起来和 $ \text{OR FWT}$ 差很少，就是个高维前缀和。

预计得分： $34$ 。

算法四

为何写起来和 $\text{FWT}$ 差很少呢？

vfk的论文里说，$\text{FWT}$ 本质上是一个高维的 $\text{DFT}$ 。而什么又是高维 $\text{DFT}$ 呢？

注意， $\text{FFT}$ 是快速求 $\text{DFT}$ 的算法，而 $\text{DFT}$ 是一个变换，不要搞混。

先来回顾一下 $\text{DFT}$ 。

咱们知道 $\text{DFT}$ 是求出了一个多项式在单位根处的点值表示。为何是主单位根呢？由于 $\text{DFT}$ 最核心的原理是求和引理 $$ \frac 1 n \sum_{k=0}^{n-1}\omega^{vk} = [v \mod n = 0] $$ 其中 $v$ 是定值。这个式子能够用等比数列求和公式来证实。

$v \mod n \not = 0$ 时，有 $w^v \not= 1$： $$ \begin {align} &\frac 1 n \sum_{k=0}^{n-1}\omega^{vk}\ =&\frac 1 n \frac {1-\omega^{vn}} {1-\omega^v} \ =& 0 \end {align} $$ $v \mod n = 0$ 时，有 $\omega^v = 1$： $$ \begin {align} &\frac 1 n \sum_{k=0}^{n-1}\omega^{vk}\ =&1 \end {align} $$ 回顾一下咱们多项式乘法卷积的式子： $$ \begin {align} c_r &= \sum_{p + q = k} a_pb_q\ &= \sum_{p,q}[p+q=r]a_pb_q \end {align} $$ 咱们先悄咪咪的将 $[p+q=r]$ 变成 $[(p+q)\mod n=r]$

因而咱们就能够用求和引理了！ $$ \begin {align} &[(p+q)\mod n=r]\ =&[(p+q-r)\mod n=0] \ =& \frac 1 n \sum_{k=0}^{n-1}\omega^{(p+q-r)k} \ =& \frac 1 n \sum_{k=0}^{n-1}\omega^{-rk}\omega^{pk}\omega^{qk} \end {align} $$

而后咱们继续化式子 $$ \begin {align} c_r &= \sum_{p + q = k} a_pb_q\ &= \sum_{p,q}[p+q=r]a_pb_q \ &= \sum_{p,q}[(p+q)\mod n=r]a_pb_q \ &= \sum_{p,q}\frac 1 n \sum_{k=0}^{n-1}\omega^{-rk}\omega^{pk}\omega^{qk}a_pb_q \ &= \frac 1 n \sum_{k=0}^{n-1}\omega^{-rk}\sum_{p,q}\omega^{pk}a_p\omega^{qk}b_q \ &= \frac 1 n \sum_{k=0}^{n-1}\omega^{-rk}\sum_{p}\omega^{pk}a_p\sum_{q}\omega^{qk}b_q \end {align} $$ 因而咱们只须要定义 $\text{DFT}$ 为 $$ \hat a_r = \sum_{q = 1}^{n-1}\omega^{rq}a_q $$ 定义 $\text{IDFT}$ 为 $$ a_r = \frac 1 n \sum_{q = 1}^{n-1}\omega^{-rq}\hat a_q $$ 因而咱们就有了一个卷积性变换 $$ \text{DFT}(fg) = \text{DFT}(f) \cdot \text{DFT}(g) \ fg = \text{IDFT}(\text{DFT}(f) \cdot \text{DFT}(g)) $$

$O(n^2)\to O(n\log n)$ 。

慢着，咱们刚刚悄咪咪的将 $[p+q=r]$ 变成 $[(p+q)\mod n=r]$ 真的没问题嘛？

固然有啊！由于原理就是这个求和引理，因此咱们的 $\text{DFT}$ 是算的是 $\mod n$ 意义下的循环卷积。因而乎，当两个下标之和大于等于 $n$ 的时候，他们的值的贡献会加到前面去。因此咱们得作一个长度为 $2n$ 的 $\text{DFT}$，才能获得咱们想要的结果。

而什么又是高维 $\text{DFT}$ 呢？

考虑高维前缀和，咱们能够对每一个方向依次作一遍一维的前缀和。

高维 $\text{DFT}$ 就把求一维前缀和换成一维变换。

高维 $\text{DFT}$ 至关于对下标进行按位分治，而后分层进行 $\text{DFT}$ 。就是说，枚举当前的位，而后对下标除了这一位以外都相同的每组数作一遍 $\text{DFT}$ 就行了。

而异或 $\text{FWT}$ 就是作了一个每维长度为 $2$ 的高维 $\text{DFT}$ 。

为何呢？

由于异或运算和 $\mod 2$ 意义下的加法等价，每一维本质上就是一个 $\text{DFT}$ 。想一想看，是否是？

虽然代码和 $\text{FFT}$ 几乎如出一辙，可是原理是不一样的。咱们刚才说 $\text{FFT}$ 是个算法，它实际上是在用分治法快速算 $\text{DFT}$ 。而 $\text{FWT}$ 实际上每维是暴力算出来的 $\text{DFT}$ 。

长度为 $2$ 的 $\text{DFT}$ ，须要 $2$ 次单位根。不就是 $-1$ 嘛！

而后 $$ \hat a_r = \sum_{q = 1}^{n-1}\omega^{rq}a_q $$ 然而咱们只有两个数，因此 $$ \hat a_0 = a_0 + a_1 \ \hat a_1 = a_0 - a_1 $$ 这样也能解释为何回来的时候要除 $2$ 。由于是在作 $\text{IDFT}$ 。

如下全部涉及到的 $\text{DFT}$ 都是高维 $\text{DFT}$ ，可是实现方法是同样的，你若是能构造出一维上的 $\text{DFT}$ ，那么你只须要按位分治，而后去掉这一位进行分组，将去掉这一位以后下标同样的数分到同一组里，而后对每组分别作一维的 $\text{DFT}$ 就行了。

回到这个题上。

Subtask5 给了一个奇怪的 $A$ 。经过人类智慧稍微构造一下，咱们就能找到这样一个 $\text{DFT}$ ： $$ \hat a_0 = a_2 \ \hat a_1 = a_0 + a_1 + a_2\ \hat a_2 = a_0-a_1+a_2 $$ 而后继续作就好辣！

预计得分： $43$ 。

算法五

你不必知道什么是循环群。你只须要知道它和模 $m$ 意义下的加法是同构的就行了。

不过咱们仍是得须要知道什么是阶和生成元。

幂的定义题面里给了。

这个阶和数论里的阶是差很少的，若 $i^{j+1} = i$ 即 $i^j = \epsilon $ ，那么 $i$ 的阶就是 $j$ 。记为 $ord(i)$ 。

生成元和原根差很少，就是一个数 $i$ 它的 $[0,ord(i))$ 次幂能遍历全部元素。其实就是那个阶最大的。

而后 $\mod m$ 意义下的加法的幂其实就是乘法，单位元 $\epsilon$ 是 $0$，生成元是全部与 $m$ 互质的数，例如 $1$ ，它的阶为 $m$ 。

而后咱们暴力求出这个 $[0,m)$ 里每一个数在 $A$ 下的阶，里面确定有个元素的阶等于 $m$ 。咱们把它当成 $1$ ，而后让它本身与本身运算，得出其余的数的映射关系。而后咱们就把下标运算转换成了 $\mod m$ 的加法，直接上 $\text{DFT}$ 而后转换回来就好辣！

预计得分： $53$ 。

算法六

后面的部分分我不会。就算看懂了题解我也不会实现。

LCA 的标程密密麻麻的看不懂，论文里用的伪代码，根本无从下手。

（P.S. LCA 的提交之因此在 LOJ 少女附中的榜里是最慢的是由于他交的实际上是个暴力，并且还用 std::cin std::cout 读写 $5e5$ 级别的数据。）

好在这个题暴力可过。因而咱们就能愉快的水过去了。

$\text{DFT}$ 本质上是一个“线性变换”。

什么意思？

至关于你给那个数列乘了一个矩阵，每一个数只会变成一堆其余的数再乘上某个系数后的和。

咱们找到那个矩阵不就好了？

咱们多项式乘法的 $\text{DFT}$ 和 $\text{FWT}$ 的 $\text{DFT}$ 的矩阵是单位根的范德蒙德矩阵。

因而咱们大胆猜测不用证实咱们这个矩阵里放的也是单位根的几回幂。

设这个矩阵为 $T$ 。

变换就能写成 $$ \hat a_r = \sum_{p=0}^{n-1}T_{rp}a_p $$

咱们化一化式子，看看它须要知足什么性质。 $$ \begin {align} \hat c_r &= \hat a_r \cdot \hat b_r \ \sum_{t=0}^{n-1}T_{rt}c_t &= \sum_{p=0}^{n-1}T_{rp}a_p \sum_{q=0}^{n-1}T_{rq}b_q \ \sum_{t=0}^{n-1}T_{rt}\sum_{p,q}a_pb_q[p \circ q=t] &= \sum_{p=0}^{n-1}T_{rp}a_p \sum_{q=0}^{n-1}T_{rq}b_q \ \sum_{p,q}a_pb_qT_{r(p \circ q)} &= \sum_{p=0}^{n-1}T_{rp}a_p \sum_{q=0}^{n-1}T_{rq}b_q \ \sum_{p,q}a_pb_qT_{r(p \circ q)} &= \sum_{p=0}^{n-1}\sum_{q=0}^{n-1}T_{rp}a_p T_{rq}b_q \ \sum_{p,q}a_pb_qT_{r(p \circ q)} &= \sum_{p,q}a_pb_q T_{rp}T_{rq} \ T_{r(p \circ q)} &= T_{rp}T_{rq} \end {align} $$

因而 $T$ 一定每一行内都知足 $T_{r(p \circ q)} = T_{rp}T_{rq}$。

而后根据题面里说的“循环率”，元素 $i$ 的 $ord(i)+1$ 次幂得是 $i$ 本身。因而 $T_{rq}^{ord(q)+1} = T_{rq}$ 。

因此根据这个性质咱们 $m^m$ 枚举，在每个位置 $i$ 填上 $ord(i)$ 次单位根的 $[0,ord(i)]$ 次幂就行了。而后填一个数检查一下是否当前行知足 $T_{r(p \circ q)} = T_{rp}T_{rq}$，填满一行以后存到矩阵里。

搜出来这个矩阵以后，给它求一个逆，就是逆变换的矩阵辣！正确性显然。

啥？你不会矩阵求逆？出门右转洛谷模板区。挺好写的。就一个高斯消元。

然而会有一个问题，若是求出来的矩阵没有逆怎么办？

可是根据 LCA 的证实，好像搜出来的矩阵各行一定都线性无关。就是必定有逆。

搜出来矩阵和逆矩阵直接 $\text{(I)DFT}$ 就好了。随便搜搜就能过。跑的挺快的。

代码超级好写。

顺便提一下，正解使用的是 Subtask7 和 Subtask8的方法科学构造出那个变换矩阵。想学的能够学一学。代码不长。

预计得分：$100$。

放一个部分分比较全的代码吧，能够参考一下。

#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
#include <vector>
inline int read() {
	register int ret, cc;
	while (!isdigit(cc = getchar())){}ret = cc-48;
	while ( isdigit(cc = getchar()))  ret = cc-48+ret*10;
	return ret;
}
const int MOD = 232792561;
const int G = 71;
inline int add(int a, int b) { return (a += b) >= MOD ? a -= MOD : a; }
inline int mul(int a, int b) { return 1ll * a * b % MOD; }
inline int qpow(int a, int p) {
	int ret = 1;
	for ( ; p; a = mul(a, a), p >>= 1)
		if (p & 1) ret = mul(a, ret);
	return ret;
}
int N, M, X, U;
long long K;


namespace SubTask1 {
	const int MAXN = 30;
	const int MAXU = 1200;
	int A[MAXN][MAXN];
	int P[MAXU];
	std::vector<int> state[MAXU];
	int trans[MAXU][MAXU];
	int f[MAXN][MAXU];

	inline std::vector<int> decode(int S) {
		std::vector<int> ret(N);
		for (int i = 0; i < N; ++i) ret[i] = S % M, S /= M;
		return ret;
	}
	inline int encode(const std::vector<int>& S) {
		int ret = 0;
		for (int i = N - 1; i >= 0; --i) ret = ret * M + S[i];
		return ret;
	}
	inline std::vector<int> transform(const std::vector<int>& lhs, const std::vector<int>& rhs) {
		std::vector<int> ret(N);
		for (int i = 0; i < N; ++i) ret[i] = A[lhs[i]][rhs[i]];
		return ret;
	}
	void Main() {
		for (int i = 0; i < M; ++i) for (int j = 0; j < M; ++j)	A[i][j] = read();
		for (int i = 0; i < U; ++i) P[i] = read();
		for (int i = 0; i < U; ++i) state[i] = decode(i);
		for (int i = 0; i < U; ++i) for (int j = 0; j < U; ++j) trans[i][j] = encode(transform(state[i], state[j]));
		for (int i = 0; i < U; ++i) f[0][i] = P[i];
		for (int i = 0; i < K; ++i)
			for (int j = 0; j < U; ++j)
				for (int k = 0; k < U; ++k)
					f[i+1][trans[j][k]] = add(f[i+1][trans[j][k]], mul(f[i][j], P[k]));
		for (int i = 0; i < U; ++i) printf("%d\n", f[K][i]);
	}
}
namespace SubTask23 {
	const int MAXU = 1 << 20 | 1;
	int P[MAXU];
	inline void FWT_OR(int* a, int n, int opt) {
		for (int i = 1; i < n; i <<= 1) 
			for (int j = 0, p = i << 1; j < n; j += p)
				for (int k = 0; k < i; ++k)
					(a[j + k + i] += (opt * a[j + k] + MOD) % MOD) %= MOD;
	}
	inline void FWT_AND(int* a, int n, int opt) {
		for (int i = 1; i < n; i <<= 1)
			for (int j = 0, p = i << 1; j < n; j += p)
				for (int k = 0; k < i; ++k)
					(a[j + k] += (opt * a[j + k + i] + MOD) % MOD) %= MOD;
	}
	inline void FWT_XOR(int* a, int n, int opt) {
		for (int i = 1; i < n; i <<= 1)
			for (int j = 0, p = i << 1; j < n; j += p)
				for (int k = 0; k < i; ++k) {
					int x = a[j + k], y = a[j + k + i];
					a[j + k] = (x + y) % MOD;
					a[j + k + i] = (x - y + MOD) % MOD;
				}
		if (opt == -1) {
			int inv = qpow(n, MOD - 2);
			for (int i = 0; i < n; ++i) a[i] = 1ll * a[i] * inv % MOD;
		}
	}
	inline void FWT_NXOR(int* a, int n, int opt) {
		FWT_XOR(a, n, opt);
		if (opt == -1) std::reverse(a, a + n);
	}
	void Main() {
		int type = 0;
		for (int i = 0; i < 4; ++i) type = type << 1 | read();
		for (int i = 0; i < U; ++i) P[i] = read();
		if ((~K&1) && type == 9) type ^= 15;
		switch (type) {
			case 9: FWT_NXOR(P, U, 1); break;
			case 6: FWT_XOR(P, U, 1);  break;
			case 7: FWT_OR(P, U, 1);   break;
			case 1: FWT_AND(P, U, 1);  break;
		}
		for (int i = 0; i < U; ++i) P[i] = qpow(P[i], (K + 1) % (MOD - 1));
		switch (type) {
			case 9: FWT_NXOR(P, U, -1); break;
			case 6: FWT_XOR(P, U, -1);  break;
			case 7: FWT_OR(P, U, -1);   break;
			case 1: FWT_AND(P, U, -1);  break;
		}
		for (int i = 0; i < U; ++i) printf("%d\n", P[i]);
	}
}
int P[500010];
int A[30][30];
int W[30][30];
std::vector<int> group;
const std::vector<int> subtask4({0, 1, 2, 1, 1, 2, 2, 2, 2});
const std::vector<int> subtask5({0, 1, 0, 1, 0, 1, 0, 1, 2});
namespace SubTask4 {
	void Main() {
		for (int i = 1; i <= U; i *= M) {
			for (int j = 0; j < U; ++j) {
				int p = j % i, q = j / i;
				if (q % 3) P[j] = add(P[j], P[(q-1)*i+p]);
			}
		}
		for (int i = 0; i < U; ++i) P[i] = qpow(P[i], (K+1)%(MOD-1));
		for (int i = 1; i <= U; i *= M) {
			for (int j = U-1; j >= 0; --j) {
				int p = j % i, q = j / i;
				if (q % 3) P[j] = add(P[j], MOD-P[(q-1)*i+p]); 
			}
		}
		for (int i = 0; i < U; ++i) printf("%d\n", P[i]);
	}
}
namespace SubTask5 {
	void Main() {
		const int inv2 = (MOD+1)/2;
		for (int i = 1; i < U; i *= 3) {
			for (int j = 0, p = i*3; j < U; j += p) {
				for (int k = 0; k < i; ++k) {
					int x = P[j+k], y = P[j+k+i], z = P[j+k+2*i];
					P[j+k] = z;
					P[j+k+i] = add(x, add(y, z));
					P[j+k+2*i] = add(x, add(MOD-y, z));
				}
			}
		}
		for (int i = 0; i < U; ++i) P[i] = qpow(P[i], (K+1)%(MOD-1));
		for (int i = 1; i < U; i *= 3) {
			for (int j = 0, p = i*3; j < U; j += p) {
				for (int k = 0; k < i; ++k) {
					int x = P[j+k], y = P[j+k+i], z = P[j+k+2*i];
					P[j+k] = add(mul(inv2, add(y, z)), MOD-x);
					P[j+k+i] = mul(inv2, add(y, MOD-z));
					P[j+k+2*i] = x;
				}
			}
		}
		for (int i = 0; i < U; ++i) printf("%d\n", P[i]);
	}
}
namespace SubTask6 {
	inline int calc(int a, int b) { return A[a][b]; }	
	int tr[10];
	inline int trans(int x) {
		int ret = 0;
		for (int i = 1; i < U; i *= M) 
			ret += tr[x / i % M] * i;
		return ret;
	}
	inline int getroot() {
		int root = -1;
		for (int i = 0; i < M; ++i) {
			int cur = A[i][i];
			int rank = 1;
			while (cur != i) {
				++rank, cur = A[cur][i];
			}
			if (rank == M) {
				root = i;
				break;
			}
		}
		return root;
	}
	int PP[500010];
	int tmp[10];
	void Main() {
		int root = getroot();
		for (int i = 1, cur = root; i <= M; ++i, cur = A[cur][root]) tr[cur] = i % M;
		for (int i = 0; i < U; ++i) PP[trans(i)] = P[i];
		for (int i = 1; i < U; i *= M) {
			for (int j = 0, p = i * M; j < U; j += p) {
				for (int k = 0; k < i; ++k) {
					for (int r = 0; r < M; ++r) tmp[r] = PP[j+k+r*i];
					int w = 1, wm = W[M][1];
					for (int r = 0; r < M; ++r, w = mul(w, wm)) {
						int sum = 0;
						for (int t = M-1; ~t; --t) sum = add(mul(sum, w), tmp[t]);
						PP[j+k+r*i] = sum;
					}
				}
			}
		}
		for (int i = 0; i < U; ++i) PP[i] = qpow(PP[i], (K+1) % (MOD-1));
		int invm = qpow(M, MOD-2);
		for (int i = 1; i < U; i *= M) {
			for (int j = 0, p = i * M; j < U; j += p) {
				for (int k = 0; k < i; ++k) {
					for (int r = 0; r < M; ++r) tmp[r] = PP[j+k+r*i];
					int w = 1, wm = qpow(W[M][1], MOD-2);
					for (int r = 0; r < M; ++r, w = mul(w, wm)) {
						int sum = 0;
						for (int t = M-1; ~t; --t) sum = add(mul(sum, w), tmp[t]);
						PP[j+k+r*i] = mul(sum, invm);
					}
				}
			}
		}
		for (int i = 0; i < U; ++i) printf("%d\n", PP[trans(i)]);
	}
}
namespace Std_Dfs {
	bool vis[30];
	int ord[30];
	inline void getord() {
		for (int i = 0; i < M; ++i) {
			int cur = i;
			do 
				ord[i]++, cur = A[cur][i];
			while (cur != i);
		}
	}
	int T[30][30];
	int R[30][30];
	int cnt;
	int tmp[500010];
	inline bool Judge(int n) {
		for (int i = 0; i <= n; ++i) 
			for (int j = 0; j <= n; ++j) 
				if (A[i][j] <= n && mul(tmp[i], tmp[j]) != tmp[A[i][j]])
					return false;
		return true;
	}
	void Dfs(int n) {
		if (cnt == M) return;
		if (n == M) {
			bool zero = 1;
			for (int i = 0; i < M; ++i) if (tmp[i]) {
				zero = 0;
				break;
			}
			if (zero) return;
			for (int i = 0; i < M; ++i) T[cnt][i] = tmp[i];
			cnt++;
			return;
		}
		for (int i = 0; i <= ord[n]; ++i){
			tmp[n] = W[ord[n]][i];
			if (Judge(n)) Dfs(n+1);
		}
	}
	void Matrix_Inv() {
		static int C[30][30];
		memcpy(C, T, sizeof C);
		for (int i = 0; i < M; ++i) R[i][i] = 1;
		for (int i = 0; i < M; ++i) {
			int p = i;
			for (int j = i; j < M; ++j) if (C[j][i]) { p = j; break; }
			if (p != i) for (int j = 0; j < M; ++j) 
				std::swap(C[i][j], C[p][j]), std::swap(R[i][j], R[p][j]);
			for (int j = 0; j < M; ++j) if (i != j) {
				int rate = mul(C[j][i], qpow(C[i][i], MOD-2));
				for (int k = 0; k < M; ++k) {
					C[j][k] = add(C[j][k], MOD-mul(C[i][k], rate));
					R[j][k] = add(R[j][k], MOD-mul(R[i][k], rate));
				}
			}
		}
		for (int i = 0; i < M; ++i) {
			int inv = qpow(C[i][i], MOD-2);
			for (int j = 0; j < M; ++j) R[i][j] = mul(R[i][j], inv);
		}
	}
	void dft(int n, int *a, int C[30][30]) {
		if (n == 1) return;
		int b = n / M;
		for (int i = 0; i < M; ++i) dft(b, a + i * b, C);
		for (int i = 0; i < n; ++i) tmp[i] = 0;
		for (int i = 0; i < M; ++i)
			for (int j = 0; j < M; ++j)
				for (int k = 0; k < b; ++k)
					tmp[i * b + k] = add(tmp[i * b + k], mul(a[j * b + k], C[i][j]));
		for (int i = 0; i < n; ++i) a[i] = tmp[i];
	}
	void Main() {
		getord();
		Dfs(0);
		Matrix_Inv();
		dft(U, P, T);
		for (int i = 0; i < U; ++i) P[i] = qpow(P[i], (K+1)%(MOD-1));
		dft(U, P, R);
		for (int i = 0; i < U; ++i) printf("%d\n", P[i]);
	}

}
int main() {
#ifdef ARK
	freopen("F.6.0.in", "r", stdin);
	freopen("test.out", "w", stdout);
#endif
	for (int i = 1; i <= 22; ++i) {
		W[i][0] = 1;
		int rt = qpow(G, (MOD-1) / i);
		for (int j = 1; j < i; ++j) W[i][j] = mul(W[i][j-1], rt);
	}
	N = read(), M = read(), U = qpow(M, N), scanf("%lld", &K);
	for (int i = 0; i < M; ++i)
		for (int j = 0; j < M; ++j) A[i][j] = read();
	for (int i = 0; i < U; ++i) P[i] = read();
	Std_Dfs::Main();
}

总结

你过了这个题以后，抽象代数那点基础知识应该就能明白一点点了，就能看 LCA 的论文了。

而后你就能完全理解 $\text{DFT}$ 的本质和高位前缀和的思想，不仅是局限于背板子了。

同时你能达到理性颓废的巅峰，体会数学的美妙和偷税的乐趣。