多項(xiàng)式基礎(chǔ) - 學(xué)習(xí)筆記

1. FFT

前置知識(shí)：復(fù)數(shù)

基礎(chǔ)定義

令 \(i\) 為方程 \(x^2 = -1\) 的解

我們定義，形如 \(z = \bm{a+bi}\)，其中 \(a, b \in \mathbb{R}\) 的數(shù) \(z\)，稱為復(fù)數(shù)；對(duì)于 \(z = a+bi\)，我們稱 \(a\) 為其實(shí)部，\(bi\) 為其虛部

定義復(fù)數(shù) \(z_1 = z_2\)，當(dāng)且僅當(dāng) \(\begin{cases} a_1 = a_2 \\ b_1 = b_2 \end{cases}\)；由此：

• 復(fù)數(shù) \(z = a+bi\) 唯一對(duì)應(yīng)平面直角坐標(biāo)系中一點(diǎn) \((a, b)\)；換句話說(shuō)，復(fù)數(shù)集與平面直角坐標(biāo)系中的點(diǎn)一一對(duì)應(yīng)
• 復(fù)數(shù) \(z = a+bi\) 也唯一對(duì)應(yīng)向量 \(\overrightarrow{OZ}\)，其中點(diǎn) \(Z\) 的坐標(biāo)為 \((a, b)\)

我們稱對(duì)應(yīng)復(fù)數(shù)集的平面直角坐標(biāo)系為復(fù)平面，x 軸為實(shí)軸，y 軸為虛軸

類似向量，定義復(fù)數(shù) \(z = a+bi\) 的模為 \(\bm{|z| = \sqrt{a^2+b^2}}\)

四則運(yùn)算

令 \(z_1 = a_1+b_1i\)，\(z_2 = a_2+b_2i\)，定義：

• 復(fù)數(shù)加法 \(z_1+z_2 = (a_1+a_2)+(b_1+b_2)i\)
• 復(fù)數(shù)減法 \(z_1-z_2 = (a_1-a_2)+(b_1-b_2)i\)
• 復(fù)數(shù)乘法 \(z_1z_2 = (a_1+b_1i)(a_2+b_2i) = a_1a_2+b_1b_2i^2+(a_1b_2+b_1a_2)i\)
  帶入 \(i^2 = -1\)，可得 \(z_1z_2 = (a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i\)
• 復(fù)數(shù)除法 \(\displaystyle \frac{z_1}{z_2} = \frac{a_1+b_1i}{a_2+b_2i}\)
  為使分母化為實(shí)數(shù)，上下同乘 \(a_2-b_2i\)，可得 \(\displaystyle \frac{z_1}{z_2} = \frac{(a_1+b_1i)(a_2-b_2i)}{(a_2+b_2i)(a_2-b_2i)} = \frac{a_1a_2+b_1b_2}{a_2^2+b_2^2} + \frac{a_2b_1-a_1b_2}{a_2^2+b_2^2}i\)

手寫(xiě)結(jié)構(gòu)體代碼如下：

#define db double
struct cp{
	db x, y; cp(db xx = 0, db yy = 0){ x = xx; y = yy; }
	cp operator + (const cp& b){ return cp(x+b.x, y+b.y); }
	cp operator - (const cp& b){ return cp(x-b.x, y-b.y); }
	cp operator * (const cp& b){ return cp(x*b.x-y*b.y, x*b.y+y*b.x); }
	cp operator / (const cp& b){ db t = (b.x*b.x+b.y*b.y); return cp((x*b.x+y*b.y)/t, (y*b.x-x*b.y)/t); }
}

輻角與輻角主值

非零復(fù)數(shù) \(z = a+bi\) 對(duì)應(yīng)點(diǎn) \((a, b)\)，我們用極坐標(biāo)將其表示為 \((r, \theta)\)，其中 \(\bm{r = |z|}\)

記向量與實(shí)軸的夾角 \(\theta\) 為 \(z\) 的輻角，記為 \(\theta = \text{Arg}\ z\)，顯然有 \(\displaystyle \tan \theta = \frac{b}{a}\)

注意到任一非零復(fù)數(shù) \(z\) 都有無(wú)窮多個(gè)輻角，我們用 \(\text{arg}\ z\) 表示其中滿足 \(-\pi < \theta < \pi\) 的特定 \(\theta\)，稱為輻角主值

由此，可以引出復(fù)數(shù)相乘的重要規(guī)律：復(fù)數(shù)相乘，模長(zhǎng)相乘，輻角相加 (非常重要)

證明：設(shè) \(z_1 = a_1+b_1i\)，\(z_2 = a_2+b_2i\)：

• 先證模長(zhǎng)相乘，\(\text{LHS} = \sqrt{a_1^2+b_1^2} \times \sqrt{a_2^2+b_2^2}\)，\(\text{RHS} = \sqrt{(a_1a_2-b_1b_2)^2+(a_1b_2+a_2b_1)^2}\)，展開(kāi)，顯然有 \(\text{LHS} = \text{RHS}\)
• 再證輻角相加，設(shè) \(z_1 = \overrightarrow{OA}\)，\(z_2 = \overrightarrow{OB}\)，相乘結(jié)果為 \(\overrightarrow{OC}\)；取 \(P(1, 0)\)，考慮證 \(\triangle OAP\) 與 \(\triangle OCB\) 相似
  • 由模長(zhǎng)相乘，有 \(\overline{OA} \times \overline{OB} = \overline{OC}\)；又顯然有 \(\overline{OP} \times \overline{OB} = \overline{OB}\)，那么只需證 \(\overline{AP} \times \overline{OB} = \overline{CB}\)
  • 套下兩點(diǎn)距離公式暴力展開(kāi)即可，不難證得相似
  • 那么有 \(\angle COP = \angle BOP + \angle BOC = \angle BOP + \angle AOP\)，得證

單位根

我們定義滿足 \(\bm{|z| = 1}\) 的復(fù)數(shù) \(z\) 構(gòu)成單位圓

稱方程 \(x^n = 1\) 在復(fù)數(shù)意義下的解為 \(n\) 次單位根；這樣的解必恰有 \(n\) 個(gè)，且 \(n\) 次單位根將單位圓 \(n\) 等分

證明：

• 由模長(zhǎng)相乘，\(|z| < 1\) 時(shí)會(huì)越乘越小，\(|z| > 1\) 時(shí)會(huì)越乘越大；于是單位根必然在單位圓上
• 由輻角相加，有 \(n \theta = 2k \pi\ (k \in \mathbb{Z})\)
  • 取 \(\displaystyle \theta = \frac{1(2 \pi)}{n}, \frac{2(2 \pi)}{n}, \frac{3(2 \pi)}{n}, \cdots, \frac{(n-1)(2 \pi)}{n}\)，至少有這 \(n\) 個(gè)解
  • 由代數(shù)基本定理，有且僅有這 \(n\) 個(gè)解；顯然它們將單位圓 \(n\) 等分，證畢

我們記這些單位根為 \(\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}\) (以下均省略指數(shù)上的 \(\bmod n\))，它們有如下性質(zhì)：

