構造 DFT :某少女附中的體育課 題解

構造 DFT :某少女附中的體育課 題解

題意

神仙題。一題能當十題寫。git

**一句話題意:**給出一個下標系統 $A$,求序列 $P$ 在該下標系統下的 $N$ 維卷積的 $K+1$ 次冪。算法

先來解釋一下什麼意思。spa

通常來講,卷積是一個以下的形式 $$ c_r = \sum_{p \circ q = r} a_pb_q $$ 當 $\circ$ 爲 $+$ 的時候,就是咱們熟悉的多項式乘法。code

那若是是位運算呢?就是 $\text{FWT}$ 乾的事情了。orm

那若是是其餘的東西呢?這個題就是給了你一個 $m \times m$ 的矩陣 $A$ ,表示一個運算表。ci

而後你有一個長度爲 $m^n$ 的序列 $P$ ,每一個下標是一個 $n$ 位 $m$ 進制數,下標之間的運算就是按位進行矩陣 $A$ 下的運算。你須要在這個意義下算出 $P$ 的 $K+1$ 次冪。get

爲了方便,如下設 $N = m^n$ 。數學

算法一

直接暴力計算卷積便可。string

預計得分: $7$。it

算法二

咱們注意到下標運算系統其實就是同或。直接上同或 $\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}$ 就好了。隨便搜搜就能過。跑的挺快的。

代碼超級好寫。

順便提一下,正解使用的是 Subtask7Subtask8的方法科學構造出那個變換矩陣。想學的能夠學一學。代碼不長。

預計得分:$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}$ 的本質和高位前綴和的思想,不僅是侷限於背板子了。

同時你能達到理性頹廢的巔峯,體會數學的美妙和偷稅的樂趣。

相關文章
相關標籤/搜索