算法學習筆記：多項式

多項式基礎

為了與冪次相對應，下文中系數均采用 \(\text{0-based}\) 下標計數。

多項式的表示方法

我們習慣于將一個 \(n\) 次多項式 \(f(x)\) 表示為 \(f(x)=\sum_{i=0}^{n-1}a_ix^i\)。那么有沒有其他表示方法呢？

我們不妨考慮代入 \(n\) 個互不相同的 \(x_i\) 得到的 \(n\) 個點值 \((x_i,y_i)\)，很自然地猜測它們可以唯一確定這個 \(n\) 次多項式。可以發現 \(n\) 個點值對應一個關于系數的 \(n\) 元一次方程，顯然有唯一解，因此這 \(n\) 個點值可以唯一確定這個 \(n\) 次多項式。我們稱這種表示方法為點值表示法。而我們常用的 \(f(x)=\sum_{i=0}^{n-1}a_ix^i\) 的表示方法則被稱作系數表示法。

既然這兩種表示法都可以唯一確定一個多項式，那么必然可以相互轉化。具體地，系數表示法轉化為點值表示法的過程就是求值，點值表示法轉化為系數表示法的過程就是插值。

離散卷積

在這里較為不嚴謹地給出離散卷積的概念。

給定函數 \(f,g\)，形如

\[h(n)=\sum_{i\otimes j=n}f(i)g(n-i) \]

的式子就是 \(f\) 與 \(g\) 的卷積，其中 \(\otimes\) 是某種運算。可以將其寫作 \(h=f*g\)。

考慮多項式乘法，給出多項式 \(F(x),G(x)\)，令 \(f(n)=[x^n]F(x)\)，\(g(n)=[x^n]G(x)\) 則

\[[x^n]f(x)g(x)=\sum_{i+j=n}f(i)g(j) \]

因此多項式乘法本質上就是其系數序列的加法卷積。

拉格朗日插值

我們簡單研究一下插值。一個顯然的方法就是 \(O(n^3)\) 高斯消元，而下面介紹的拉格朗日插值可以做到 \(O(n^2)\) 插值。

現在給出 \(n\) 個點值 \((x_i,y_i)\)，保證 \(x_i\) 互不相同，我們想要找到對應的多項式 \(f(x)\) 的系數表示法。

拉格朗日插值的思想就是構造 \(n\) 個 \(n\) 次多項式 \(g_i(x)\)，滿足 \(g_i(x_i)=y_i\)，\(g_i(x_j)=0(j\neq i)\)，那么所求 \(f(x)\) 即為 \(\sum_{i=0}^{n-1}g_i(x)\)。考慮如何構造 \(g_i(x)\)。顯然 \(x_j(j\neq i)\) 給出了 \(g_i(x)\) 的 \(n-1\) 個根，因此它必然包含 \(\prod_{j\neq i}(x-x_j)\) 的因式，而我們還要求 \(g_i(x_i)=y_i\)，調整系數，乘上 \(\frac{y_i}{g_i(x_i)}\) 即可。綜上

\[g_i(x)=y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \]

因此

\[f(x)=\sum_{i=0}^{n-1}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \]

實現時可以 \(O(n^2)\) 求出 \(\prod_{j=0}^{n-1}(x-x_j)\) 的各項系數，對于每個 \(i\) 直接 \(O(n)\) 除掉 \(x-x_i\)，再調整系數，最后對所有 \(i\) 系數相加即可。時間復雜度為 \(O(n^2)\)。

Luogu P4781 【模板】拉格朗日插值只要求單點插值，直接帶入拉格朗日插值公式計算即可，依然是 \(O(n^2)\) 的。

橫坐標是連續整數的快速單點插值

若 \(x_i\) 為連續整數，則我們可以給出 \(O(n)\) 的單點插值的方法。

設 \(x_i=k+i\)，則拉格朗日插值公式變為

\[f(x)=\sum_{i=0}^{n-1}y_i\prod_{j\neq i}\frac{x-j-k}{i-j} \]

注意到分子是經典的去除一個單點的形式，考慮預處理前后綴積，即令 \(pre_i=\prod_{j=0}^i(x-j-k)\)，\(suf_i=\prod_{j=i}^{n-1}(x-j-k)\)。而對于分母，\(j<i\) 部分的乘積為 \(i!\)，\(j>i\) 部分的乘積為 \((-1)^{n-i-1}(n-i-1)!\)。容易得到

\[f(x)=\sum_{i=0}^{n-1}\frac{y_ipre_{i-1}suf_{i+1}}{(-1)^{n-i-1}i!(n-i-1)!} \]

預處理階乘逆元即可 \(O(n)\) 計算。

快速傅里葉變換

快速傅里葉變換（FFT）是多項式的基石。這里默認讀者具有良好的線性代數與復數知識儲備。

離散傅里葉變換

離散傅里葉變換（DFT）是傅里葉變換在時域和頻域上都呈離散的形式，將信號的時域采樣變換為其 DTFT 的頻域采樣。在 OI 中，其恰好能用來將多項式的系數表示法轉化為點值表示法。

序列 \(\{x_i\}_{i=0}^{n-1}\) 的 DFT 為

\[X_i=\sum_{j=0}^{n-1}x_j\omega_n^{-ij} \]

若將 \(\{x_i\}_{i=0}^{n-1}\) 視為多項式 \(f(x)\) 的系數序列，則我們有

\[X_i=f(\omega_n^{-ij}) \]

通常以符號 \(\mathcal{F}\) 表示該變換，即 \(\hat{x}=\mathcal{F}x\)。而它的逆離散傅里葉變換為

\[x_i=\frac{1}{n}\sum_{j=0}^{n-1}X_j\omega_n^{ij} \]

也記作 \(x=\mathcal{F}^{-1}\hat{x}\)。

容易發現，離散傅里葉變換可以看作一個 \(n-1\) 次多項式在 \(n\) 個單位根處求值。

分治實現 FFT

FFT 是一種高效實現 DFT 的算法，不過其中的一些細節又與 DFT 略有區別，后文將會提及。

我們的目標就是，在可接受的時間復雜度下，快速求出一個 \(n-1\) 次多項式的 \(n\) 個點值。

FFT 的核心思想就是分治。我們考慮 \(f(x)\) 拆成一個偶函數 \(f_e(x)\) 和一個奇函數 \(f_o(x)\) 之和，也就是令 \(f_e(x)=a_0x^0+a_2x^2+\cdots\)，\(f_o(x)=a_1x^1+a_3x^3+\cdots\)，此時

\[f(x)=f_e(x)+f_o(x) \]

這樣并不能得到很好的復雜度，因為 \(f_e(x)\) 和 \(f_o(x)\) 依然是 \(n-1\) 次的多項式。但注意到項數減半，考慮換元，令 \(t=x^2\)，則 \(f_e(x)=a_0t^0+a_2t^1+\cdots\)，\(f_o(x)=x(a_1t^0+a_3t^1+\cdots)\)，換句話說

\[f(x)=f_e(x^2)+xf_o(x^2) \]

此時 \(f_e(x)\) 和 \(f_o(x)\) 降成了 \(\frac{n}{2}\) 次多項式，我們做到了規模減半。

接下來，由偶函數和奇函數的性質，考慮代入 \(\frac{n}{2}\) 對互為相反數的 \((x_i,-x_i)\)，因為這樣會使得它們的 \(f_e(x^2)\) 和 \(f_o(x^2)\) 值是相同的。而我們需要保證子問題也能找到這樣成對的相反數，也即要找到 \(n\) 個互不相同的 \(x_i\)，它們的 \(n\) 次冪相等，顯然所有 \(n\) 次單位根就滿足這個性質，因此我們代入 \(\omega_n^k\) 和 \(-\omega_n^k=\omega_n^{k+\frac{n}{2}}\)：