• \(\omega_{n}^{0} = \omega_{n}^{n} = 1\)
• \(\omega_{n}^{k} = (\omega_{n}^{1})^k\)，由模長(zhǎng)相乘輻角相加，顯然
• \(\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i+j}\)，先轉(zhuǎn) \(i\) 次再轉(zhuǎn) \(j\) 次即為轉(zhuǎn) \(i+j\) 次
• \(\bm{\omega_{n}^{k} = \omega_{2n}^{2k}}\)，即將每段細(xì)分為兩份 (重要)
• \(\bm{\omega_{2n}^{k+n} = -\omega_{2n}^{k}}\)，由 \(\omega_{2n}^{k+n} = \omega_{2n}^{k} \times \omega_{2n}^{n}\)，多乘 \(\omega_{2n}^{n}\) 相當(dāng)于多轉(zhuǎn)半圈，關(guān)于原點(diǎn)對(duì)稱 (重要)

DFT 與 IDFT 思想

P.S. 下面使用 \(F_{i}\) 表示多項(xiàng)式 \(F(x)\) 的 \(i\) 次項(xiàng)系數(shù)

兩多項(xiàng)式相乘(卷積)定義為 \(H(x) = \displaystyle \sum_{i=0} \sum_{j=0} F_iG_jx^{i+j}\)

• 令 \(n\) 次多項(xiàng)式 \(F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\)；這稱為多項(xiàng)式的系數(shù)表達(dá)

  系數(shù)表達(dá)較為常用；但系數(shù)表達(dá)直接算卷積是 \(O(n^2)\) 的

• 由 \(n\) 個(gè)等式可以解出 \(n\) 元方程的解，我們可以用 \(n+1\) 個(gè)點(diǎn)唯一確定 \(F(x)\)；這稱為多項(xiàng)式的點(diǎn)值表達(dá)

  這種表達(dá)方法的好處是：

  • 若我們知道 \(F(x)\) 與 \(G(x)\) 在 \(x_0\) 處的點(diǎn)值 \(F(x_0)\) 與 \(G(x_0)\)，直接就有 \(H(x_0) = F(x_0) \times G(x_0)\)

    換句話說(shuō)，若我們知道 \(F(x)\) 與 \(G(x)\) 的點(diǎn)值表達(dá)，可以 \(O(n)\) 算出 \(H(x)\) 的點(diǎn)值表達(dá)

  • 系數(shù)表達(dá)允許我們帶入具有特殊性質(zhì)的值，這有助于創(chuàng)造好的性質(zhì)

由上，考慮使用點(diǎn)值表達(dá)解決問(wèn)題

那么，現(xiàn)在只缺將系數(shù)表達(dá)轉(zhuǎn)化為點(diǎn)值表達(dá)(求值)，以及將點(diǎn)值表達(dá)轉(zhuǎn)化為系數(shù)表達(dá)(插值)的算法

樸素系數(shù)轉(zhuǎn)點(diǎn)值算法即為 DFT，樸素點(diǎn)值轉(zhuǎn)系數(shù)算法即為 IDFT

FFT 原理

將樸素 DFT 算法優(yōu)化，可得 FFT

FFT 的主要思想是利用分治來(lái)優(yōu)化時(shí)間復(fù)雜度

設(shè)有 \(n = 2^k\ (k \in \mathbb{N^*})\)，對(duì)于項(xiàng)數(shù)不足 \(2^k\) (次數(shù)不足 \(2^k-1\) ) 的多項(xiàng)式 ，我們?cè)谄?strong>高次項(xiàng)補(bǔ)零即可

將多項(xiàng)式 \(F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\) 根據(jù)次數(shù)奇偶拆成兩半，分別為：

• \(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2}\)
• \(a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1}\)

我們根據(jù)這兩部分構(gòu)造兩個(gè)多項(xiàng)式 \(FL(x)\) 與 \(FR(x)\)，具體的：

• \(FL(x) = a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}\)
• \(FR(x) = a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1}\)

那么顯然有關(guān)系式：\(\displaystyle \bm{F(x) = FL(x^2)+xFR(x^2)}\)

考慮帶入 \(\omega_{n}^{k}\)，發(fā)掘性質(zhì)：

\(\displaystyle F(\omega_{n}^{k}) = FL((\omega_{n}^{k})^2)+\omega_{n}^{k}FR((\omega_{n}^{k})^2)\)

\(\displaystyle = FL(\omega_{n}^{2k})+\omega_{n}^{k}FR(\omega_{n}^{2k})\)

\(\displaystyle = FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})\)

由此，只要知道 \(FL(x)\) 與 \(FR(x)\) 在 \(\displaystyle \omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) 處的取值，就可以得到 \(F(x)\) 在 \(\displaystyle \omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{\frac{n}{2}-1}\) 處的取值

對(duì)于后半段，考慮帶入 \(\displaystyle \omega_{n}^{k+\frac{n}{2}}\)：

\(\displaystyle F(\omega_{n}^{k+\frac{n}{2}}) = FL((\omega_{n}^{k+\frac{n}{2}})^2)+\omega_{n}^{k+\frac{n}{2}}FR((\omega_{n}^{k+\frac{n}{2}})^2)\)

\(\displaystyle = FL(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}FR(\omega_{n}^{2k+n})\)

\(\displaystyle = FL(\omega_{n}^{2k})+\omega_{n}^{k+\frac{n}{2}}FR(\omega_{n}^{2k})\)

\(\displaystyle = FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})\)

由此，只要知道 \(FL(x)\) 與 \(FR(x)\) 在 \(\displaystyle \omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) 處的取值，就可以得到 \(F(x)\) 在 \(\displaystyle \omega_{n}^{\frac{n}{2}}, \omega_{n}^{\frac{n}{2}+1}, \cdots, \omega_{n}^{n-1}\) 處的取值

拼起來(lái)，我們已經(jīng)可以得到 \(F(x)\) 在 \(\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}\) 處的所有值了

總結(jié)下兩個(gè)關(guān)鍵式子：

