多项式 学习笔记

· · 个人记录

和计算几何一样是东西很多很杂而且都是板子的一个工具性专题。

东西太多了，只学了 zaky 板子里包含的东西。

题单可以看 这个 和 卡老师的遗产。

拉格朗日插值

用途：\mathcal O(n^2) 点值表示转多项式。

题意：给 n 个点 (x_i,y_i)，求一个唯一的度数不大于 n-1 的多项式过这 n 个点。

朴素的做法是设出这 n 个系数 a_{1\cdots n}，然后待定系数法用高斯消元求解。

复杂度是 \mathcal O(n^3) 的，不够优秀。

拉格朗日插值的导出有两种形式，一种在模意义下用 CRT 构造，一种纯代数，这里只讲后者。

我们对每个点 P_i(x_i,y_i)，构造其在 x 轴上的投影 P_i^\prime=(x_i,0)。

假如我们构造了 n 个多项式 f_{1\cdots n}，使得 f_i 过 P_i 和 P_j^\prime(i\ne j)，那么由答案的唯一性，\sum f_i 就是答案。

这个构造是不难的。由于 n-1 个零点已经给出，我们设：

f_i=a\prod_{i\ne j}(x-x_j)

而 f_i(x_i)=y_i，可以得到：

a=\dfrac{y_i}{\prod_{i\ne j}(x_i-x_j)}

因此：

f_i=y_i\cdot \dfrac{\prod_{i\ne j}(x-x_j)}{\prod_{i\ne j}(x_i-x_j)}=y_i\cdot \prod_{i\ne j}\dfrac{x-x_j}{x_i-x_j}

答案即为：

F=\sum_{i=1}^ny_i\cdot \prod_{j\ne i}\dfrac{x-x_j}{x_i-x_j}

这个式子可以方便地代入一个值 k，求出 F(k)。

复杂度 \mathcal O(n^2)，注意所有分母要乘起来之后取逆元，否则会变成 \mathcal O(n^2\log P)。

fo(i,1,n) {
    ll cur = y[i], inv = 1;
    fo(j,1,n) if(i != j) {
        cur = cur * (k - x[j]) % P;
        inv = inv * (x[i] - x[j]) % P;
    } ans = (ans + cur * power(inv, P - 2) % P) % P;
} ans = (ans % P + P) % P;

注意如果处理不当，拉格朗日插值得到的答案很可能是负数，一定要 ans = (ans % P + P) % P;！！！

考虑一些特殊的应用：

这时式子变成：

\sum_{i=1}^n y_i\prod_{j\ne i}\dfrac{k-j}{i-j}

我们预处理 \text{Pre}_x=\prod_{i=1}^x(k-i), \text{Suf}_x=\prod_{i=x}^n(k-i) 和 n 以内所有数的阶乘，原式可以化成：

\sum_{i=1}^ny_i\cdot \dfrac{\text{Pre}_{i-1}\times \text{Suf}_{i+1}}{(i-1)!\times (n-i)!\times (-1)^{n-i}}

这样就可以 \mathcal O(n) 求了。

容易发现这个拉格朗日插值法没有办法支持动态插一个点。

我们整理一下式子：

F=\sum_{i=1}^ny_i\cdot \prod_{j\ne i}\dfrac{x-x_j}{x_i-x_j}=\prod_{i=1}^n(x-x_i)\sum_{i=1}^n \dfrac{y_i}{x-x_i}\prod_{j\ne i}\dfrac{1}{x_i-x_j}

其实每个部分都不难维护了。设 g=\prod_{i=1}^n(x-x_i),h_i=\dfrac{y_i}{\prod_{j\ne i}(x_i-x_j)}。

每次加入一个点可以 \mathcal O(n) 维护 g,h，得到：

F=g\sum_{i=1}^n\dfrac{h_i}{(k-x_i)}

这被称为重心拉格朗日插值法。

其实不是很难，我们计算 c_i=y_i \prod_{j\ne i}\frac{1}{x_i-x_j}，那么：

F=\sum_{i=1}^ny_i\cdot \prod_{j\ne i}\dfrac{x-x_j}{x_i-x_j}=\sum_{i=1}^nc_i\times \dfrac{\prod_{j=1}^n (x-x_j)}{x-x_i}

