NTT

單位根的定義就不說了。

顯然有：

 

\[\omega_n^k=\cos {2\pi k \over n}+i\sin {2\pi k \over n} \]

 

帶入這個能夠直接證得：

 

\[\omega_{2n}^{2k}=\cos {2\pi \cdot 2k \over 2n}+i\sin {2\pi \cdot 2k \over2 n}=\omega_n^k \]

 

咱們用圖像理性理解可得（就是繞着原點把向量轉了180度）：

 

\[\omega_n^{k+{n\over 2}}=-\omega_n^k \]

 

依然根據圖像可得（就是繞着原點旋轉了360度）:

 

\[\omega_n^{k+n}=\omega_n^k \]

 

由歐拉公式得：

 

\[\omega_n^k=\cos {2\pi k \over n}+i\sin {2\pi k \over n}=e^{i\cdot{2\pi k \over n}} \]

 

因此有：

 

\[\omega_n^k=e^{i\cdot{2\pi k \over n}}=(e^{i\cdot{2\pi \over n}})^k=(\omega_n^1)^k \]

 

咱們先帶入\(\omega_n^1,\omega_n^2,...,\omega_n^n\)到多項式\(A\)中，求出\(A(\omega_n^1),A(\omega_n^2),...,A(\omega_n^n)\)。

爲了方便，假設長度\(n\)爲\(2^k\)(高位不夠的添0便可)

把A下標奇偶分類：

 

\[A_1(x)=a_0+a_2 x+a_4 x^2+... \]

 

 

\[A_2(x)=a_1+a_3 x+a+5 x^2+... \]

 

顯然有：

 

\[A(x)=A_1(x^2)+xA_2(x^2) \]

 

 

\[\therefore A(\omega_{n}^{k})=A_1(\omega_{n}^{2k})+\omega_{n}^{k}A_2(\omega_{n}^{2k}) \]

 

 

\[=A_1(\omega_{n\over 2}^{k})+\omega_{n}^{k}A_2(\omega_{n\over 2}^{k}) \]

 

 

\[A(\omega_{n}^{k+{n\over 2}})=A_1(\omega_{n}^{2k+n})-\omega_{n}^{k}A_2(\omega_{n}^{2k+n}) \]

 

 

\[=A_1(\omega_{n\over 2}^{k})-\omega_{n}^{k}A_2(\omega_{n\over 2}^{k}) \]

 

這兩個式子只有後面一項是相反的，能夠遞歸求解。

因而給出代碼：

inline void FFT(complex<double> *a, int len) {
	if (!len) return ; complex<double> a1[len], a2[len];
	for (int i = 0; i < len; ++i) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
	FFT(a1, len >> 1); FFT(a2, len >> 1);
	complex<double> w(cos(PI / len), sin(PI / len)), wk(1, 0);
	for (int i = 0; i < len; ++i, wk *= w)
		a[i] = a1[i] + wk * a2[i], a[i + len] = a1[i] - wk * a2[i];
}

考慮怎麼從點值多項式轉換到係數多項式。

咱們欽定\(y_i=A(\omega_n^i)\)，在有一多項式\(C\)，知足:

 

\[C(x)=\sum y_i x^i \]

 

則咱們帶入\(\omega_n^{-k}\)，獲得：

 

\[C(\omega_n^{-k})=c_k=\sum_{i=0}^{n-1} y_i (\omega_n^{-k})^i \]

 

 

\[=\sum_{i=0}^{n-1} [\sum_{j=0}^{n-1} a_j(\omega_n^{i})^j ](\omega_n^{-k})^i \]

 

 

\[=\sum_{i=0}^{n-1} \sum_{j=0}^{n-1} a_j(\omega_n^{j})^i(\omega_n^{-k})^i \]

 

 

\[=\sum_{i=0}^{n-1} \sum_{j=0}^{n-1} a_j(\omega_n^{j-k})^i \]

 

 

\[=\sum_{i=0}^{n-1} a_i \sum_{j=0}^{n-1}(\omega_n^{i-k})^j \]

 

設:

 