• \(\bm{\displaystyle F(\omega_{n}^{k}) = FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})}\)
• \(\bm{\displaystyle F(\omega_{n}^{k+\frac{n}{2}})= FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})}\)

據(jù)此遞歸求解即可，時(shí)間復(fù)雜度 \(O(n \log n)\)

FFT 實(shí)現(xiàn)細(xì)節(jié)

求單位根

小技巧：使用 acos(-1) 可以取到精確的 \(\pi\) 值

由三角函數(shù)相關(guān)知識(shí)，\(\displaystyle \bm{\omega_{n}^{1} = (\cos(\frac{2\pi}{n}), \sin(\frac{2\pi}{n}))}\)，從 \((1, 0)\) 開(kāi)始一直乘即可

IDFT 的實(shí)現(xiàn)

結(jié)論：將 DFT 過(guò)程中的 \(\omega_{n}^{i}\) 換為 \(\omega_{n}^{-i}\)，做完后除以 \(n\) 即得到系數(shù)表達(dá)

證明：設(shè)點(diǎn)值序列為 \(G(x)\)，則有 \(\displaystyle G(k) = \sum_{i=0}^{n-1} F_i(\omega_{n}^{k})^i\)

現(xiàn)在要證 \(\displaystyle n \times F_k = \sum_{i=0}^{n-1} G_i(\omega_{n}^{-k})^i\)；考慮暴力帶入：
\(\displaystyle \sum_{i=0}^{n-1} G_i(\omega_{n}^{-k})^i = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} F_j(\omega_{n}^{i})^j (\omega_{n}^{-k})^i = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} F_j(\omega_{n}^{ij-ik}) = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} F_j(\omega_{n}^{i(j-k)})\)

考慮對(duì) \(j\) 與 \(k\) 的大小分類討論：

• \(j = k\)，原式為 \(\displaystyle \sum_{i=0}^{n-1} F_k(\omega_{n}^{0}) = n \times F_k\)

• \(j \ne k\)，對(duì)每個(gè) \(j\)，不妨設(shè) \(t = j-k\)，原式為 \(\displaystyle \sum_{i=0}^{n-1} F_{t+k}(\omega_{n}^{it})\)

  由等比數(shù)列求和，\(\displaystyle \sum_{i=0}^{n-1} \omega_{n}^{it} = \frac{\omega_{n}^{nt}-\omega_{n}^{0}}{\omega_{n}^{t}-1} = \frac{\omega_{n}^{0}-\omega_{n}^{0}}{\omega_{n}^{t}-1} = 0\)

  那么每個(gè) \(j \ne k\) 的 \(j\) 都沒(méi)有貢獻(xiàn)

綜上，總貢獻(xiàn)為 \(n \times F_k\)，證畢

void fft(cp *f, int l, bool flg){ // 0->IDFT, 1->DFT
	if (l == 1) return;
	cp *fl = f, *fr = f+l/2;
	for (int i=0; i<l; i++) tmp[i] = f[i];
	for (int i=0; i<l/2; i++) fl[i] = tmp[(i<<1)], fr[i] = tmp[(i<<1|1)];
	fft(fl, l/2, flg); fft(fr, l/2, flg); cp w1(cos(2.0*Pi/l), sin(2.0*Pi/l)), w(1, 0);
	if (flg == 0) w1.y = -w1.y;
	for (int i=0; i<l/2; i++){ tmp[i] = fl[i]+w*fr[i]; tmp[i+l/2] = fl[i]-w*fr[i]; w = w*w1; }
	for (int i=0; i<l; i++) f[i] = tmp[i];
}

蝶形運(yùn)算優(yōu)化

每次分治，我們將多項(xiàng)式 \(F(x)\) 按照次數(shù)奇偶分成左右兩半，使用了大量?jī)?nèi)存移動(dòng)

若我們一開(kāi)始就知道 \(F_i\) 最終被分到哪個(gè)地方，就可以省去分開(kāi)的過(guò)程，直接倍增合并，達(dá)到優(yōu)化常數(shù)的目的

不妨設(shè) \(p(n, k)\) 表示項(xiàng)數(shù)為 \(2^n\) 時(shí)，\(k\) 次項(xiàng)最終被分到哪里：

• 若 \(n = 0\)，顯然 \(p(n, k) = k\)，原地不動(dòng)
• 若 \(n > 0\)，有 \(p(n, k) = ((k\ \texttt{&}\ 1)\ \texttt{<<}\ (n-1))+p(n-1, k\ \texttt{>>}\ 1)\)，奇數(shù)項(xiàng)放右邊，偶數(shù)項(xiàng)放左邊

展開(kāi)可得：\(\displaystyle p(n, k) = \sum_{i=0}^{n-1} (((k\ \texttt{>>}\ i)\ \texttt{&}\ 1)\ \texttt{<<}\ (n-i+1))\)

這顯然是將 \(k\) 二進(jìn)制反轉(zhuǎn)

于是我們預(yù)處理數(shù)組 \(tr_i\) 表示 \(i\) 的二進(jìn)制反轉(zhuǎn)，在 FFT 前對(duì)每組 \((i, tr_i)\) 對(duì)應(yīng)的系數(shù)交換一次即可

下面考慮如何預(yù)處理 \(tr_i\)；我們將數(shù) \(i\) 看作最低位與其他位兩部分組成

• \((tr_{(i\ \texttt{>>}\ 1)})\ \texttt{>>}\ 1\) 可求出其他位的二進(jìn)制反轉(zhuǎn)；多右移一次是因?yàn)榉崔D(zhuǎn)前 \(i\ \texttt{>>}\ 1\) 時(shí)在最高位多帶進(jìn)來(lái)了個(gè) \(0\)
• 若 \(i\ \texttt{&}\ 1\) 為 \(1\)，反轉(zhuǎn)后在最高位補(bǔ)上 \(1\) 即可

預(yù)處理代碼為 for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0));

合并過(guò)程優(yōu)化

蝶形運(yùn)算優(yōu)化后，注意到項(xiàng)數(shù)為 \(2^n\) 時(shí)，對(duì)于 \(k\) 次項(xiàng)：

• \(\displaystyle FL(\omega_{\frac{n}{2}}^{k})\) 存儲(chǔ)在數(shù)組下標(biāo)為 \(k\) 的位置，\(\displaystyle F(\omega_{n}^{k})\) 應(yīng)該存儲(chǔ)在數(shù)組下標(biāo)為 \(k\) 的位置
• \(\displaystyle FR(\omega_{\frac{n}{2}}^{k})\) 存儲(chǔ)在數(shù)組下標(biāo)為 \(k+\frac{n}{2}\) 的位置，\(\displaystyle F(\omega_{n}^{k})\) 應(yīng)該存儲(chǔ)在數(shù)組下標(biāo)為 \(k+\frac{n}{2}\) 的位置