\[\begin{align*} f(\omega_n^k)&=f_e(\omega_n^{2k})+\omega_n^kf_o(\omega_n^{2k}) \\ f(\omega_n^{k+\frac{n}{2}})&=f_e(\omega_n^{2k})-\omega_n^kf_o(\omega_n^{2k}) \end{align*} \]

再利用 \(\omega_{n}^{2k}=\omega_{\frac{n}{2}}^k\)：

\[\begin{align*} f(\omega_n^k)&=f_e(\omega_{\frac{n}{2}}^k)+\omega_n^kf_o(\omega_{\frac{n}{2}}^k) \\ f(\omega_n^{k+\frac{n}{2}})&=f_e(\omega_{\frac{n}{2}}^k)-\omega_n^kf_o(\omega_{\frac{n}{2}}^k) \end{align*} \]

我們以 \(n=1\) 作為遞歸邊界，在每個子問題中對系數 \(a_i\) 按 \(i\) 的奇偶性左右分組，再遞歸分治處理，最后統計答案，就可以分治實現 FFT 了！

\(T(n)=T(\frac{n}{2})+O(n)\)，由主定理，時間復雜度為 \(O(n\log{n})\)。

注意上述過程都是在 \(n=2^k(k\in\mathbb{N})\) 的前提下進行的。因此具體實現時，要把 \(n\) 向上補成 \(2^k\)。

我們還需要支持復數運算，可以手寫結構體，也可以使用 STL complex。

分治實現的代碼：

void FFT(complex<double> *a, int n) {
    if (n == 1) return;
    int mid = n >> 1;
    for (int i = 0; i < n; i += 2)
        tmp[i >> 1] = a[i], tmp[(i >> 1) + mid] = a[i + 1];
    for (int i = 0; i < n; ++i) a[i] = tmp[i];
    FFT(a, mid, tp), FFT(a + mid, mid, tp);
    w[0] = { 1, 0 }, w[1] = { cos(PI * 2 / n), sin(PI * 2 / n) };
    for (int i = 2; i < mid; ++i) w[i] = w[i - 1] * w[1];
    for (int i = 0; i < mid; ++i)
        tmp[i] = a[i] + w[i] * a[i + mid],
        tmp[i + mid] = a[i] - w[i] * a[i + mid];
    for (int i = 0; i < n; ++i) a[i] = tmp[i];
}

倍增實現 FFT/蝴蝶變換

分治實現 FFT 需要遞歸處理，效率較慢，我們希望找到非遞歸的實現方式。

考察 FFT 對系數位置的變換。對于一個規模為 \(n\) 的問題，位置 \(p\) 上的系數會在 \(p\) 為偶數時被分到左側，為奇數時分到右側。進一步地，位置 \(p\) 在第 \(i\) 次遞歸的方向決定了它最終下標二進制表示下第 \(\log_2n-i\) 位的值，這等價于 \(p\) 最終下標二進制表示下第 \(\log_2n-i\) 位的取值，就是 \(p\) 二進制表示下第 \(i\) 位的取值。因此，系數 \(a_p\) 在 \(\log_2n\) 次遞歸后的下標就是 \(p\) 的二進制表示反轉的結果，不妨記作 \(rev_p\)。

顯然 \(rev_p\) 可以 \(O(n)\) 遞推：

\[rev_p=\left\lfloor\frac{rev_{\left\lfloor\frac{p}{2}\right\rfloor}}{2}\right\rfloor+(p\bmod{2})\frac{n}{2} \]

我們在一開始就做一次位逆序置換，即令 \(a_i\leftarrow a_{rev_i}\)。然后考慮如何自底向上合并。假設我們已經完成了所有規模為 \(n\) 的變換，而想要轉而完成所有規模為 \(2n\) 的變換，回到先前的分治式子

\[\begin{align*} f(\omega_{2n}^k)&=f_e(\omega_n^k)+\omega_{2n}^kf_o(\omega_n^k) \\ f(\omega_{2n}^{k+n})&=f_e(\omega_n^k)-\omega_{2n}^kf_o(\omega_n^k) \end{align*} \]

此時，\(f_e(\omega_n^k)\) 和 \(f_o(\omega_n^k)\) 的值分別存在下標為 \(k\) 和 \(k+n\) 的位置，而我們想求出的 \(f(\omega_{2n}^k)\) 和 \(f(\omega_{2n}^{k+n})\) 恰好將分別存在下標為 \(k\) 和 \(k+n\) 的位置，所以我們直接在原數組上覆寫即可。具體來說，令 \(x=a_k\)，\(y=\omega_{2n}^ka_{k+n}\)，然后令 \(a_k\leftarrow x+y\)，\(a_{k+n}\leftarrow x-y\) 就完成了一對值的求解。這個過程就稱為蝴蝶變換。

我們理一下算法的過程：先對系數數組做位逆序置換，然后枚舉問題規模 \(2k(1\leq k<n)\)，枚舉每個子問題 \(2ki(0\leq 2ki<n)\)，枚舉子問題中的對應位置 \(x,y(x=2ki+j,y=2ki+j+k)\)，記 \(tx=a_x\)，\(ty=\omega_{2k}^ja_y\)，令 \(a_x\leftarrow tx+ty\)，\(a_y\leftarrow tx-ty\) 即可。

由 \(rev\) 的對稱性，做位逆序置換時只需在 \(i<rev_i\) 時 \(\operatorname{swap}(a_i,a_{rev_i})\)。

時間復雜度還是 \(O(n\log{n})\) 的，但相比分治實現效率高了不少。

蝴蝶變換的代碼：

void FFT(complex<double> *a, int n) {
    int mid = n >> 1;
    for (int i = 1; i < n; ++i) {
        rev[i] = (rev[i >> 1] >> 1) + (i & 1 ? mid : 0);
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    for (int k = 1; k < n; k <<= 1) {
        w[0] = { 1, 0 }, w[1] = { cos(PI / k), sin(PI / k) };
        for (int i = 2; i < k; ++i) w[i] = w[i - 1] * w[1];
        for (int i = 0; i < n; i += k << 1) for (int j = 0; j < k; ++j) {
            complex<double> x = a[i | j], y = w[j] * a[i | j | k];
            a[i | j] = x + y, a[i | j | k] = x - y;
        }
    }
}

多項式求值和線性代數

在講解快速傅里葉逆變換之前，我們先講一下多項式求值和線性代數之間的關系。

其實多項式求值很顯然可以用矩陣乘法描述。例如，把多項式 \(f(x)=\sum_{i=0}^{n-1}a_ix^i\) 對 \(x_0,x_1,\cdots,x_{n-1}\) 求值，就可以表示為

\[\begin{bmatrix} x_0^0 & x_0^1 & \cdots & x_0^{n-1} \\ x_1^0 & x_1^1 & \cdots & x_1^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n-1}^0 & x_{n-1}^1 & \cdots & x_{n-1}^{n-1} \end{bmatrix} \times \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}= \begin{bmatrix} f(x_0) \\ f(x_1) \\ \vdots \\ f(x_{n-1}) \end{bmatrix} \]

左側的矩陣就是范德蒙德矩陣。

那么 FFT 也可以用這種方式表示：

\[\begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots & (\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \times \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}= \begin{bmatrix} f(x_0) \\ f(x_1) \\ \vdots \\ f(x_{n-1}) \end{bmatrix} \]

快速傅里葉逆變換

設 FFT 中的范德蒙德矩陣為 \(T\)，即