暴力卷积预处理出 \prod_{j=1}^n (x-x_j)，每次除以一个 x-x_i 再乘系数即可。

复杂度 \mathcal O(n^2)。

一般解题过程中，如果可以证明答案是一个 n 次多项式，那么我们就可以找好算的 n+1 个点然后用拉格朗日插值法求出答案。

显然这是一个 k+1 次多项式，我们需要插入 k+2 个点。

求这 k+2 个点可以线性筛，也可以 \mathcal O(k\log P) 暴力。

求完了就是坐标连续插值的板子了。

快速傅里叶变换

用途：\mathcal O(n\log n) 计算两多项式相乘（卷积）。

离散傅里叶变换

所谓离散傅里叶变换（Discrete Fourier Transform, DFT），在多项式意义下就是把系数表示的多项式转为点值表示。

反之，IDFT 就是把点值表示的多项式转为系数表示。

为什么计算两多项式相乘要转成点值表示？因为多项式相乘可以直接转化为点值相乘。具体地：

F\times G=H\Rightarrow \forall x, F(x)G(x)=H(x)

十分漂亮。

由于点值表示的多项式对点的位置没有特殊要求，所以我们会选取若干特殊的 x 值来加速 DFT 和 IDFT 的过程。

我们取一个 n 是最小的比乘积度数大的 2 的整次幂，用所有的 \omega_n^k （n 次单位根）来做这些 x 值进行 DFT。

这个东西有什么好处呢？我们发现，如果把最终答案系数 a_{0\cdots n-1} 和计算出的点值 y_{0\cdots n-1} 看做向量，有：

\begin{bmatrix} 1 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ 1 & \omega_n^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1)^{n-1} \\ 1 & \omega_n^2 & (\omega_n^2)^2 & \cdots & (\omega_n^2)^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} &(\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\a_2\\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} y_0\\ y_1\\y_2\\ \vdots \\ y_{n-1} \end{bmatrix}

这意味着我们只需要求出左边那个矩阵（这被称为傅里叶矩阵 F_n）的逆矩阵，再左乘于答案上即可。

事实上，我们将 F_n 中的每一位取反得到 G_n，展开可以得到 G_nF_n=nI，因此 F_n^{-1}=\dfrac{G_n}{n}。

这意味着，我们直接代入 (\omega_n^0,\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-n+1}) 再进行一遍 DFT 所得到的向量正好是答案的 n 倍，除以 n 即可。

分治优化

怎么进行 DFT 的过程呢？我们考虑分治：

对于一个 n-1 次多项式 A(x)，我们把 A(x) 中的每一项按次数的奇偶性拆分成 A_1(x^2)+xA_2(x^2)。

具体地，设 A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1}，我们有：

A_1(x)=a_0+a_2x+\cdots +a_{n-2}x^{\frac{n}{2}-1} A_2(x)=a_1+a_3x+\cdots +a_{n-1}x^{\frac{n}{2}-1}

容易发现 A_1,A_2 都是 \dfrac{n}{2}-1 次多项式，分治下去求出所有 A_1(\omega_{\frac{n}{2}}^k),A_2(\omega_{\frac{n}{2}}^k) 的值。

然后，我们对 k<\dfrac{n}{2}，有：