因此新開(kāi)臨時(shí)數(shù)組去記錄是非常浪費(fèi)的，只需要一個(gè)臨時(shí)變量即可

關(guān)鍵代碼 (可以結(jié)合后面的完整代碼理解)：

cp tmp = w*f[k+nl];
f[k+nl] = f[k]-tmp;
f[k] = f[k]+tmp; w = w*w1;

三次變兩次

注意到復(fù)數(shù)有 \(x\) 與 \(y\) 兩維，可以同時(shí)存儲(chǔ)兩份信息；我們考慮利用這一性質(zhì)

記復(fù)多項(xiàng)式 \(H(x) = F(x)+G(x)i\)，則有：

\(\displaystyle H^2(x) = (F(x)+G(x)i)^2 = F^2(x)-G^2(x)+2F(x)G(x)i\)

那么只需要一次 DFT，一次 IDFT 即可，最終記得多除 \(2\)

最終代碼

P3803 【模板】多項(xiàng)式乘法（FFT）

#include<bits/stdc++.h>
using namespace std;
#define db double
const db Pi = acos(-1);
const int MAXN = 4e6+5;

int n, m, l;
int tr[MAXN];

struct cp{
	db x, y; cp(db xx = 0, db yy = 0){ x = xx; y = yy; }
	cp operator + (const cp& b){ return cp(x+b.x, y+b.y); }
	cp operator - (const cp& b){ return cp(x-b.x, y-b.y); }
	cp operator * (const cp& b){ return cp(x*b.x-y*b.y, x*b.y+y*b.x); }
	cp operator / (const cp& b){ db t = (b.x*b.x+b.y*b.y); return cp((x*b.x+y*b.y)/t, (y*b.x-x*b.y)/t); }
}f[MAXN];

void fft(cp *f, int l, bool flg){ // 0->IDFT, 1->DFT
	for (int i=0; i<l; i++) if (i < tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){
		int nl = (i>>1);
		cp w1(cos(2*Pi/i), sin(2*Pi/i));
		if (flg == 0) w1.y = -w1.y;
		for (int j=0; j<l; j+=i){
			cp w(1, 0);
			for (int k=j; k<j+nl; k++){
				cp tmp = w*f[k+nl];
				f[k+nl] = f[k]-tmp;
				f[k] = f[k]+tmp; w = w*w1;
			}
		}
	}
}

int main(){
	scanf("%d%d", &n, &m);
	for (int i=0; i<=n; i++) scanf("%lf", &f[i].x);
	for (int i=0; i<=m; i++) scanf("%lf", &f[i].y);
	for (l=1; l<=n+m; l<<=1); 
	for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0));
	fft(f, l, 1); for (int i=0; i<l; i++) f[i] = f[i]*f[i]; fft(f, l, 0);
	for (int i=0; i<=n+m; i++) printf("%d ", (int)(f[i].y/l/2+0.49));
	return 0;
}

2. NTT

NTT 原理

FFT 需要依賴復(fù)數(shù)單位根，容易引發(fā)精度誤差

大多數(shù)問(wèn)題都在模意義下進(jìn)行，因此，我們希望為 FFT 找一個(gè)模意義下的替代；于是考慮引入原根 \(g\) 代替單位根

不過(guò)，對(duì)于模數(shù) \(p\) 的原根 \(g\)，其階數(shù)為 \(\varphi(p)=p-1\)；現(xiàn)在對(duì)于 \(n = 2^k\ (k \in \mathbb{N^*})\)，我們希望找到階數(shù)為 \(n\) 的數(shù) \(x\)，以此保證 \(x^0, x^1, \cdots, x^{n-1} (\bmod p)\) 互不相同；顯然直接用 \(g\) 是不行的

注意到 \(g^k\) 的階數(shù)為 \(\displaystyle \frac{p-1}{\gcd(k, p-1)}\)，取 \(\displaystyle k = \frac{p-1}{n}\) 即可

對(duì) \(g^k\) 階數(shù)的證明：可將 \(g^1, g^2, \cdots, g^{p-1}\) 看成環(huán)，每次乘 \(g^k\) 相當(dāng)于環(huán)上走 \(k\) 步

• \(k\ \nmid\ (p-1)\)，顯然可以走到環(huán)上所有位置，階數(shù)為 \(\displaystyle p-1 = \frac{p-1}{\gcd(k, p-1)}\)
• \(k\ |\ (p-1)\)，不妨設(shè) \(\displaystyle k = \frac{p-1}{i}\)，則共能走到 \(0, k, 2k, 3k, \cdots, (i-1)k\ (\bmod p)\) 共 \(i\) 個(gè)位置，階數(shù)為 \(\displaystyle i = \frac{p-1}{k} = \frac{p-1}{\gcd(k, p-1)}\)

記 \(\displaystyle t = g^{\frac{p-1}{n}}\) (或 \(t_n\))，下面說(shuō)明 \(t\) 具有用到的所有單位根的性質(zhì)：

• \(\omega_{n}^{k} = \omega_{n}^{k \bmod n}\)：在小環(huán)上轉(zhuǎn)圈與在單位圓上轉(zhuǎn)圈類似，多轉(zhuǎn) \(n\) 次就轉(zhuǎn)回來(lái)了，易得 \(t^{k} = t^{k \bmod n}\)；以下省略指數(shù)上取模
• \(\omega_{n}^{k} = (\omega_{n}^{1})^k\)，更進(jìn)一步，\(\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i+j}\)：顯然 \(t^i \times t^j = t^{i+j}\)
• \(\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}\) 互不相同：小環(huán)大小為 \(n\)，顯然 \(t^{0}, t^{1}, \cdots, t^{n-1}\) 互不相同
• \(\displaystyle \omega_{n}^{k+\frac{n}{2}} = -\omega_{n}^{k}\)：
  • 即證 \(\displaystyle t^{k+\frac{n}{2}} \equiv -t^{k}\ (\bmod p)\)，移項(xiàng)，即 \(\displaystyle t^{k}(t^{\frac{n}{2}}+1) \equiv 0\ (\bmod p)\)；下面證 \(\displaystyle t^{\frac{n}{2}} \equiv -1\ (\bmod p)\)
  • 顯然有 \(t^n \equiv 1\ (\bmod p)\) (帶入即證)，即 \(\displaystyle (t^{\frac{n}{2}})^2 \equiv 1\ (\bmod p)\)；那么 \(\displaystyle t^{\frac{n}{2}} \bmod p\) 的取值只能為 \(\pm 1\)
  • 若 \(\displaystyle t^{\frac{n}{2}} \equiv 1\ (\bmod p)\)，\(t\) 的階數(shù)就不是 \(n\)，而是 \(\displaystyle \frac{n}{2}\) 了，顯然矛盾；那么 \(\displaystyle t^{\frac{n}{2}} \equiv -1\ (\bmod p)\)，證畢