\[T=\begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots & (\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \]

由前面的矩陣乘法表示，容易想到逆變換就是左乘上這個矩陣的逆矩陣 \(T^{-1}\)。

考慮利用拉格朗日插值公式求逆。觀察拉格朗日插值公式，容易得出

\[(T^{-1})_{ij}=[x^i]\prod_{k\neq j}\frac{x-\omega_n^k}{\omega_n^j-\omega_n^k} \]

分子和分母都是

\[g(x)=\frac{\prod_{i=0}^{n-1}(x-\omega_n^i)}{x-\omega_n^j} \]

的形式，考慮研究其性質。由代數基本定理，\(\prod_{i=0}^{n-1}(x-\omega_n^i)=x^n-1\)，因此化為

\[g(x)=\frac{x^n-1}{x-\omega_n^j} \]

模擬短除法，可以得到

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

因此 \([x^i]\prod_{k\neq j}(x-\omega_n^k)=(\omega_n^j)^{n-1-i}\)，\(g(\omega_n^j) =n(\omega_n^j)^{n-1}\)，于是

\[(T^{-1})_{ij}=\frac{(\omega_n^j)^{n-1-i}}{n(\omega_n^j)^{n-1}}=\frac{\omega_n^{-i}}{n}=\frac{\overline{\omega_n^i}}{n} \]

即我們所求的逆矩陣

\[T^{-1}=\begin{bmatrix} \left(\overline{\omega_n^0}\right)^0 & \left(\overline{\omega_n^0}\right)^1 & \cdots & \left(\overline{\omega_n^0}\right)^{n-1} \\ \left(\overline{\omega_n^1}\right)^0 & \left(\overline{\omega_n^1}\right)^1 & \cdots & \left(\overline{\omega_n^1}\right)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ \left(\overline{\omega_n^{n-1}}\right)^0 & \left(\overline{\omega_n^{n-1}}\right) & \cdots & \left(\overline{\omega_n^{n-1}}\right)^{n-1} \end{bmatrix} \]

所以 IFFT 的過程和 FFT 沒有很大的區別。我們只需要將 FFT 中代入的單位根變為原來的共軛，并在最后對每一項系數除以 \(n\) 即可。

為了方便，我們可以把 FFT 與 IFFT 放在一個函數里，用 \(tp=1\) 表示是 FFT，用 \(tp=-1\) 表示是 IFFT，那么這個 \(tp\) 就恰好可以乘到單位根的虛部上來轉成共軛。

FFT/IFFT 的代碼：

void FFT(complex<double> *a, int n, int tp = 1) {
    int mid = n >> 1;
    for (int i = 1; i < n; ++i) {
        rev[i] = (rev[i >> 1] >> 1) + (i & 1 ? mid : 0);
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    for (int k = 1; k < n; k <<= 1) {
        w[0] = { 1, 0 }, w[1] = { cos(PI / k), tp * sin(PI / k) };
        for (int i = 2; i < k; ++i) w[i] = w[i - 1] * w[1];
        for (int i = 0; i < n; i += k << 1) for (int j = 0; j < k; ++j) {
            complex<double> x = a[i | j], y = w[j] * a[i | j | k];
            a[i | j] = x + y, a[i | j | k] = x - y;
        }
    }
    if (tp == -1) for (int i = 0; i < n; ++i) a[i] /= n;
}

DFT 和 FFT

回顧前面的內容，可以總結出 DFT 和 FFT 之間的幾點差別：

FFT 要求 \(n\) 是 \(2\) 的整數次冪，DFT 則不做要求。
DFT 與 FFT 代入單位根的順序相反。

快速數論變換

數論變換（NTT）是 DFT 在數論基礎上的實現，快速數論變換（FNTT）則是 FFT 在數論基礎上的實現。更具體地，FNTT 是在模 \(p\) 意義下的整數域進行的 FFT。

對次數為 \(n-1\) 的多項式進行變換，我們需要找到 \(n\) 次單位根 \(a\) 滿足 \(\delta_p(a)=n\)，這要求模數 \(p\) 為一個質數，且 \(2^k\mid p-1\)，其中 \(2^k\leq\) 最大的可能 NTT 長度。

根據原根的性質，\(\delta_p(g)=p-1\)，即 \(g\) 的 \(0\sim p-2\) 次冪互不相同，因此 \(g^k\) 和 \(\omega_{p-1}^k(0\leq k<p-1)\) 形成的域是同構的。進一步地，\(\omega_{n}\) 等價于 \(g^{\frac{p-1}{n}}\)，其 \(0\sim n-1\) 次冪互不相同，滿足我們做 FFT 時所需的性質。

因此我們只需要將 FFT 時的 \(\omega_n\) 換成 \(g^{\frac{p-1}{n}}\)，\(\overline{\omega_n}\) 換成 \((g^{-1})^\frac{p-1}{n}\) 就變成了 NTT 啦！由于沒有了浮點數運算，NTT 的效率比 FFT 高不少。

下面列出一些 NTT 常見大模數：

\(998244353=7\times 17\times 2^{23}+1\)，有原根 \(3\)。
\(1004535809=479\times 2^{21}+1\)，有原根 \(3\)。
\(469762049=7\times 2^{26}+1\)，有原根 \(3\)。

NTT 的代碼（\(p=998244353\)）：

const int g1 = 3, g2 = (MOD + 1) / g1;
void NTT(int *a, int n, int tp = 1) {
    int mid = n >> 1;
    for (int i = 1; i < n; ++i) {
        rev[i] = (rev[i >> 1] >> 1) + (i & 1 ? mid : 0);
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    for (int k = 1; k < n; k <<= 1) {
        w[0] = 1, w[1] = qpow(tp == 1 ? g1 : g2, (MOD - 1) / (k << 1));
        for (int i = 2; i < k; ++i) w[i] = 1ll * w[i - 1] * w[1] % MOD;
        for (int i = 0; i < n; i += k << 1) for (int j = 0; j < k; ++j) {
            int x = a[i | j], y = 1ll * w[j] * a[i | j | k] % MOD;
            a[i | j] = (x + y) % MOD, a[i | j | k] = (x - y + MOD) % MOD;
        }
    }
    if (tp == -1) {
        int iv = qpow(n, MOD - 2);
        for (int i = 0; i < n; ++i) a[i] = 1ll * a[i] * iv % MOD;
    }
}

FFT/NTT 的應用

多項式乘法

Luogu P3803 【模板】多項式乘法（FFT）

題意：給出一個次數為 \(n\) 的多項式 \(f(x)\) 和一個次數為 \(m\) 的多項式 \(g(x)\)，求它們的卷積。\(1\leq n,m\leq 10^6\)，\(0\leq [x^i]f(x),[x^i]g(x)\leq 9\)。

下面視 \(n,m\) 同階。按照卷積的定義暴力做顯然是 \(O(n^2)\) 的。不妨從點值表示法的角度考慮，容易發現兩個多項式的卷積的點值表示即兩個多項式的點值對應相乘。所以我們用 FFT 把 \(f(x),g(x)\) 轉成點值表示法，\(O(n)\) 點值相乘，再 IFFT 轉回系數表示法即可。時間復雜度為 \(O(n\log{n})\)。

進一步地，我們發現卷積后的系數不超過 \(nV^2\)，其中 \(V\) 為多項式系數的值域，所以我們可以把 FFT 改成對 \(998244353\) 取模的 NTT，從而優化常數。

代碼中的函數calc_len 和 clr 在后面的代碼中會用到：

int calc_len(int n) { return 1 << (int)ceil(log2(n)); }
void clr(int *a, int n) { memset(a, 0, n << 2); }

cin >> n >> m; ++n, ++m;
for (int i = 0; i < n; ++i) cin >> f[i];
for (int i = 0; i < m; ++i) cin >> g[i];
int len = calc_len(n + m - 1);
NTT(f, len), NTT(g, len);
for (int i = 0; i < len; ++i) f[i] = 1ll * f[i] * g[i] % MOD;
NTT(f, len, -1);
for (int i = 0; i < n + m - 1; ++i) cout << f[i] << " \n"[i == n + m - 2];

分治 FFT/NTT

半在線卷積

Luogu P4721 【模板】分治 FFT

題意：給定長度為 \(n-1\) 的序列 \(g\)，求出一個長度為 \(n\) 的序列 \(f\)，滿足

\[f_i= \begin{cases} 1, & i=0 \\ \sum_{j=1}^if_{i-j}g_j,& 0<i<n \end{cases} \]

答案對 \(998244353\) 取模。\(2\leq n\leq 10^5\)。

對于一類形如 \(f_i=\sum_{j=1}^if_{i-j}g_j\) 的遞推式，\(f_i\) 的值依賴于 \(f[0,i)\) 的值，我們稱其為半在線卷積問題。

這是一個在線問題，考慮用 CDQ 分治將其轉化為離線問題。具體來說，假設我們求出了 \(f[l,mid]\)，現在想要計算其對 \(f[mid+1,r]\) 的貢獻。對于 \(f_i(i\in[mid+1,r])\)，容易計算得到 \(f[l,mid]\) 對它的貢獻為 \(\sum_{j=l}^{mid}f_{j}g_{i-j}\)，這就變成了很簡單的卷積問題了，我們對 \(f[l,mid]\) 和 \(g[0,r-l]\) 做卷積即可得到貢獻序列。

代碼：

void solve(int l, int r) {
    if (l == r) return;
    static int a[N], b[N];
    int mid = (l + r) >> 1;
    solve(l, mid);
    copy(f + l, f + mid + 1, a), copy(g, g + r - l + 1, b);
    int len = calc_len(r - l + 1);
    NTT(a, len), NTT(b, len);
    for (int i = 0; i < len; ++i) a[i] = 1ll * a[i] * b[i] % MOD;
    NTT(a, len, -1);
    for (int i = mid + 1; i <= r; ++i) f[i] = (f[i] + a[i - l]) % MOD;
    clr(a, len), clr(b, len);
    solve(mid + 1, r);
}

單項式乘積

給出 \(n\) 個一次多項式 \(f_i(x)\)，我們希望計算它們的乘積多項式 \(\prod_{i=1}^n{f_i(x)}\)。

如果我們從左往右 NTT 合并，復雜度會是 \(O(n^2\log{n})\) 的。

考慮分治計算，令 \(g_{l,r}(x)=\prod_{i=l}^rf_i(x)\)，于是顯然 \(g_{l,r}(x)=g_{l,mid}(x)g_{mid+1,r}(x)\)，遞歸合并即可。\(T(n)=2T(\frac{n}{2})+O(n\log{n})\)，由主定理，時間復雜度為 \(O(n\log^2{n})\)。

乍一看這似乎不太合理，為什么把從左到右合并改成分治合并就能得到更優的時間復雜度呢？究其原因，從左到右合并就是不斷求一個 \(k\) 次多項式和一個單項式的乘積的過程，NTT 很浪費時間，甚至不如根據定義暴力合并。而分治則平衡了左右兩側的多項式的規模，使得復雜度得到了優化。

多項式牛頓迭代

泰勒級數和麥克勞林級數

這里簡單介紹學習牛頓迭代的一些基礎知識。讀者需要有簡單的微積分基礎。

對于一個光滑的函數 \(f(x)\)，我們有

\[f(x)=\sum_{n=0}^{+\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n \]

上式稱為 \(f(x)\) 在 \(x=a\) 處的泰勒展開式。特別地，在 \(x=0\) 處的泰勒展開式被稱為麥克勞林展開式。

這里列出一些常見的麥克勞林展開式，后文會用到一部分：

\[\begin{align*} (1+x)^m&=\sum_{n=0}^{+\infty}\binom{m}{n}x^n\\ \ln(1-x)&=-\sum_{n=1}^{+\infty}\frac{x^n}{n} \\ \ln(1+x)&=-\sum_{n=1}^{+\infty}\frac{(-1)^nx^n}{n} \\ \exp(x)&=-\sum_{n=0}^{+\infty}\frac{x^n}{n!} \end{align*} \]

一般的牛頓迭代

一般的牛頓迭代按以下方式求出方程 \(f(x)=0\) 的根：

設我們得到的近似根為 \(x_0\)。
將 \(f(x)\) 在 \(x=x_0\) 處做泰勒展開。
取前兩項的系數來近似成 \(f(x)\)，即 \(f(x)=f(x_0)+f'(x_0)(x-x_0)=0\)。
解得 \(x=x_0-\frac{f(x_0)}{f'(x_0)}\)，將其作為新的近似根繼續迭代下去。

容易發現，一般的牛頓迭代的本質上就是不斷求切線的過程。

多項式牛頓迭代

我們可以用牛頓迭代的方式解決這樣一類問題：給出關于多項式 \(f(x)\) 的函數 \(g(f(x))\)，求一個多項式 \(f(x)\) 使得 \(g(f(x))\equiv 0\pmod{x^n}\) 成立。

與一般的牛頓迭代不同，由于多項式的一些特性，用牛頓迭代求解這類問題并無精度問題。

假設我們求出了 \(f_0(x)\) 使得 \(g(f_0(x))\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)，我們把 \(g(f(x))\) 在 \(f(x)=f_0(x)\) 處做泰勒展開

\[g(f(x))\equiv\sum_{n=0}^{+\infty}\frac{g^{(n)}(f_0(x))}{n!}(f(x)-f_0(x))^n\equiv 0\pmod{x^n} \]

注意到，\(f_0(x)\) 和 \(f(x)\) 的前 \(\left\lceil\frac{n}{2}\right\rceil\) 項相同，即

\[\begin{align*} f(x)-f_0(x)\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}} \Leftrightarrow (f(x)-f_0(x))^2\equiv 0\pmod{x^n} \end{align*} \]

所以對于 \(n\geq 2\)，\(\frac{g^{(n)}(f_0(x))}{n!}(f(x)-f_0(x))^n\equiv 0\pmod{x^n}\)，于是我們的泰勒級數依然只需保留前兩項：

\[g(f(x))\equiv g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\equiv 0\pmod{x^n} \]

類似解得

\[f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n} \]

我們只需要在初始時給出 \(\bmod{x^1}\) 的解就能進行多項式牛頓迭代了。

有一些細節需要注意：

由于 \(g(f_0(x))\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)，\(g'(f_0(x))\) 中次數 \(\geq\left\lceil\frac{n}{2}\right\rceil\) 的項會被 \(\bmod{x^n}\) 截斷，所以 \(g'(f_0(x))\) 只需要做到 \(\bmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\) 的精度。
\(g'(f_0(x))\) 是將 \(f_0(x)\) 視作主元，而非 \(x\)，即分母要求的是 \(\frac{\textw0obha2h00g(f_0(x))}{\textw0obha2h00f_0(x)}\)。

顯然多項式牛頓迭代的精度成平方增長，這通常允許我們使用倍增求解。讀者可以結合接下來的一些例子體會。

多項式乘法逆

Luogu P4238 【模板】多項式乘法逆

題意：給定一個 \(n-1\) 次多項式 \(f(x)\)，求出 \(g(x)\) 使得 \(f(x)g(x)\equiv 1\pmod{x^n}\)，系數對 \(998244353\) 取模。\(1\leq n\leq 10^5\)。

這里給出兩種求解方法：

方法一：平方法

顯然 \(\bmod{x^1}\) 的解就是 \(a_0^{-1}\bmod{998244353}\)。接下來假設我們求出了 \(g_0(x)\) 使得 \(g_0(x)f(x)\equiv 1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)，來快樂地推一下式子：

\[\begin{align*} (g(x)-g_0(x))^2 & \equiv 0 & \pmod{x^n} \\ g(x)^2-2g(x)g_0(x)+g_0(x)^2 & \equiv 0 & \pmod{x^n} \\ g(x)-2g_0(x)+f(x)g_0(x)^2 & \equiv 0 & \pmod{x^n} \\ g(x) & \equiv 2g_0(x)-f(x)g_0(x)^2 & \pmod{x^n} \end{align*} \]

那么倍增求解到最小的 \(k\) 使得 \(2^k\geq n\) 即可。\(T(n)=T(\frac{n}{2})+O(n\log{n})\)，由主定理，時間復雜度為 \(O(n\log{n})\)。

方法二：牛頓迭代

我們相當于要讓 \(h(g(x))\equiv g(x)f(x)-1\equiv 0\)。直接代入多項式牛頓迭代的公式求解：

\[\begin{align*} g(x)& \equiv g_0(x)-\frac{f(x)g_0(x)-1}{f(x)} \\ & \equiv g_0(x)-g_0(x)(f(x)g_0(x)-1) \\ & \equiv 2g_0(x)-f(x)g_0(x)^2\pmod{x^n} \end{align*} \]

這和前面推的式子是一樣的，同樣 \(O(n\log{n})\) 倍增求解即可。

回顧多項式牛頓迭代的推導過程，讀者可以體會平方法本質上與牛頓迭代的相似性。

倍增法的細節

倍增法有一些很重要的細節，不加注意就會出現問題：

注意考慮好 NTT/FFT 的變換長度。
記得清空！！！記得清空！！！記得清空！！！

代碼：

void poly_inv(const int *a, int *res, int n) {
    static int ta[N], tr[N], b[N];
    b[0] = qpow(a[0], MOD - 2);
    for (int k = 1; k >> 1 < n; k <<= 1) {
        clr(ta, k << 1), clr(tr, k << 1);
        copy(a, a + min(n, k), ta), copy(b, b + min(n, k), tr);
        NTT(ta, k << 1), NTT(tr, k << 1);
        for (int i = 0; i < k << 1; ++i)
            b[i] = (2 - 1ll * ta[i] * tr[i] % MOD + MOD) % MOD * tr[i] % MOD;
        NTT(b, k << 1, -1);
        clr(b + k, k);
    }
    copy(b, b + n, res);
}

多項式開根

Luogu P5205 【模板】多項式開根 | Luogu P5277 【模板】多項式開根（加強版）

題意：給定一個 \(n-1\) 次多項式 \(f(x)\)，求出 \(g(x)\) 使得 \(g(x)^2\equiv f(x)\pmod{x^n}\)，系數對 \(998244353\) 取模。非加強版保證 \(a_0=1\)，加強版則只保證 \(a_0\) 是模 \(998244353\) 意義下的二次剩余。\(1\leq n\leq 10^5\)。

多項式開根同樣可以平方法做，讀者可以自行推導。這里直接給出牛頓迭代做法。

構造函數 \(h(g(x))=g(x)^2-f(x)\)，代入公式計算：

\[\begin{align*} g(x)& \equiv g_0(x)-\frac{g_0(x)^2-f(x)}{2g_0(x)} \\ & \equiv \frac{g_0(x)^2+f(x)}{2g_0(x)}\pmod{x^n} \end{align*} \]

對分母求逆即可計算上式。對于初值，若保證 \(a_0=1\)，直接取 \(g_0(x)\equiv 1\pmod{x^1}\) 即可，否則使用 Cipolla 求出模意義下開根的最小解即可。時間復雜度 \(O(n\log{n})\)。

加強版代碼：

void poly_sqrt(const int *a, int *res, int n) {
    int k;
    static tmpa[N], tmpr[N], inv[N];
    res[0] = cipolla(a[0]);
    chk_min(res[0], MOD - res[0]);
    for (k = 1; k >> 1 < n; k <<= 1) {
        clr(tmpa, k << 1), clr(tmpr, k << 1), clr(inv, k << 1);
        for (int i = 0; i < k << 1; ++i) tmpa[i] = tmpr[i] = inv[i] = 0;
        copy(a, a + k, tmpa), copy(res, res + k, tmpr);
        poly_inv(tmpr, inv, k), NTT(inv, k << 1);
        NTT(tmpa, k << 1);
        for (int i = 0; i < k << 1; ++i) res[i] = 1ll * tmpa[i] * inv[i] % MOD;
        NTT(res, k << 1, -1);
        for (int i = 0; i < k; ++i) res[i] = 1ll * (res[i] + tmpr[i]) % MOD * iv2 % MOD;
        for (int i = k; i < k << 1; ++i) res[i] = 0;
    }
    for (int i = n; i < k; ++i) res[i] = 0;
}

多項式除法

Luogu P4512 【模板】多項式除法

題意：給定一個 \(n\) 次多項式 \(f(x)\) 和一個 \(m\) 次多項式 \(g(x)\)，求出一個 \(n-m\) 次多項式 \(q(x)\) 和一個 \(m-1\) 次多項式 \(r(x)\)，使得 \(f(x)\equiv g(x)q(x)+r(x)\pmod{x^n}\)。\(1\leq m<n\leq 10^5\)。

這題的做法還是比較牛的。

若保證能整除，直接求逆即可，因此考慮去除 \(f(x)\) 的余多項式。前面的題中，我們經常用 \(\bmod{x^k}\) 來消去無用的項，這種做法可以消掉高次項。但余多項式顯然在低次項的位置，考慮反轉系數。令 \(f_{R}(x)=x^nf(\frac{1}{x})\)，我們嘗試推一下式子：

\[\begin{align*} f(x) & \equiv g(x)q(x)+r(x) & \pmod{x^n} \\ f(\frac{1}{x}) & \equiv g(\frac{1}{x})q(\frac{1}{x})+r(\frac{1}{x})& \pmod{x^n} \\ x^nf(\frac{1}{x})& \equiv x^mg(\frac{1}{x})x^{n-m}q(\frac{1}{x})+x^nr(\frac{1}{x})& \pmod{x^n} \\ f_R(x) & \equiv g_R(x)q_R(x)+x^{n-m+1}r_R(\frac{1}{x}) & \pmod{x^n} \\ q_R(x) & \equiv \frac{f_R(x)}{g_R(x)} & \pmod{x^{n-m+1}} \end{align*} \]

那么我們就可以多項式求逆求出 \(q_R(x)\)，反轉系數后得到 \(q(x)\)，最后 \(r(x)\equiv f(x)-g(x)q(x)\) 即可。時間復雜度 \(O(n\log{n})\)。

實現細節還是比較多的。代碼：

void poly_div(const int *a, const int *b, int *q, int *r, int n, int m) {
    static int ar[N], br[N], inv[N], tq[N];
    reverse_copy(a, a + n, ar), reverse_copy(b, b + m, br);
    for (int i = n - m + 1; i < m; ++i) br[i] = 0;
    poly_inv(br, inv, n - m + 1);
    int len = calc_len(n * 2 - m);
    NTT(ar, len), NTT(inv, len);
    for (int i = 0; i < len; ++i) ar[i] = 1ll * ar[i] * inv[i] % MOD;
    NTT(ar, len, -1);
    reverse_copy(ar, ar + n - m + 1, q), reverse_copy(ar, ar + n - m + 1, tq);
    clr(ar, len), clr(br, len), clr(inv, len);
    len = calc_len(n);
    copy(b, b + m, br);
    NTT(br, len), NTT(tq, len);
    for (int i = 0; i < len; ++i) r[i] = 1ll * br[i] * tq[i] % MOD;
    NTT(r, len, -1);
    for (int i = 0; i < m - 1; ++i) r[i] = (a[i] - r[i] + MOD) % MOD;
    clr(br, len), clr(tq, len);
}

多項式對數函數

Luogu P4725 【模板】多項式對數函數（多項式 ln）

題意：給定一個 \(n-1\) 次多項式 \(f(x)\)，求出一個多項式 \(g(x)\) 使得 \(g(x)\equiv \ln(f(x))\pmod{x^n}\)，系數對 \(998244353\) 取模。保證 \(a_0=1\)，\(1\leq n\leq 10^5\)。

\(\ln\) 這個函數比較棘手，但注意到其導數具有很好的性質，兩邊同時求導再積分就做完了：

\[g(x)\equiv\int g'(x)\textw0obha2h00x\equiv\int\frac{f'(x)}{f(x)}\textw0obha2h00x\pmod{x^n} \]

讀者可能會想到求導/積分會帶來常數項的損失，不過題目欽定了 \(a_0=1\)，所以我們無需擔心這一點。

那這又導出了另一個問題：如果不保證 \(a_0=1\) 呢？我們聲稱

\(f(x)\) 的對數多項式存在當且僅當 \(a_0=1\)。

考慮用麥克勞林級數對 \(\ln(f(x))\) 展開：

\[\begin{align*} \ln(1-x) & =-\sum_{n=1}^{+\infty}\frac{x^n}{n} \\ \ln(f(x)) & =-\sum_{n=1}^{+\infty}\frac{(1-f(x))^n}{n} \end{align*} \]

考察其常數項，即

\[[x^0]\ln(f(x))=-\sum_{n=1}^{+\infty}\frac{(1-a_0)^n}{n} \]

顯然該級數收斂當且僅當 \(a_0=1\)。得證。

代碼：

void poly_diff(const int *a, int *res, int n) {
    for (int i = 1; i < n; ++i) res[i - 1] = 1ll * a[i] * i % MOD;
    res[n - 1] = 0;
}
void poly_intg(const int *a, int *res, int n) {
    for (int i = 1; i < n; ++i) res[i] = 1ll * a[i - 1] * qpow(i, MOD - 2) % MOD;
    res[0] = 0;
}
void poly_ln(const int *a, int *res, int n) {
    static int df[N], inv[N];
    poly_inv(a, inv, n), poly_diff(a, df, n);
    int len = calc_len(n << 1);
    NTT(inv, len), NTT(df, len);
    for (int i = 0; i < len; ++i) df[i] = 1ll * df[i] * inv[i] % MOD;
    NTT(df, len, -1), poly_intg(df, res, n);
    clr(inv, len), clr(df, len);
}

多項式指數函數

Luogu P4726 【模板】多項式指數函數（多項式 exp）

題意：給定一個 \(n-1\) 次多項式 \(f(x)\)，求出一個多項式 \(g(x)\) 使得 \(g(x)\equiv e^{f(x)}\pmod{x^n}\)，系數對 \(998244353\) 取模。保證 \(a_0=0\)，\(1\leq n\leq 10^5\)。

\(\exp\) 函數怎么求導/積分都是本身，怎么辦？注意到其與 \(\ln\) 互為反函數，對等式兩邊同時取 \(\ln\)，得到

\[\ln(g(x))\equiv f(x)\Leftrightarrow \ln(g(x))-f(x)\equiv 0\pmod{x^n} \]

這個形式直接上牛頓迭代：

\[\begin{align*} g(x) & \equiv g_0(x)-\frac{\ln(g_0(x))-f(x)}{\frac{1}{g_0(x)}} \\ & \equiv g_0(x)(1-\ln(g_0(x))+f(x))\pmod{x^n} \end{align*} \]

倍增 \(O(n\log{n})\) 做即可。

和多項式對數函數類似，考慮 \(\exp(f(x))\) 的麥克勞林級數展開即可證明

\(f(x)\) 的指數多項式存在當且僅當 \(a_0=0\)。

代碼：

void poly_exp(const int *a, int *res, int n) {
    static int ta[N], tr[N], b[N], lnr[N];
    b[0] = 1;
    for (int k = 1; k >> 1 < n; k <<= 1) {
        clr(ta, k << 1), clr(tr, k << 1);
        copy(a, a + min(n, k), ta), copy(b, b + min(n, k), tr);
        poly_ln(tr, lnr, k);
        for (int i = 0; i < k << 1; ++i) ta[i] = (ta[i] - lnr[i] + MOD) % MOD;
        ta[0] = (ta[0] + 1) % MOD;
        NTT(ta, k << 1), NTT(tr, k << 1);
        for (int i = 0; i < k << 1; ++i) b[i] = 1ll * ta[i] * tr[i] % MOD;
        NTT(b, k << 1, -1);
        clr(b + k, k);
    }
    copy(b, b + n, res);
}

多項式冪函數

Luogu P5245 【模板】多項式快速冪 | Luogu P5273 【模板】多項式冪函數（加強版）

題意：給定一個 \(n-1\) 次多項式 \(f(x)\) 和 \(k\)，求出一個多項式 \(g(x)\) 使得 \(g(x)\equiv f(x)^k\pmod{x^n}\)，系數對 \(998244353\) 取模。\(1\leq n\leq 10^5\)，\(0\leq k\leq 10^{10^5}\)。非加強版保證 \(a_0=1\)，加強版則不做保證。

看到冪函數，一個思路就是兩邊同時取 \(\ln\) 再取 \(\exp\)，得到

\[g(x)\equiv e^{k\ln(f(x))}\pmod{x^n} \]

不過 \(k\) 很大，如何處理？考慮右側式子的麥克勞林級數展開：

\[\begin{align*} e^{k\ln(f(x))}&=-\sum_{n=0}^{+\infty}\frac{(k\ln(f(x)))^n}{n!} \\ &=-\sum_{n=0}^{+\infty}\frac{k^n\ln(f(x))^n}{n!} \end{align*} \]

因此我們可以令 \(k\leftarrow k\bmod{998244353}\)，而不影響答案正確性。

我們需要求 \(\ln(f(x))\)，而這要求 \(a_0=1\)，因此上述做法只能通過非加強版。

非加強版的代碼：

void poly_pow(int *a, int *res, int n, int b) {
    static int lna[N];
    poly_ln(a, lna, n);
    for (int i = 0; i < n; ++i) lna[i] = 1ll * lna[i] * b % MOD;
    poly_exp(lna, res, n);
    clr(lna, n);
}

cin >> n >> str + 1;
for (int i = 1; str[i]; ++i) k = (1ll * k * 10 + (str[i] ^ 48)) % MOD;
for (int i = 0; i < n; ++i) cin >> a[i];
poly_pow(a, ans, n, k);
for (int i = 0; i < n; ++i) cout << ans[i] << " \n"[i == n - 1];

既然上述做法要求 \(a_0=1\)，考慮對 \(f(x)\) 的每一項系數除掉 \(a_0\)，即

\[f(x)^k\equiv \left(\frac{f(x)}{a_0}\right)^ka_0^k\pmod{x^n} \]

這就做完了嗎？并沒有。加強版還不保證 \(a_0\neq 0\)。那考慮給 \(f(x)\) 降次再升次，即找到最小的 \(t\) 使得 \(a_t\neq 0\)，然后轉化成

\[f(x)^k\equiv \left(\frac{f(x)}{x^t}\right)^kx^{tk}\pmod{x^n} \]

也就是先把系數整體左移 \(t\) 個位置，沿用前面除掉常數項系數的做法，最后再整體右移 \(tk\) 個位置，這題就做完啦！

但還有很多細節。首先計算 \(a_0^k\) 需要利用歐拉定理，即

\[a_0^k\equiv a_0^{k\bmod{\varphi(998244353)}}\equiv a_0^{k\bmod{998244352}}\pmod{998244353} \]

需要和多項式的冪次做區分，即要存兩個 \(k\)，其中 \(k_1=k\bmod{998244353}\)，作為多項式的冪次，\(k_2=k\bmod{998244352}\)，作為計算 \(a_0^k\) 的冪次。

其次，我們需要特判偏移量 \(tk\geq n\) 的情況。由于 \(n\leq 10^5\)，所以直接把 \(k\) 在數位上截斷低 \(6\) 位判斷即可，即令 \(k_3=k\bmod{10^6}\)，判斷此時是否有 \(k_3\geq n\)，然后再去判斷是否有 \(tk_1\geq n\)。

代碼：

void poly_pow(int *a, int *res, int n, int b1, int b2) {
    int t = 0;
    static int ta[N], tr[N], lna[N];
    for (int i = 0; i < n && !a[i]; ++i, ++t);
    if (1ll * t * b1 >= n) { clr(res, n); return; }
    int a0 = a[t], iv = qpow(a0, MOD - 2);
    for (int i = 0; i < n; ++i) ta[i] = 1ll * a[i + t] * iv % MOD;
    poly_ln(ta, lna, n);
    for (int i = 0; i < n; ++i) lna[i] = 1ll * lna[i] * b1 % MOD;
    poly_exp(lna, tr, n);
    t *= b1;
    int pw = qpow(a0, b2);
    clr(res, t);
    for (int i = t; i < n; ++i) res[i] = 1ll * tr[i - t] * pw % MOD;
}

ios::sync_with_stdio(false), cin.tie(nullptr);
cin >> n >> str + 1;
for (int i = 1; str[i]; ++i) {
    k1 = (1ll * k1 * 10 + (str[i] ^ 48)) % MOD;
    k2 = (1ll * k2 * 10 + (str[i] ^ 48)) % (MOD - 1);
}
for (int i = 1; str[i] && i <= 6; ++i) k3 = k3 * 10 + (str[i] ^ 48);
for (int i = 0; i < n; ++i) cin >> a[i];
if (!a[0] && k3 >= n) {
    for (int i = 0; i < n; ++i) cout << 0 << " \n"[i == n - 1];
    return 0;
}
poly_pow(a, ans, n, k1, k2);
for (int i = 0; i < n; ++i) cout << ans[i] << " \n"[i == n - 1];
return 0;

多項式多點求值

Luogu P5050 【模板】多項式多點求值

題意：給定一個 \(n\) 次多項式 \(f(x)\) 和 \(m\) 個整數 \(p_i\)，\(\forall i\in[1,m]\) 求出 \(f(p_i)\)，答案對 \(998244353\) 取模。\(1\leq n,m\leq 6.4\times 10^4\)。

挺有趣的科技啊。下面默認 \(n,m\) 同階。

方法一：多項式取模

構造多項式 \(g_i(x)=x-p_i\)，設 \(f(x)\bmod{g_i(x)}=r_i\)，可以發現此時恰好有 \(r_i=f(p_i)\)，因為考慮將 \(f(x)\) 對 \(g_i(x)\) 進行多項式除法，得到 \(f(x)=g_i(x)q_i(x)+r_i\)，代入 \(x=p_i\)，此時 \(g_i(x)=0\)，所以 \(r_i=f(p_i)\)。容易將其進一步推廣：

令 \(g_{l,r}(x)=\prod_{i=l}^r(x-p_i)\)，\(f(x)\bmod{g_{l,r}(x)}=r_{l,r}(x)\)，則 \(\forall i\in[l,r],f(p_i)=r_{l,r}(p_i)\)。

于是我們得到了一個較為顯然的做法。首先容易 \(O(n\log^2{n})\) 分治 NTT 計算出 \(g_{l,r}(x)\)，然后，我們依然分治解決問題，設當前遞歸到 \([l,r]\)，并得到了多項式 \(f_{l,r}(x)\)，我們只需將 \(f_{l,r}(x)\) 分別對 \(g_{l,mid}(x)\) 和 \(g_{mid+1,r}(x)\) 取模，得到 \(f_{l,mid}(x)\) 和 \(f_{mid+1,r}(x)\)，繼續遞歸求解即可。當遞歸到邊界 \([l,l]\) 時，由上述推論，顯然此時 \(f_{l,l}(x)\) 僅有常數項且它就是 \(f(p_l)\) 的值。初始時只需令 \(f_{1,m}(x)=f(x)\bmod{g_{1,m}(x)}\)。取模保證了遞歸到 \([l,r]\) 時的多項式 \(f_{l,r}(x)\) 是一個 \(r-l\) 次多項式，復雜度得到了保證。由主定理，這部分復雜度還是 \(O(n\log^2{n})\) 的。

由于多項式取模的大常數，通常上述做法需要經過刻意卡常才能通過本題。

方法二：轉置原理

咕咕咕。

代碼：

void poly_mult(const int *a, const int *b, int *res, int n, int m) {
    int len = calc_len(n);
    static int ta[N], tb[N];
    copy(a, a + n, ta), reverse_copy(b, b + m, tb);
    NTT(ta, len), NTT(tb, len);
    for (int i = 0; i < len; ++i) ta[i] = 1ll * ta[i] * tb[i] % MOD;
    NTT(ta, len, -1);
    copy(ta + m - 1, ta + n, res);
    clr(ta, len), clr(tb, len);
}

namespace PolyEval {
#define ls(p) (p << 1)
#define rs(p) (p << 1 | 1)
    int inv[N];
    int len[N << 2], *g[N << 2], *q[N << 2];
    void eval_build(int p, int l, int r, const int *a) {
        if (l == r) {
            len[p] = 1;
            g[p] = new int[2], q[p] = new int[2];
            g[p][0] = 1, g[p][1] = (-a[l] + MOD) % MOD;
            return;
        }
        int mid = (l + r) >> 1;
        eval_build(ls(p), l, mid, a), eval_build(rs(p), mid + 1, r, a);
        len[p] = len[ls(p)] + len[rs(p)];
        g[p] = new int[len[p] + 1], q[p] = new int[len[p] + 1];
        poly_mul(g[ls(p)], g[rs(p)], g[p], len[ls(p)] + 1, len[rs(p)] + 1);
    }
    void eval_solve(int p, int l, int r, int *ans) {
        if (l == r) { ans[l] = q[p][0]; return; }
        int mid = (l + r) >> 1;
        poly_mult(q[p], g[rs(p)], q[ls(p)], r - l + 1, len[rs(p)] + 1);
        eval_solve(ls(p), l, mid, ans);
        poly_mult(q[p], g[ls(p)], q[rs(p)], r - l + 1, len[ls(p)] + 1);
        eval_solve(rs(p), mid + 1, r, ans);
    }
    void solve(const int *a, const int *b, int n, int m, int *ans) {
        static int tmp[N];
        eval_build(1, 1, m, b);
        poly_inv(g[1], inv, m + 1), reverse(inv, inv + m + 1);
        poly_mul(a, inv, tmp, n, m + 1), copy(tmp + n, tmp + n + m, q[1]);
        eval_solve(1, 1, m, ans);
        for (int i = 1; i <= m; ++i) ans[i] = (1ll * ans[i] * b[i] % MOD + a[0]) % MOD;
    }
#undef ls
#undef rs
}

多項式快速插值

Luogu P5158 【模板】多項式快速插值

題意：給出 \(n\) 個點 \((x_i,y_i)\)，求出一個 \(n-1\) 次多項式 \(f(x)\) 使得 \(\forall i\in[1,n],f(x_i)\equiv y_i\pmod{998244353}\)，系數對 \(998244353\) 取模。\(1\leq n\leq 10^5\)。

回顧拉格朗日插值公式：

\[f(x)=\sum_{i=1}^{n}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \]

不妨將分母提出，得到

\[f(x)=\sum_{i=1}^{n}\frac{y_i}{\prod_{j\neq i}{(x_i-x_j)}}\prod_{j\neq i}(x-x_j) \]

設 \(g(x)=\prod_{i=1}^n(x-x_i)\)，那么

\[\prod_{j\neq i}{(x_i-x_j)}=\frac{g(x_i)}{x_i-x_i} \]

而函數 \(\frac{g(x)}{x-x_i}\) 是連續可導的，考慮取極限并運用洛必達法則：

\[\frac{g(x_i)}{x_i-x_i}=\lim_{x\rightarrow x_i}\frac{g(x)}{x-x_i}=g'(x_i) \]

于是插值公式化為

\[f(x)=\sum_{i=1}^{n}\frac{y_i}{g'(x_i)}\prod_{j\neq i}(x-x_j) \]

這是一個比較典的可以分治計算的形式，具體來說，令

\[g_{l,r}(x)=\prod_{i=l}^r(x-x_i) \\ h_{l,r}(x)=\sum_{i=l}^{r}\frac{y_i}{g'(x_i)}\prod_{\substack{j\neq i\\l\leq j\leq r}}(x-x_j) \]

即 \(h_{l,r}(x)\) 表示 \((x_l,y_l),\cdots,(x_r,y_r)\) 插值得到的多項式。我們考慮缺口位于左半區間還是右半區間，推一下式子：

\[\begin{align*} h_{l,r}(x) & =\sum_{i=l}^{r}\frac{y_i}{g'(x_i)}\prod_{\substack{j\neq i\\l\leq j\leq r}}(x-x_j) \\ & = \left(\prod_{i=mid+1}^r(x-x_i)\right)\sum_{i=l}^{mid}\frac{y_i}{g'(x_i)}\prod_{\substack{j\neq i\\l\leq j\leq mid}}(x-x_j)+\left(\prod_{i=l}^{mid}(x-x_i)\right)\sum_{i=mid+1}^{r}\frac{y_i}{g'(x_i)}\prod_{\substack{j\neq i\\mid+1\leq j\leq r}}(x-x_j) \\ & = g_{mid+1,r}(x)h_{l,mid}(x)+g_{l,mid}(x)h_{mid+1,r}(x) \end{align*} \]

把 \(g'(x)\) 對 \(x_i(i\in[1,n])\) 進行多點求值，然后分治 \(O(n\log^2{n})\) 解決即可。

代碼：

namespace PolyInter {
#define ls(p) (p << 1)
#define rs(p) (p << 1 | 1)
    int val[N];
    int len[N << 2], *g[N << 2], *h[N << 2];
    void inter_build(int p, int l, int r, const int *a) {
        if (l == r) {
            len[p] = 1;
            g[p] = new int[2], h[p] = new int[2];
            g[p][0] = -a[l] + MOD, g[p][1] = 1;
            return;
        }
        int mid = (l + r) >> 1;
        inter_build(ls(p), l, mid, a), inter_build(rs(p), mid + 1, r, a);
        len[p] = len[ls(p)] + len[rs(p)];
        g[p] = new int[len[p] + 1], h[p] = new int[len[p] + 1];
        poly_mul(g[ls(p)], g[rs(p)], g[p], len[ls(p)] + 1, len[rs(p)] + 1);
    }
    void inter_solve(int p, int l, int r, const int *y) {
        if (l == r) { h[p][0] = 1ll * y[l] * qpow(val[l], MOD - 2) % MOD; return; }
        int mid = (l + r) >> 1;
        inter_solve(ls(p), l, mid, y), inter_solve(rs(p), mid + 1, r, y);
        static int tmp[N];
        poly_mul(h[ls(p)], g[rs(p)], tmp, len[ls(p)], len[rs(p)] + 1);
        poly_mul(h[rs(p)], g[ls(p)], h[p], len[rs(p)], len[ls(p)] + 1);
        for (int i = 0; i < len[p]; ++i) h[p][i] = (h[p][i] + tmp[i]) % MOD;
        clr(tmp, len[p]);
    }
    void solve(const int *x, const int *y, int n) {
        static int df[N];
        inter_build(1, 1, n, x);
        poly_diff(g[1], df, n + 1);
        PolyEval::solve(df, x, n + 1, n, val);
        inter_solve(1, 1, n, y);
        clr(df, n);
    }
#undef ls
#undef rs
}

例題

Luogu P5395 第二類斯特林數·行

題意：給定 \(n\)，對于所有整數 \(i\in[0,n]\)，求出 \({n\brace i}\)，答案對 \(167772161=5\times 2^{25}+1\) 取模。\(1\leq n\leq 2\times 10^5\)。

這里不加推導地給出第二類斯特林數的通項公式，其由二項式反演容易得到：

\[{n\brace m}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!} \]

遇到類似 \(i\) 和 \(m-i\) 這類和為定值的形式，考慮構造成加法卷積的形式。對這題而言，

\[\begin{align*} {n\brace m} & =\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!} \\ & =\sum_{i=0}^m\frac{i^n}{i!}\times\frac{(-1)^{m-i}}{(m-i)!} \end{align*} \]

令 \(f_i=\frac{i^n}{i!},g_i=\frac{(-1)^i}{i!}\)，則 \({n\brace m}=\sum_{i=0}^mf_ig_{m-i}\)。這就是很顯然的卷積形式了，構造多項式 \(f(x)\) 滿足 \([x^i]f(x)=f_i\)，和 \(g(x)\) 滿足 \([x^i]g(x)=g_i\)，NTT 計算它們的卷積即可。時間復雜度 \(O(n\log{n})\)。

AtCoder ABC196F Substring 2

題意：給定 \(0/1\) 串 \(s,t\)，求最少修改多少次 \(t\) 可以使其成為 \(s\) 的子串。\(1\leq |t|\leq|s|\leq 10^6\)。

令 \(n=|s|,m=|t|\)，枚舉把 \(t\) 改成子串 \(s[k,k+m-1]\)，此時的修改次數為

\[\sum_{i=0}^{m-1}t_i\oplus s_{k+i} \]

操作數是 \(0/1\)，此時我們有經典結論 \(x\oplus y=(x-y)^2=x-2xy+y\)，代回上面的式子

\[\sum_{i=0}^{m-1}t_i+\sum_{i=0}^{m-1}s_{k+i}-2\sum_{i=0}^{m-1}t_is_{k+i} \]

\(\sum_{i=0}^{m-1}t_i\) 顯然是一個定值，\(\sum_{i=0}^{m-1}s_{k+i}\) 可以前綴和處理，于是重點在于計算 \(\sum_{i=0}^{m-1}t_is_{k+i}\)。式子中出現了乘法，但不是我們所熟知的和為定值的形式，而是差為定值，對于這種形式，有經典套路：將其中一個多項式反轉。對這題而言，令 \(s'_i=s_{n-1-i}\)，則 \(\sum_{i=0}^{m-1}t_is_{k+i}=\sum_{i=0}^{m-1}t_is_{n-1-k-i}'\)，此時和為定值 \(n-1-k\)，就可以 NTT 直接做了。

posted @ 2025-03-01 21:53 P2441M 閱讀(89) 評論(0) 收藏舉報

刷新頁面返回頂部

P2441M

算法學習筆記：多項式

多項式基礎

多項式的表示方法

離散卷積

拉格朗日插值

橫坐標是連續整數的快速單點插值

快速傅里葉變換

離散傅里葉變換

分治實現 FFT

倍增實現 FFT/蝴蝶變換

多項式求值和線性代數

快速傅里葉逆變換

DFT 和 FFT

快速數論變換

FFT/NTT 的應用

多項式乘法

分治 FFT/NTT

半在線卷積

單項式乘積

多項式牛頓迭代

泰勒級數和麥克勞林級數

一般的牛頓迭代

多項式牛頓迭代

多項式乘法逆

方法一：平方法

方法二：牛頓迭代

倍增法的細節

多項式開根

多項式除法

多項式對數函數

多項式指數函數

多項式冪函數

多項式多點求值

方法一：多項式取模

方法二：轉置原理

多項式快速插值

例題

Luogu P5395 第二類斯特林數·行

AtCoder ABC196F Substring 2

公告