A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_{\frac{n}{2}}^k) A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k})=A_1(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_2(\omega_{\frac{n}{2}}^k)

于是就可以 \mathcal O(n\log n) 递归做了。

常数优化与实现

但是原始写法常数太大，我们优化一下。

注意细节：

https://duck.ac/problem/1002 代码：

#include<bits/stdc++.h>
using namespace std;
#define fo(v,a,b) for(int v = a; v <= b; v++)
#define fr(v,a,b) for(int v = a; v >= b; v--)
#define cl(a,v) memset(a, v, sizeof(a))

const int N = 3000000;

typedef double db;
typedef complex<db> cp;

inline int read()
{
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}

int n, L, len1, len2;

void FFT(cp *a, cp *omg) {
    fo(i,0,n-1) {
        int t = 0;
        fo(j,0,L-1) if(i >> j & 1) t |= (1 << (L - 1 - j));
        if(t > i) swap(a[i], a[t]);
    }
    for(int l = 2; l <= n; l <<= 1) {
        int m = l / 2;
        for(cp *p = a; p != a + n; p += l)
            fo(i,0,m-1) {
                cp t = p[i + m] * omg[n / l * i];
                p[i + m] = p[i] - t, p[i] += t;
            }
    }
}

cp a[N], b[N], c[N], omg[N], inv[N];

void poly_multiply(unsigned *v1, int len1, unsigned *v2, int len2, unsigned *v3)
{
    fo(i,0,len1) a[i].real(v1[i]);
    fo(i,0,len2) b[i].real(v2[i]);

    n = 1; while(n <= len1 + len2) n <<= 1, L++;
    fo(i,0,n-1) {
        omg[i] = cp(cos(2.0 * M_PI * i / n), sin(2.0 * M_PI * i / n));
        inv[i] = conj(omg[i]);
    }

    FFT(a, omg), FFT(b, omg);
    fo(i,0,n-1) c[i] = a[i] * b[i];
    FFT(c, inv);
    fo(i,0,len1 + len2) v3[i] = (unsigned) floor(c[i].real() / n + 0.5);
}

FFT 的一个主要的运用是求离散卷积，即求：

H(i)=\sum_{j=0}^i F(j)G(i-j)

把式子分成左右两边分别计算即可。

和卷积与差卷积

有的时候要计算的卷积不是传统的 H(i)=\sum_j F(j)G(i-j)，而是反一反的形式：

H(i)=\sum_{j=i}^n F(j)G(j-i)

这被叫做差卷积。

这时可以尝试把 F 倒一倒，构造 F^\prime(i)=F(n-i)，那么：

H(i)=\sum_{j=i}^n F^\prime(n-j)G(j-i)

再把 H 倒一倒，构造 H^\prime(i)=H(n-i)，就有 H^\prime = F^\prime \times G 了。

类似地还有和卷积：

H(i)=\sum_{j=0}^{n-i} F(j)G(i+j)

这会我们有 H^\prime =F\times G^\prime，当然翻另一个也都是可以的。

快速数论变换

用途：带模数的卷积。

单位根与原根

NTT 是 FFT 在 \Z_p 环上的一个应用，解决了带模数的卷积问题。

具体地，我们给 FFT 在 \mathbb C 中要用到的单位根 \omega_n 在 \Z_p 中找个对应：原根。

比如，998244353=2^{23}\times 7 \times 17+1，其最小原根为 3。

这意味着 \omega_{n}=g^{\frac{P-1}{n}}，需要满足 n\mid (P - 1) 才能有可以用的原根。

而 2^{23}\mid (998244353-1)，可以满足绝大多数的要求。

实现时可以预处理好所有蝴蝶变换的值，减小常数。

现在开始，我们要尝试对多项式板子做一个封装了。

这个板子在文章的最后会展示，所以后面的东西都不会放代码。

分治 FFT / NTT

一个技巧。

一种形式是：我们要计算 n 个多项式的乘积，每个只有两项。

由于对于长度分别为 a,b 的两个多项式，卷积复杂度是 f(a+b) 的，所以可以分治把复杂度降到 \mathcal O(n\log^2 n)。

另一种形式是：我们有自己向自己的卷积形式的递推，如P4721 【模板】分治 FFT。

做法也很显然，就是 CDQ 分治。

void solve(int l, int r) {
    if(l == r) return ;
    int mid = (l + r) >> 1;
    solve(l, mid);

    A.reset(), A.resize(mid - l + 1);
    B.reset(), B.resize(r - l + 1);
    fo(i,l,mid) A[i - l] = f[i];
    fo(i,0,r-l) B[i] = g[i];
    A = A * B;
    fo(i,mid+1,r) f[i] = (f[i] + A[i - l]) % P;

    solve(mid + 1, r);
}

任意模数 NTT

做法是不难的，过程就是取若干模数进行 NTT（一般三个），保证它们的乘积大于答案。

然后 CRT 合并，容易发现可以确定出唯一答案，这样边做边取模即可。

这个合并可以直接 __int128 搞，也可以先合并前两个再快速乘搞出答案。

多项式求逆

考虑分治。

由于 m\ge n 时，\bmod x^m 的逆元一定也是 \bmod x^n 的逆元，所以我们找到一个 \ge n 的 2 的整幂次处理其逆元。

如果我们已经求出 F(x)G_0(x)\equiv 1\pmod{x^{\frac{n}{2}}}，那么就是 x^{\frac{n}{2}}\mid F(x)G_0(x) - 1。

平方得到：x^n\mid (F(x)G_0(x) - 1)^2，做简单的代数推导可以得到：

F(x)(2\times G_0(x)-F(x)G_0^2(x))\equiv 1\pmod{x^n}

因此新的逆元满足：

G(x)=G_0(x)(2-F(x)G_0)

可以直接递归 / 递推计算。

复杂度的话，由于我们每次只需要保留 F 和 G_0 的前 \dfrac{n}{2} 项，所以复杂度满足：

T(n)=T\left(\dfrac{n}{2}\right)+\Theta(n\log n)

得到 T(n)=\Theta(n\log n)。

注意实现时可以少做几次 NTT 减小常数。

多项式指数函数 / 对数函数

泰勒级数与麦克劳林级数

我们先大致了解一下泰勒级数和麦克劳林级数，由于是大致了解所以不会写出什么严谨过程。

我们知道，对任意一个可导的函数 f(x)，它在 x=a 处的切线方程可以写作 y=f(a)+f^\prime(a)(x-a)。

类似地，我们找到一个最接近 f(a) 的二次切线，可以写出其方程为 y=f(a)+f^\prime(a)(x-a)+\dfrac{f^{\prime\prime}(a)}{2}(x-a)^2。

将其推广，我们找到了一个用 N 次多项式近似 f 的方法：

P_N(x)=\sum_{n=0}^N\dfrac{f^{(n)}(a)}{n!}(x-a)^n

这被称为泰勒多项式，实际使用时可以取前几项作为近似。

而如果把 N\to \infty，我们就得到了泰勒级数。而麦克劳林级数是取 a=0 的一个特例，形式是：

\sum_{n=0}^{\infty}\dfrac{f^{(n)}(0)}{n!}x^2

如果上面这个式子和 f(x) 相等，那就好了，但是并不总是这样。

判断泰勒级数在 x 处收敛于 f(x) 只需证明其余项 R_N(x)\to 0，其中余项的定义为：

R_N(x)=\dfrac{f^{(N+1)}(c)}{(N+1)!}(x-a)^{N+1}

其中 c 在 a,x 之间。

比如，我们可以给 e^x 做一个麦克劳林展开：

\sum_{n=0}^{N}\dfrac{e^0x^n}{n!}=\sum_{n=0}^{N}\dfrac{x^n}{n!}

那么

R_N(x)=\dfrac{e^c}{(N+1)!}x^{N+1}

由于 e^c\le \max(e^x,1)，所以可以看做一个常数，这样由于 \dfrac{x^n}{n!} 收敛于 0，我们证明了：

e^x=\sum_{n=0}^{N}\dfrac{x^n}{n!}=1+x+\dfrac{x^2}{2}+\dfrac{x^3}{6}+\cdots

对于 \ln 的处理就没有那么简单了，因为 \ln x 在 0 处没有定义。

我们考虑曲线救国，对 \ln (1-x) 进行麦克劳林展开。

归纳可得：

\ln^{(n)} (1-x)=(-1)^{n-1}\dfrac{(n-1)!}{(x-1)^n}

所以有：

\ln (1-x)=\sum_{i=1}^\infty \dfrac{(-1)^{i-1}(i-1)!}{(-1)^ii!}x^i=-\sum_{i=1}^\infty\dfrac{x^i}{i}

形式出奇地简单！

具体的证收敛过程不会，反正在 |x|<1 或 x=-1 时收敛。

在多项式上的推广

我们考虑把这个定义推广到多项式上。对于一个多项式 f(x)，我们定义：

\exp f(x)=e^{f(x)}=\sum_{i=0}^\infty \dfrac{f^i(x)}{i!} \ln (1-f(x))=-\sum_{i=1}^\infty\dfrac{f^i(x)}{i}

代入 g(x)=1-f(x)，可以得到多项式 \ln 的形式：

\ln g(x)=-\sum_{i=1}^{\infty}\dfrac{(1-g(x))^i}{i}

对数函数在模意义下收敛当且仅当 [x^0]f(x)=1，还是不会证。

牛顿迭代法

求 \exp 可以用牛顿迭代法，它解决的问题是：

给定多项式 g(x)，求 f(x) 使 x^n\mid g(f(x))。

可以分治。n=1 时直接取 f(x)=-[x^0]g(x)。

设我们已经找到 \bmod x^{\frac{n}{2}} 意义下的解 f_0(x)，将 g(f(x)) 在 f_0(x) 处泰勒展开：

\sum_{i=0}^n \dfrac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i\equiv 0\pmod{x^n}

由于 i\ge 2 时的项 \bmod x^n=0，所以我们只需要关注：

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

整理一下：

f(x)\equiv f_0(x)-\dfrac{g(f_0(x))}{g^\prime(f_0(x))}\pmod{x^n}

这个形式被叫作牛顿迭代法，作 \log n 次即可得解。

多项式 exp / ln 的求法

考虑求多项式 \ln，即解：

g(f(x))=\exp (B(x))-A(x)\equiv 0\pmod{x^n}

由于 [x^0]A(x)=1，所以 n=1 时直接 f(x)=0。

有递推式：

f(x)\equiv f_0(x)-\dfrac{\exp (B(x))-A(x)}{\exp(B(x))}\equiv f_0(x)-1-\dfrac{A(x)}{\exp(B(x))}\pmod{x^n}

好像走不通？

其实有一个更好用的性质：\ln 的导数。由链式法则：

\ln^{\prime} f(x)=\dfrac{f^\prime(x)}{f(x)}

因此：

\int \dfrac{f^\prime(x)}{f(x)}\text{d}x=\ln f(x)

可以 \mathcal O(n\log n) 计算。

再考虑求多项式 \exp，这时牛迭就能派上用场了：

g(f(x))=\ln f(x)-h(x)\equiv 0\pmod {x^n}

在保证 [x^0]h(x)=0 时易得 f_0(x)=1，递推式是：

f(x)\equiv f_0(x)-\dfrac{\ln f_0(x)-h(x)}{\frac{1}{f_0(x)}}\equiv f_0(x)(1-\ln f_0(x)+h(x))\pmod{x^n}

多项式开方

依然采用牛顿迭代法，令 g(f(x))=f^2(x)-h(x)，我们有递推式：

f(x)=f_0(x)-\dfrac{f_0^2(x)-h(x)}{2f_0(x)}=\dfrac{1}{2}\left(f_0(x)+\dfrac{h(x)}{f_0(x)}\right)

多项式快速幂

由于 (a+b)^P\equiv a^P+b^P\pmod{P}，所以我们可以得到：

f^P(x)\equiv \sum_{i=0}^na_i^Px^{iP}\equiv a_0^P\pmod{P}

这意味着 a_0=1 时我们可以直接求 f^{k\ \bmod\ P}(x)，这样我们有：

f^k(x)=\exp(k\cdot \ln f(x))

否则怎么办呢？我们可以考虑：

f^k(x)=\left(\dfrac{f(x)}{[x^0]f(x)}\right)^k\times ([x^0]f(x))^k

这适用于 [x^0]f(x)\ne 0 的情况。

板子

由于 \texttt{z}\color{red}\texttt{houkangyang} 常用的板子只有这些，那就只有这些吧。

#include<bits/stdc++.h>
using namespace std;
#define fo(v,a,b) for(int v = a; v <= b; v++)
#define fr(v,a,b) for(int v = a; v >= b; v--)
#define cl(a,v) memset(a, v, sizeof(a))

typedef long long ll;

const ll P = 998244353, gP = 3;
const int N = 600000;

ll power(ll a, ll x) {
    ll res = 1;
    while(x) {
        if(x & 1LL) res = res * a % P;
        a = a * a % P, x >>= 1;
    } return res;
}
inline ll inv(ll x) { return power(x, P - 2);  }

int Len, lg[N], Trans[2 * N], TransPos, Beg[30];
ll omgVal[N], invVal[N], numInv[N];

struct Poly
{
    vector<ll> v;

    // Basic Operations
    int siz() {  return (int)v.size();  } 
    ll &operator[](int x) {  return v[x];  }
    ll vv(int x) {  return x < 0 || x >= (int)v.size() ? 0 : v[x]; }
    void reset() {  v.clear();  }
    void resize(int x) {  v.resize(x);  }
    Poly(vector<ll> vec = {0}) { v = vec; }

    // Linear Operations
    friend Poly operator*(Poly A, ll c) {
        for(ll &x : A.v) x = x * c % P;
        return A;
    }
    friend Poly operator*(ll c, Poly A) {
        return A * c;
    }
    friend Poly operator+(Poly A, Poly B) {
        Poly res; res.resize(max(A.siz(), B.siz()));
        fo(i,0,res.siz()-1) res[i] = (A.vv(i) + B.vv(i)) % P;
        return res;
    }
    friend Poly operator-(Poly A, Poly B) {
        Poly res; res.resize(max(A.siz(), B.siz()));
        fo(i,0,res.siz()-1) res[i] = (A.vv(i) + P - B.vv(i)) % P;
        return res;
    }
    Poly& operator+=(Poly A) {  return (*this) = (*this) + A;  }
    Poly& operator-=(Poly A) {  return (*this) = (*this) - A;  }
    Poly& operator*=(Poly A) {  return (*this) = (*this) * A;  }

    // Convolution
    void NTT(ll *omg) {
        int n = v.size(), L = lg[n];
        if(n != (1 << L)) L++, v.resize(n = (1 << L));
        assert(n <= Len);
        fo(i,0,n-1) if(Trans[Beg[L] + i] > i) swap(v[i], v[Trans[Beg[L] + i]]);
        for(int l = 2; l <= n; l <<= 1) {
            int m = l / 2, rate = Len / l;
            for(vector<ll>::iterator p = v.begin(); p != v.end(); p += l)
                fo(i,0,m-1) {
                    ll cur = p[i + m] * omg[rate * i] % P;
                    p[i + m] = (p[i] + P - cur) % P, p[i] = (p[i] + cur) % P;
                }
        }
    }
    friend Poly operator*(Poly iA, Poly iB) {
        Poly A = iA, B = iB;
        int cur = A.siz() + B.siz() - 1, mx = 1;
        while(mx <= cur) mx <<= 1;
        A.resize(mx), B.resize(mx), A.NTT(omgVal), B.NTT(omgVal);
        ll iv = inv(mx); fo(i,0,mx-1) A.v[i] = A.v[i] * B.v[i] % P;
        A.NTT(invVal); fo(i,0,mx-1) A.v[i] = A.v[i] * iv % P;
        return A.resize(mx), A;
    }

    // Inversion
    Poly Inv() { // 0.2 s
        Poly F, G; G.resize(1), G[0] = inv(v[0]);
        for(int m = 2, pm = 4; m / 2 < (int)v.size(); m <<= 1, pm <<= 1) {
            F.reset(), F.resize(pm), G.resize(pm);
            fo(i,0,m-1) F[i] = this->vv(i);
            G.NTT(omgVal), F.NTT(omgVal);
            fo(i,0,pm-1) G[i] = G[i] * (P + 2 - G[i] * F[i] % P) % P;
            G.NTT(invVal); ll invm = inv(pm);
            fo(i,0,pm-1) G[i] = G[i] * invm % P;
            G.resize(m);
        } return G.resize(v.size()), G;
    }

    // Derivative & Integral
    Poly deriv() {
        Poly res; res.resize(v.size() - 1);
        fo(i,1,(int)v.size()-1) res[i - 1] = v[i] * i % P;
        return res;
    }
    Poly integ() {
        Poly res; res.resize(v.size() + 1);
        fo(i,0,(int)v.size()-1) res[i + 1] = v[i] * numInv[i + 1] % P;
        return res;
    }

    // exp & ln
    Poly Ln() { // 0.2 s
        assert(v[0] == 1);
        Poly res = ((this -> deriv()) * (this -> Inv())).integ();
        return res.resize(v.size()), res;
    }
    Poly Exp() { // 0.5 s
        assert(v[0] == 0); Poly res({1});
        for(int m = 2, pm = 4; m / 2 < (int)v.size(); m <<= 1, pm <<= 1) {
            Poly f; res.resize(m), f.resize(m);
            fo(i,0,m-1) f[i] = this->vv(i);
            res = res * (Poly({1}) - res.Ln() + f);
            res.resize(m);
        } return res.resize(v.size()), res;
    }

    // Square Root : 0.7 s
    Poly Sqrt() {
        assert(v[0] == 1); Poly res({1}); ll Inv2 = (P + 1) / 2;
        for(int m = 2, pm = 4; m / 2 < (int)v.size(); m <<= 1, pm <<= 1) {
            Poly f; res.resize(pm), f.resize(pm);
            fo(i,0,m-1) f[i] = this->vv(i);
            res = (res + f * res.Inv()) * Inv2;
            res.resize(m);
        } return res.resize(v.size()), res;
    }

    // Quick Pow : 0.7 s
    Poly Shift(int x) { // x^{n} -> x^{n+k}
        if(v.size() + x <= 0) return Poly();
        Poly res; res.resize(v.size() + x);
        fo(i,0,(int)v.size()-1) if(i + x >= 0) res[i + x] = v[i];
        return res;
    }
    Poly Power_1(ll k) {
        Poly res = (k * (this->Ln()));
        res = res.Exp();
        return res;
    }
    Poly Power_not0(ll k, ll k2) { // the real k mod P, k mod (P - 1)
        Poly res(v); ll val = inv(v[0]);
        res = res * val, res = res.Power_1(k);
        return res = res * power(v[0], k2), res;
    }
    Poly Power(ll k, ll k2, bool fl = true) { // fl is (if the real k < P)
        if(v[0]) return Power_not0(k, k2);
        Poly res({0}); if(!fl) return res.resize(v.size()), res;
        int c = -1; fo(i,1,(int)v.size()-1) if(v[i]) {  c = i; break;  }
        assert(c != -1);
        if(c * k >= (int)v.size())
            return res = Poly(), res.resize(v.size()), res;
        res = (this->Shift(-c)).Power_not0(k, k2), res = res.Shift(c * k);
        return res.resize(v.size()), res;
    }
} ;

void Pinit(int x) {
    Len = 1; while(Len <= x) Len <<= 1;
    fo(i,2,Len) lg[i] = lg[i >> 1] + 1;
    fo(j,0,lg[Len]) {
        Beg[j] = TransPos;
        fo(i,0,(1<<j)-1) {
            fo(now,0,j-1) if(i >> now & 1) Trans[TransPos] |= (1 << (j - 1 - now));
            TransPos++;
        }
    }
    omgVal[0] = invVal[0] = numInv[1] = 1;
    ll Sg = power(gP, (P - 1) / Len), Sinvg = inv(Sg);
    fo(i,1,Len-1) {
        omgVal[i] = omgVal[i - 1] * Sg % P;
        invVal[i] = invVal[i - 1] * Sinvg % P;
    }
    fo(i,2,Len) numInv[i] = numInv[P % i] * (P - P / i) % P;
}

int n; Poly A;

int main()
{
    cin >> n, A.resize(n), Pinit(4 * n + 2);

    return 0;
}
/*
 在 Z_{998244353} 中的一个多项式板子，包含：
 - 线性运算
 - NTT & 卷积
 - 求逆
 - 导数 & 积分
 - ln & exp
 - 开平方
 - 快速幂
 上面标了对于 n = 10^5 在洛谷上的表现（含 IO 耗时），不得不说常数还是不小的。
*/

终于完结啦！！！

（真的吗？）

坑点一览

场上打 vector 的多项式板子还是有些困难的，把寄过的点都写这吧。