• \(\omega_{2n}^{2k} = \omega_{n}^{k}\)：直接帶入，\(\displaystyle t_{2n} = (g^{\frac{p-1}{2n}})^2 = g^{\frac{p-1}{n}} = t_{n}\)

NTT 實(shí)現(xiàn)細(xì)節(jié)

我們直接將 FFT 中的單位根 \(\omega_{n}^{k}\) 換成 \(t^k\) 即可，求 \(t^1\) 可以直接快速冪，IDFT 中遇到的 \(t^{-1}\) 即為逆元

NTT 要求 \(n\ |\ (p-1)\)，換句話說(shuō)，\(p-1\) 中需要含有大量 \(2\) 的冪；同時(shí)可能需要記住常見(jiàn)模數(shù)的原根

P3803 NTT 代碼：

#include<bits/stdc++.h>
using namespace std;
#define db double
#define ll long long
const db Pi = acos(-1);
const int MAXN = 4e6+5;
const ll MOD = 998244353, G = 3, INVG = 332748118;

int n, m, l;
int tr[MAXN], f[MAXN], g[MAXN];

ll QuickPow(ll a, ll b){
	ll res = 1;
	while (b){
		if (b&1) res = (res*a)%MOD;
		a = (a*a)%MOD;
		b >>= 1;
	} return res;
}

void ntt(int *f, int l, bool flg){ // 0->IDFT, 1->DFT
	for (int i=0; i<l; i++) if (i < tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){
		int nl = (i>>1), w1 = QuickPow((flg == 1 ? G : INVG), (MOD-1)/i);
		for (int j=0; j<l; j+=i){
			ll w = 1;
			for (int k=j; k<j+nl; k++){
				int tmp = (w*f[k+nl])%MOD;
				f[k+nl] = f[k]-tmp; if (f[k+nl] < 0) f[k+nl] += MOD;
				f[k] = f[k]+tmp; if (f[k] >= MOD) f[k] -= MOD; 
				w = (w*w1)%MOD;
			}
		}
	}
}

int main(){
	scanf("%d%d", &n, &m);
	for (int i=0; i<=n; i++) scanf("%d", &f[i]);
	for (int i=0; i<=m; i++) scanf("%d", &g[i]);
	for (l=1; l<=n+m; l<<=1); int inv = QuickPow(l, MOD-2);
	for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0));
	ntt(f, l, 1); ntt(g, l, 1); for (int i=0; i<l; i++) f[i] = (1ll*f[i]*g[i])%MOD; ntt(f, l, 0);
	for (int i=0; i<=n+m; i++) printf("%d ", (1ll*f[i]*inv)%MOD);
	return 0;
}

這里也放個(gè)封裝版本：

void pre(int l){ for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0)); }
void ntt(int *g, int l, bool flg){ // 0->IDFT, 1->DFT, l>n
	pre(l); static int f[MAXN] = {0};
	for (int i=0; i<l; i++) f[i] = g[i]; for (int i=0; i<l; i++) if (i<tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){ int nl = (i>>1), w1 = QuickPow((flg?G:INVG), (MOD-1)/i);
		for (int j=0; j<l; j+=i){ ll w = 1;
			for (int k=j; k<j+nl; k++){
				int tmp = (w*f[k+nl])%MOD;
				f[k+nl] = f[k]-tmp; if (f[k+nl]<0) f[k+nl]+=MOD;
				f[k] = f[k]+tmp; if (f[k]>=MOD) f[k]-=MOD; w = (w*w1)%MOD;}}}
	if (flg == 0){ int inv = QuickPow(l, MOD-2); for (int i=0; i<l; i++) g[i] = (1ll*f[i]*inv)%MOD; }
	else{ for (int i=0; i<l; i++) g[i] = f[i]; }
	for (int i=0; i<(l<<1); i++) f[i] = 0;           
}
void mulp(int *f, int *g, int l){ for (int i=0; i<l; i++) f[i] = (1ll*f[i]*g[i])%MOD; } // l > n
void times(int *f, int *g, int n, int m){ // max{len} = n, mod x^m
	static int tmp[MAXN] = {0};
	int l = 1; for (l=1; l<=n+n; l<<=1); for (int i=0; i<l; i++) tmp[i] = g[i];
	ntt(f, l, 1); ntt(tmp, l, 1); mulp(f, tmp, l); ntt(f, l, 0);
	for (int i=0; i<l; i++) tmp[i] = 0; for (int i=m; i<l; i++) f[i] = 0;
}

函數(shù)多了很容易搞混，包括后面，實(shí)現(xiàn)時(shí)千萬(wàn)注意傳的參是次數(shù)還是 \(n = 2^k\ (k \in \mathbb{N^*})\)！

3. 多項(xiàng)式全家桶

多項(xiàng)式四則運(yùn)算

多項(xiàng)式加減法

多項(xiàng)式的加減定義為：

• \(\displaystyle F(x)+G(x) = \sum_{i=0} (F_i+G_i)x^i\)
• \(\displaystyle F(x)-G(x) = \sum_{i=0} (F_i-G_i)x^i\)

系數(shù)表達(dá)直接相加 / 相減即可，可以 \(O(n)\) 算出

多項(xiàng)式乘法

使用上文所述 FFT 或 NTT 即可，可以 \(O(n \log n)\) 算出

多項(xiàng)式乘法逆元

當(dāng) \(F(x)G(x) = 1\) 時(shí)，我們稱 \(F(x)\) 與 \(G(x)\) 互為乘法逆元；其可幫助我們將(整)除轉(zhuǎn)化為乘

考慮次數(shù)倍增，不斷擴(kuò)大次數(shù)上界，直到求出完整的逆元

不妨設(shè) \(R_{\triangle}(x)\) 為 \(F(x)\) 的前 \(\displaystyle \frac{n}{2}\) 位逆元，考慮如何求出 \(R(x)\) 表示 \(F(x)\) 的前 \(n\) 位逆元

顯然有：\(\displaystyle R(x)-R_{\triangle}(x) \equiv 0\ (\bmod x^{\frac{n}{2}})\)；考慮平方求得 \(\bmod x^n\) 下的逆元：