\[S(\omega_n^k)=\sum_{i=0}^{n-1} (\omega_n^k)^i={(\omega_n^k)^{n}-1\over \omega_n^k-1} \]

 

當\(k\neq 0\)時爲0，不然爲\(n\)

則：

 

\[\sum_{j=0}^{n-1}(\omega_n^{i-k})^j=S(\omega_n^{i-k}) \]

 

即當\(i=k\)時爲\(n\)，因此：

 

\[c_k=\sum_{i=0}^{n-1} a_i \sum_{j=0}^{n-1}(\omega_n^{i-k})^j=na_k \]

 

咱們驚訝的發現這樣對\(C\)作一次FFT以後點值除以n就是多項式的係數了。

代碼結合一下：

inline void FFT(complex<double> *a, int len, int flag) {
	if (!len) return ; complex<double> a1[len], a2[len];
	for (int i = 0; i < len; ++i) a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
	FFT(a1, len >> 1, flag); FFT(a2, len >> 1, flag);
	complex<double> w(cos(PI / len), flag * sin(PI / len)), wk(1, 0);
	for (int i = 0; i < len; ++i, wk *= w)
		a[i] = a1[i] + wk * a2[i], a[i + len] = a1[i] - wk * a2[i];
}

發現遞歸版的碼會T，手玩一下發現實際上奇偶變換後下標的操做至關於二進制反過來，能夠改爲非遞歸來模擬，本身對着碼手玩看看就明白。

inline void FFT(complex *a, int type) {
	for (int i = 0; i < lim; ++i)
		if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; mid < lim; mid <<= 1) {
		complex wn; wn = complex(cos(pi / mid), type * sin(pi / mid));
		for (int j = 0; j < lim; j += mid << 1) {
			complex bas; bas = complex(1, 0);
			for (int k = 0; k < mid; ++k, bas = bas * wn) {
				complex x = a[j + k], y = bas * a[j + mid + k];
				a[j + k] = x + y;
				a[j + mid + k] = x - y;
			}
		}
	}
}

解釋一下，第一層循環枚舉的是遞歸的層數，即當前合併的兩個多項式的長度。第二層就是枚舉當前要合併多項式的起點，第三層就是枚舉的具體的那一個係數。這麼說不是很清楚，仍是本身造樣例跟着代碼手玩一下就明白了。

而NTT呢？設模數爲p，g是p的原根，則不須要證實的給出，\(\omega_n^1\)等價於\(g^{p-1\over n} \bmod p\)。把上面的碼代碼裏的wn換成這個就好了。通常p=998244353，此時g=3。

給個板子：

struct poly {
	int n;
	vector<ll> x;
	
	inline void NTT(int flag) {
		for (int i = 0; i < n; ++i)
			if (i < rev[i]) swap(x[i], x[rev[i]]);
		for (int mid = 1; mid < n; mid <<= 1) {
			ll wn = power(flag == 1 ? G : Gi, (mod - 1) / (mid << 1));
			for (int j = 0; j < n; j += mid << 1) {
				ll bas = 1;
				for (int k = 0; k < mid; ++k, bas = (bas * wn) % mod) {
					ll xx = x[j + k], y = (bas * x[j + mid + k]) % mod;
					x[j + k] = (xx + y) % mod;
					x[j + mid + k] = ((xx - y) % mod + mod) % mod;
				}
			} cerr << endl;
		}
	}
	
};

inline int max_(int a, int b) {
	return a > b ? a : b;
}

inline poly mul(poly A, poly B) {
	poly a, b; a = A; b = B;
	int tmp = a.n + b.n; a.n = 1;
	int L = 0;
	while (a.n <= tmp) a.n <<= 1, ++L;
	for (int i = 0; i <= a.n; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L - 1);
	b.n = a.n;
	a.NTT(1); b.NTT(1);
	for (int i = 0; i < a.n; ++i) a.x[i] = (a.x[i] * b.x[i]) % mod;
	a.NTT(-1);
	const ll inv = power(a.n, mod - 2);
	a.n = tmp;
	for (int i = 0; i <= a.n; ++i) a.x[i] = (a.x[i] * inv) % mod;
	return a;
}