\(\displaystyle (R(x)-R_{\triangle}(x))^2 \equiv 0\ (\bmod x^{n})\)

\(\displaystyle R^2(x)-2R(x)R_{\triangle}(x)+R_{\triangle}^2(x) \equiv 0\ (\bmod x^{n})\)，乘個(gè) \(F(x)\) 消去 \(R(x)\)：

\(\displaystyle F(x)(R^2(x)-2R(x)R_{\triangle}(x)+R_{\triangle}^2(x)) \equiv 0\ (\bmod x^{n})\)

\(\displaystyle R(x)-2R_{\triangle}(x)+R_{\triangle}^2(x)F(x) \equiv 0\ (\bmod x^{n})\)

\(\displaystyle R(x) \equiv 2R_{\triangle}(x)-R_{\triangle}^2(x)F(x)\ (\bmod x^{n})\)

根據(jù)上式倍增即可；邊界條件為 \(\bmod x^1\) 時(shí)就是求數(shù)字逆元

void invp(int *f, int n){ // len = n, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}, tmp[MAXN] = {0}; w[0] = QuickPow(f[0], MOD-2);
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; }
		for (int j=0; j<i; j++) tmp[j] = f[j];
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(tmp, (i<<1), 1); mulp(w, tmp, (i<<1)); ntt(w, (i<<1), 0); 
        // 三個(gè)多項(xiàng)式相乘，次數(shù)超過(guò) i，開(kāi)到 (i<<1) 避免循環(huán)卷積溢出破壞前面幾位
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = r[j]-w[j]; if (w[j] < 0) w[j] += MOD; }
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = tmp[i] = 0;
}

多項(xiàng)式帶余除法

給定 \(n\) 次多項(xiàng)式 \(F(x)\) 與 \(m\) 次多項(xiàng)式 \(G(x)\)，令 \(F(x) = Q(x)G(x)+R(x)\)，求 \(n-m\) 次商多項(xiàng)式 \(Q(x)\) 與 \(m-1\) 次余多項(xiàng)式 \(R(x)\)

我們希望利用 \(\bmod x^k\) 的方式消去 \(R(x)\)；但這種方法只適用于消高次項(xiàng)，不適用于消低次項(xiàng)，于是我們考慮將多項(xiàng)式反轉(zhuǎn)

定義 \(n\) 次多項(xiàng)式 \(F(x)\) 的反轉(zhuǎn)為 \(\displaystyle F^{T}(x) = x^nF(\frac{1}{x})\)，容易發(fā)現(xiàn)即為將其系數(shù)前后反轉(zhuǎn)

例：\(F(x) = 3+x+4x^2\)，則有 \(\displaystyle F^T(x) = x^2F(\frac{1}{x}) = 3x^2+x^2 \cdot \frac{1}{x}+x^2 \cdot \frac{4}{x^2} = 4+x+3x^2\)

那么由 \(F(x) = Q(x)G(x)+R(x)\)，考慮同乘 \(x^n\) 轉(zhuǎn)化為反轉(zhuǎn)形式：

\(\displaystyle x^nF(\frac{1}{x}) = (x^{n-m}Q(\frac{1}{x}))(x^{m}G(\frac{1}{x}))+x^{n-m+1}(x^{m-1}R(\frac{1}{x}))\)

\(\displaystyle F^T(x) = Q^T(x)G^T(x)+x^{n-m+1}R^T(x)\)，此時(shí)對(duì) \(x^{n-m+1}\) 取模消去 \(R^T(x)\)：

\(\displaystyle F^T(x) = Q^T(x)G^T(x) (\bmod x^{n-m+1})\)

\(\displaystyle Q^T(x) = F^T(x)(G^T(x))^{-1} (\bmod x^{n-m+1})\)

先全部反轉(zhuǎn)，求逆后相乘可得 \(Q^T(x)\)，反轉(zhuǎn)回來(lái)即得 \(Q(x)\)；最后將 \(Q(x)\) 與 \(G(x)\) 乘起來(lái)，與 \(F(x)\) 作差即得 \(R(x)\)

反轉(zhuǎn)可以直接用自帶的 reverse 函數(shù)

void modp(int *f, int *g, int n, int m){ // len_f = n, len_g = m, len_q = n-m, len_r = m-1
	static int g_[MAXN] = {0}, f_[MAXN] = {0}; int l = n-m+1;
	reverse(g, g+m+1); for (int i=0; i<l; i++) g_[i] = g[i]; reverse(g, g+m+1);
	reverse(f, f+n+1); for (int i=0; i<l; i++) f_[i] = f[i]; reverse(f, f+n+1);
	invp(g_, l); times(g_, f_, l, l); reverse(g_, g_+l); times(g, g_, n, n); 
	for (int i=0; i<n; i++){ g[i] = f[i]-g[i]; if (g[i] < 0) g[i] += MOD; }
	for (int i=0; i<l; i++) f[i] = g_[i]; for (int i=m; i<n; i++) g[i] = 0; 
	for (int i=0; i<(l<<1); i++) g_[i] = f_[i] = 0;
}

牛頓迭代

牛頓迭代是一種常用的分析方法，解決形如 "\(G(F(x)) = 0\)，求 \(F(x)\ (\bmod x^n)\)" 的問(wèn)題

這里不加證明地給出結(jié)論：\(\displaystyle \bm{F(x) = (F_{\triangle}(x)-\frac{G(F_{\triangle}(x))}{G'(F_{\triangle}(x))})\ (\bmod x^n)}\) (特別重要)

其中，\(F_{\triangle}(x)\) 表示倍增上一層的求解結(jié)果，\(G'(x)\) 表示多項(xiàng)式 \(G(x)\) 求導(dǎo)

\(n\) 次多項(xiàng)式 \(\displaystyle F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\) 求導(dǎo)定義為：\(\displaystyle F'(x) = a_1+2a_2x+3a_3x^2+\cdots+(n-1)a_{n-1}x^{n-1}\)
\(n\) 次多項(xiàng)式 \(\displaystyle F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\) 積分定義為：\(\displaystyle \int F(x)\text{dx} = C+a_0x+\frac{1}{2}a_1x^2+\cdots+\frac{1}{n+1}a_nx^{n+1}\)
實(shí)現(xiàn)時(shí)需要線性求逆元：

void Invinit(int n){ inv[1] = 1; for (int i=2; i<=n; i++) inv[i] = (1ll*inv[MOD%i]*(MOD-MOD/i))%MOD; }
void derip(int *f, int n){ for (int i=1; i<=n; i++) f[i-1] = (1ll*f[i]*i)%MOD; f[n] = 0; } // len = n, F'(x)
void sump(int *f, int n){ Invinit(n+1); for (int i=n+1; i>=1; i--) f[i] = (1ll*f[i-1]*inv[i])%MOD; f[0] = 0; } // len = n

多項(xiàng)式開(kāi)根

已知 \(B(x)\)，希望求使 \(A^2(x)-B(x) = 0\) 的 \(A(x)\)

法一： 構(gòu)造 \(G(A(x)) = A^2(x)-B(x)\)，則其導(dǎo)數(shù) \(G'(A(x)) = 2A(x)\) (將 \(B(x)\) 看作常量)，考慮直接套牛頓迭代：

\(\displaystyle A(x) = A_{\triangle}(x)-\frac{G(A_{\triangle}(x))}{G'(A_{\triangle}(x))} = A_{\triangle}(x)-\frac{A_{\triangle}^2(x)-B(x)}{2A_{\triangle}(x)} = \frac{A_{\triangle}^2(x)+B(x)}{2A_{\triangle}(x)}\)

那么直接倍增，求逆+乘法即可

邊界條件為：

• 若 \(B_0 = 1\)，則 \(A_0 = 1\)
• 反之，需要求解同余方程 \(A_0^2 \equiv B_0 (\bmod p)\)

法二： 轉(zhuǎn)化為后文要講的多項(xiàng)式快速冪，指數(shù)視為 \(2\) 關(guān)于 \(p\) 的逆元即可

void sqrtp(int *f, int n){ // rem->f, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; } invp(r, i);
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(w, (i<<1), 0);
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = f[j]+w[j]; if (w[j] >= MOD) w[j] -= MOD; } times(w, r, i, i); r[i] = 0;
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}

多項(xiàng)式 ln 與 exp

兩者均由麥克勞林級(jí)數(shù)定義：

• \(\displaystyle \ln(F(x)) = -\sum_{i=1} \frac{(1-F(x))^i}{i}\)
• \(\displaystyle \exp(F(x)) = \sum_{i=0} \frac{F^i(x)}{i!}\)

性質(zhì)：\(\ln(F(x))\) 與 \(\exp(F(x))\) 互為逆運(yùn)算 (重要)

對(duì)于 \(\ln(F(x))\)，我們肯定是不用定義式來(lái)求的；注意到 \(\displaystyle \ln'(x) = \frac{1}{x}\)，于是對(duì) \(\displaystyle \ln(A(x)) = B(x)\) 同時(shí)求導(dǎo)：

\(\displaystyle \ln'(A(x)) = \frac{A'(x)}{A(x)} = B'(x)\) (注意復(fù)合函數(shù)求導(dǎo)法則)

那么對(duì) \(\displaystyle \frac{A'(x)}{A(x)}\) 再積分一次即得 \(B(x)\)

void lnp(int *f, int n){ // len = n
	static int g[MAXN] = {0};
	for (int i=0; i<=n; i++) g[i] = f[i]; invp(g, n);
	derip(f, n); times(f, g, n, n); sump(f, n-1);
	for (int i=0; i<=n; i++) g[i] = 0;
}

對(duì)于 \(\exp(F(x))\)，由 \(\displaystyle e^{A(x)} = e^{\ln(B(x))} = B(x)\)，考慮構(gòu)造 \(G(B(x)) = \ln(B(x))-A(x)\)，直接牛頓迭代

\(\displaystyle G(B(x)) = B_{\triangle}(x)-\frac{G(B_{\triangle}(x))}{G'(B_{\triangle}(x))} = B_{\triangle}(x)-\frac{\ln(B_{\triangle}(x))-A(x)}{\frac{1}{B_{\triangle}(x)}} = B_{\triangle}(x)(1-\ln(B_{\triangle}(x))+A(x))\)

那么直接倍增，求 \(\ln\) + 乘法即可；邊界條件為 \(A_0 = 1\) 時(shí)，\(B_0 = 1\)

void exp(int *f, int n){ // len = n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = w[j]; } lnp(r, i);
		for (int j=0; j<i; j++){ r[j] = f[j]-r[j]; if (r[j] < 0) r[j] += MOD; } 
		r[0] = (r[0]+1)%MOD; times(w, r, i, i); r[i] = 0;
 	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}

多項(xiàng)式快速冪

給定 \(A(x)\) 與 \(k\)，求 \(A^k(x)\)；考慮對(duì) \(A^k(x)\) 取 \(\ln\)，將乘法變?yōu)榧臃?/strong>，最后 \(\exp\) 回去即可

void powp(int *f, int n, int k){ lnp(f, n); for (int i=0; i<=n; i++) f[i] = (1ll*f[i]*k%MOD+MOD)%MOD; exp(f, n); } // len = n

#include<bits/stdc++.h>
using namespace std;
#define db long double
#define ll long long
const db Pi = acos(-1);
const int MAXN = 4e6+5;

int n, m, p;
int tr[MAXN], f[MAXN], g[MAXN];

struct cp{
	db x, y; cp(db xx = 0, db yy = 0){ x = xx; y = yy; }
	cp operator + (const cp& b){ return cp(x+b.x, y+b.y); }
	cp operator - (const cp& b){ return cp(x-b.x, y-b.y); }
	cp operator * (const cp& b){ return cp(x*b.x-y*b.y, x*b.y+y*b.x); }
	cp operator / (const cp& b){ db t = (b.x*b.x+b.y*b.y); return cp((x*b.x+y*b.y)/t, (y*b.x-x*b.y)/t); }
};

void pre(int l){ for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0)); }
void fft(cp *f, int l, bool flg){ // 0->IDFT, 1->DFT
	pre(l); for (int i=0; i<l; i++) if (i<tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){ int nl = (i>>1); cp w1(cos(2*Pi/i), sin(2*Pi/i)); if (flg == 0) w1.y = -w1.y;
		for (int j=0; j<l; j+=i){ cp w(1, 0);
			for (int k=j; k<j+nl; k++){
				cp tmp = w*f[k+nl]; f[k+nl] = f[k]-tmp; f[k] = f[k]+tmp; w = w*w1;}}}
	if (flg == 0) for (int i=0; i<l; i++) f[i].x /= l, f[i].y /= l;
}
void mtt(int *f, int *g, int n, int m){
	static cp p1[MAXN], p2[MAXN], q[MAXN]; 
	int l = 0; for (l=1; l<=n+n; l<<=1);
	for (int i=0; i<l; i++){
		p1[i] = cp(f[i]>>15, f[i]&32767);
		p2[i] = cp(f[i]>>15, -(f[i]&32767));
		q[i] = cp(g[i]>>15, g[i]&32767);
	} fft(p1, l, 1); fft(p2, l, 1); fft(q, l, 1);
	for (int i=0; i<l; i++){ p1[i] = p1[i]*q[i]; p2[i] = p2[i]*q[i]; } fft(p1, l, 0); fft(p2, l, 0);
	for (int i=0; i<m; i++){
		ll a = (ll)floor((p1[i].x+p2[i].x)/2+0.49)%p;
		ll b = (ll)floor((p1[i].y+p2[i].y)/2+0.49)%p;
		ll c = (ll)floor((p1[i].y-b)+0.49)%p;
		ll d = (ll)floor((p2[i].x-a)+0.49)%p;
		f[i] = ((((a<<15)+(b+c))<<15)+d)%p;
		if (f[i] < 0) f[i] += p;
	}
}

int main(){
	scanf("%d%d%d", &n, &m, &p);
	for (int i=0; i<=n; i++) scanf("%d", &f[i]);
	for (int i=0; i<=m; i++) scanf("%d", &g[i]); mtt(f, g, max(n, m), n+m+1);
	for (int i=0; i<=n+m; i++) printf("%d ", f[i]);
	return 0;
}

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MAXN = 4e6+5;
const ll MOD = 998244353, G = 3, INVG = 332748118, INV2 = 499122177;

int n, k;
int tr[MAXN], f[MAXN], g[MAXN], inv[MAXN];

int read(){
	bool f=1;ll x=0;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-') f=!f;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);x=x%MOD;ch=getchar();}
	x=(f?x:-x);return (int)x;
}

ll QuickPow(ll a, ll b){
	ll res = 1;
	while (b){
		if (b&1) res = (res*a)%MOD;
		a = (a*a)%MOD;
		b >>= 1;
	} return res;
}

void Invinit(int n){ inv[1] = 1; for (int i=2; i<=n; i++) inv[i] = (1ll*inv[MOD%i]*(MOD-MOD/i))%MOD; }
void pre(int l){ for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0)); }
void ntt(int *g, int l, bool flg){ // 0->IDFT, 1->DFT, l>n
	pre(l); static int f[MAXN] = {0};
	for (int i=0; i<l; i++) f[i] = g[i]; for (int i=0; i<l; i++) if (i<tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){ int nl = (i>>1), w1 = QuickPow((flg?G:INVG), (MOD-1)/i);
		for (int j=0; j<l; j+=i){ ll w = 1;
			for (int k=j; k<j+nl; k++){
				int tmp = (w*f[k+nl])%MOD;
				f[k+nl] = f[k]-tmp; if (f[k+nl]<0) f[k+nl]+=MOD;
				f[k] = f[k]+tmp; if (f[k]>=MOD) f[k]-=MOD; w = (w*w1)%MOD;}}}
	if (flg == 0){ int inv = QuickPow(l, MOD-2); for (int i=0; i<l; i++) g[i] = (1ll*f[i]*inv)%MOD; }
	else{ for (int i=0; i<l; i++) g[i] = f[i]; }
	for (int i=0; i<(l<<1); i++) f[i] = 0;           
}
void mulp(int *f, int *g, int l){ for (int i=0; i<l; i++) f[i] = (1ll*f[i]*g[i])%MOD; } // l > n
void times(int *f, int *g, int n, int m){ // max{len} = n, mod x^m
	static int tmp[MAXN] = {0};
	int l = 1; for (l=1; l<=n+n; l<<=1); for (int i=0; i<l; i++) tmp[i] = g[i];
	ntt(f, l, 1); ntt(tmp, l, 1); mulp(f, tmp, l); ntt(f, l, 0);
	for (int i=0; i<l; i++) tmp[i] = 0; for (int i=m; i<l; i++) f[i] = 0;
}
void invp(int *f, int n){ // len = n, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}, tmp[MAXN] = {0}; w[0] = QuickPow(f[0], MOD-2);
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; }
		for (int j=0; j<i; j++) tmp[j] = f[j];
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(tmp, (i<<1), 1); mulp(w, tmp, (i<<1)); ntt(w, (i<<1), 0); 
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = r[j]-w[j]; if (w[j] < 0) w[j] += MOD; }
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = tmp[i] = 0;
}
void sqrtp(int *f, int n){ // rem->f, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; } invp(r, i);
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(w, (i<<1), 0);
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = f[j]+w[j]; if (w[j] >= MOD) w[j] -= MOD; } times(w, r, i, i); r[i] = 0;
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}
void modp(int *f, int *g, int n, int m){ // len_f = n, len_g = m, len_q = n-m, len_r = m-1
	static int g_[MAXN] = {0}, f_[MAXN] = {0}; int l = n-m+1;
	reverse(g, g+m+1); for (int i=0; i<l; i++) g_[i] = g[i]; reverse(g, g+m+1);
	reverse(f, f+n+1); for (int i=0; i<l; i++) f_[i] = f[i]; reverse(f, f+n+1);
	invp(g_, l); times(g_, f_, l, l); reverse(g_, g_+l); times(g, g_, n, n); 
	for (int i=0; i<n; i++){ g[i] = f[i]-g[i]; if (g[i] < 0) g[i] += MOD; }
	for (int i=0; i<l; i++) f[i] = g_[i]; for (int i=m; i<n; i++) g[i] = 0; 
	for (int i=0; i<(l<<1); i++) g_[i] = f_[i] = 0;
}
void derip(int *f, int n){ for (int i=1; i<=n; i++) f[i-1] = (1ll*f[i]*i)%MOD; f[n] = 0; } // len = n, F'(x)
void sump(int *f, int n){ Invinit(n+1); for (int i=n+1; i>=1; i--) f[i] = (1ll*f[i-1]*inv[i])%MOD; f[0] = 0; } // len = n
void lnp(int *f, int n){ // len = n
	static int g[MAXN] = {0};
	for (int i=0; i<=n; i++) g[i] = f[i]; invp(g, n);
	derip(f, n); times(f, g, n, n); sump(f, n-1);
	for (int i=0; i<=n; i++) g[i] = 0;
}
void exp(int *f, int n){ // len = n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = w[j]; } lnp(r, i);
		for (int j=0; j<i; j++){ r[j] = f[j]-r[j]; if (r[j] < 0) r[j] += MOD; } 
		r[0] = (r[0]+1)%MOD; times(w, r, i, i); r[i] = 0;
 	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}
void powp(int *f, int n, int k){ lnp(f, n); for (int i=0; i<=n; i++) f[i] = (1ll*f[i]*k%MOD+MOD)%MOD; exp(f, n); } // len = n

int main(){		
	return 0;
}