FFT / NTT
域和多项式的定义：
域的定义：
如果一个包含至少两个元素的集合 \(F\) 及两个运算 \(+\) 和 \(\cdot\) 满足：
- \((F,+)\) 是交换群，如果这个群的单位元为 \(0\)，那我们记 \(F^{\times}=F∖\{0\}\)；
- \((F^{\times},\cdot)\) 也是交换群，一般这个群的单位元被称为 \(1\)。
那我们就称 \((F,+,\cdot)\) 构成一个域。在不引起歧义的情况下，域的运算符默认是 \(+\) 和 \(\cdot\)。
比如我们常见的 \(\mathbb{Q},\mathbb{R},\mathbb{C}\) 都在常规的加法和乘法下构成一个域。
如果 \(p\) 是一个素数，那么模 \(p\) 的完全剩余系 \(\{0,1,\cdots,p-1\}\) 就在模 \(p\) 意义下的加法和乘法构成一个域，这个域常被记作 \(\mathbb{F}_p\) 或 \(\mathbb{Z}_p\)。
域上多项式的定义：
如果 \((F,+,\cdot)\) 是一个域，那么 \(F\) 上的多项式指的是形如 \(a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\) 的式子，其中 \(a_0,a_1,\cdots,a_n\in F\)。这里的 \(x\) 只是一个标识符，没有实际含义。
为方便，我们记 \(a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\) 为 \(\displaystyle\sum_{k=0}^n a_kx^k\)，记 \(F\) 上的所有多项式为 \(F[x]\)。
如果 \(f(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1}\)，定义最大的使 \(a_k\) 不为 \(0\) 的 \(k\) 为 \(f(x)\) 的次数，记做 \(\operatorname{deg} f\)。如果所有的 \(a_k\) 都是 \(0\)，则称 \(f(x)\) 为 \(0\) 多项式，记作 \(f(x)=0\)，此时补充定义 \(\operatorname{deg} f=-\infty\)。
域上多项式也可以带入具体的 \(x\)，也可以让多项式相加、相减、相乘，只不过运算实在 \(F\) 上进行的。
推广后的多项式仍然保持原来的一些良好性质，比如加法具有结合和交换律，乘法对加减法有分配律。
多项式乘法：
普通的多项式乘法是每个元素挨个相乘在相加，这样很慢，时间复杂度是 \(\mathcal{O}(n^2)\) 的。下面介绍一个更为快速的做法。
离散傅里叶变换 / DFT：
考虑换一种方式思考，我们知道 \(n+1\) 个横坐标不同的点确定一个 \(n\) 次多项式，那知道 \(n+1\) 个横坐标不同的点可以确定一个 \(n\) 次域上多项式吗？这是可以的，证明如下：
待定域上多项式 \(f(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\)，即证在域上的运算满足：
第一个 \(n\times n\) 的矩阵是范德蒙德矩阵，其行列式为 \(\displaystyle\prod_{0\le j<i\le n}(x_i-x_j)\)，其显然不为 \(0\)，因为 \(x\) 互不相同，由线性方程的求解理论证毕。\(\Box\)
后面我们会简称域上多项式为多项式。
考虑两个多项式 \(f(x),g(x)\in F[x]\)，设 \(f(x)\) 为 \(n\) 次多项式，\(g(x)\) 为 \(m\) 次多项式，则 \(f(x)g(x)\) 为 \(n+m\) 次多项式，由上述引理，我们只要知道 \(f(x_0)g(x_0),f(x_1)g(x_1),\cdots,f(x_{n+m})g(x_{n+m})\) 的值就可以确定 \(f(x)g(x)\) 这个多项式。
现在有两个困难。一：如何快速的通过 \(f(x)\) 的系数求解出 \(f(x_0),f(x_1),\cdots f(x_{n+m})\) 的值，即从 \(a\) 变到 \(y\)。二：如何快速的通过 \(f(x_0)g(x_0),f(x_1)g(x_1),\cdots,f(x_{n+m})g(x_{n+m})\) 求出 \(f(x)g(x)\) 的系数，即从 \(y\) 变到 \(a\)，其实也就是“一”的逆变换。
注意到，\(x_0,x_1,\cdots,x_{n+m}\) 我们可以任取，我们肯定想要 \(x\) 构成的范德蒙德矩阵的逆矩阵好求，这让我们想到了本原单位根！
我们设 \(N=n+m\)，\(\omega\) 为 \(N\) 次本原单位根，
则：
这样我们就不需要用 \(\mathcal{O}(N^3)\) 的时间高斯消元求逆矩阵了，这非常的棒！
我们设列向量 \(a=(a_0,\cdots,a_{N-1})^{\top}\)，记 \(a\) 的离散傅里叶变换（DFT）为 \(\mathcal{F}a=P\times a\)，\(a\) 的离散傅里叶逆变换（IDFT）为 \(\mathcal{F}^{-1}a=P^{-1}\times a\)，我们有 \(\mathcal{F}^{-1}(\mathcal{F}a)=a\)。
可以发现离散傅里叶变换就是通过系数求出了函数值，离散傅里叶逆变换就是通过函数值求出了系数，且我们通过定义可以发现其满足一个性质：\(\mathcal{F}f\cdot\mathcal{F}g=\mathcal{F}(f\ast g)\)（这个式子特别重要，许多的卷积都需要构造出一种变换满足这个性质），其中 \(\cdot\) 是对应位置相乘，\(\ast\) 是卷积（在这里就是多项式乘法）。所以我们可以通过如下方式求出 \(f\ast g\)：\(f\ast g=\mathcal{F}^{-1}(\mathcal{F}(f\ast g))=\mathcal{F}^{-1}(\mathcal{F}f\cdot\mathcal{F}g)\)，\(\cdot\) 运算可以 \(\mathcal{O}(N)\) 快速求，接下来的困难就是如何快速进行离散傅里叶变换和离散傅里叶逆变换。
虽然这个想法很天才，但是计算的时间复杂度并没有降低，进行矩阵乘法的时间复杂度还是 \(\mathcal{O}(N^2)\) 的。
快速傅里叶变换 / FFT：
分治写法：
FFT 是 DFT 的一种基于分治思想的实现方式，可以将计算离散傅里叶变换的复杂度降低到 \(\mathcal{O}(n\log n)\)（这里的 \(n\) 是两个多项式的次数之和）。
为了方便分治，不妨设 \(n=2^k\)，因为在多项式后面补 \(0\) 不会有影响。考虑 \(b=\mathcal{F}a\) 的第 \(k\) 个分量，\(b_k=\displaystyle\sum_{i=0}^{n-1}\omega^{ki}a_i=\sum_{i=0}^{\frac n2-1}(\omega^2)^{ki}a_{2i}+\omega^k\sum_{i=0}^{\frac n2-1}(\omega^2)^{ki}a_{2i+1}\)。注意到这个式子其实就是 \(b_k=(\mathcal{F}a')_{k\bmod \frac n2}+\omega^k(\mathcal{F}a'')_{k\bmod \frac n2}\)，其中 \(a'=(a_0,a_2,\cdots,a_{n-2})^{\top},a''=(a_1,a_3,\cdots,a_{n-1})\)。
注意这里离散傅里叶变换的本原单位根并不相同。
\(k\) 需要模 \(\frac n2\) 的原因是因为当 \(k\ge\frac n2\) 时，\((\mathcal{F}a')_k\) 和 \((\mathcal{F}a'')_k\) 没有意义，但是 \((\omega^2)^{\frac n2}=1\)（\(\omega^2\) 是 \(\frac n2\) 次本原单位根），\((\omega^2)^k=(\omega^2)^{k\bmod\frac n2}\)。
将上面的式子变一下：\(b_k=(\mathcal{F}a')_k+\omega^k(\mathcal{F}a'')_k,b_{k+\frac n2}=(\mathcal{F}a')_k-\omega^k(\mathcal{F}a'')_k,0\le k<\frac n2\)。（因为 \(\omega^{\frac n2}=-1\)）
递归到最后一层就是直接返回，因为 \(1\) 次单位根就是 \(1\)，所以此时 \(\mathcal{F}a=a\)。
代码实现：
void fft(int n, C *a, C g) { // C 是复数类型
if (n == 1) return;
static C tmp[kMaxN], t, x, y;
for (int i = 0; i < n / 2; i++) {
tmp[i] = a[i * 2], tmp[i + n / 2] = a[i * 2 + 1];
}
copy(tmp, tmp + n, a), fft(n / 2, a, g * g), fft(n / 2, a + n / 2, g * g), t = 1;
for (int i = 0; i < n / 2; i++, t = t * g) {
x = a[i], y = t * a[i + n / 2];
a[i] = x + y, a[i + n / 2] = x - y;
}
}
逆变换也很简单，逆变换其实就是将本原单位根 \(\omega\) 变成 \(\omega^{-1}\) 再跑一遍 FFT。
void idft(int n, C *a, C g) {
fft(N, a, {g.a, -g.b});
for (int i = 0; i < N; i++) {
a[i] = a[i] / N;
}
}
最终的代码就很容易写出来了。
#include <bits/stdc++.h>
using namespace std;
const int kMaxN = 4e6 + 5;
const double pi = acos(-1);
int n, m, N;
struct C { // 复数
double a, b;
C(double x = 0, double y = 0) { a = x, b = y; }
C operator+(const C &y) { return {a + y.a, b + y.b}; }
C operator-(const C &y) { return {a - y.a, b - y.b}; }
C operator*(const C &y) { return {a * y.a - b * y.b, a * y.b + b * y.a}; }
C operator*(const double &y) { return {a * y, b * y}; }
C operator/(const C &y) { return {(a * y.a + b * y.b) / (y.a * y.a + y.b * y.b), (y.a * b - a * y.b) / (y.a * y.a + y.b * y.b)}; }
C operator/(const double &y) { return {a / y, b / y}; }
} g, a[kMaxN], b[kMaxN];
void fft(int n, C *a, C g) {
if (n == 1) return;
static C tmp[kMaxN], t, x, y;
for (int i = 0; i < n / 2; i++) {
tmp[i] = a[i * 2], tmp[i + n / 2] = a[i * 2 + 1];
}
copy(tmp, tmp + n, a), fft(n / 2, a, g * g), fft(n / 2, a + n / 2, g * g), t = 1;
for (int i = 0; i < n / 2; i++, t = t * g) {
x = a[i], y = t * a[i + n / 2];
a[i] = x + y, a[i + n / 2] = x - y;
}
}
void idft(int n, C *a, C g) {
fft(N, a, {g.a, -g.b});
for (int i = 0; i < N; i++) {
a[i] = a[i] / N;
}
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> n >> m, n++, m++;
for (int i = 0; i < n; i++) {
cin >> a[i].a;
}
for (int i = 0; i < m; i++) {
cin >> b[i].a;
}
N = pow(2, ceil(log2(n + m))), g = C(cosl(2 * pi / N), sinl(2 * pi / N));
fft(N, a, g), fft(N, b, g);
for (int i = 0; i < N; i++) {
a[i] = a[i] * b[i];
}
idft(N, a, g);
for (int i = 0; i < n + m - 1; i++) {
cout << (int)(a[i].a + 0.5) << ' '; // 这里四舍五入是为了防止精度误差
}
return 0;
}
这份代码可以通过 P3803 【模板】多项式乘法（FFT）。
倍增写法：
倍增的写法比分治写法常数更小。
我们考虑分治中 \(a\) 是怎么变的：
- \((a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)\)
- \((a_0,a_2,a_4,a_6),(a_1,a_3,a_5,a_7)\)
- \((a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)\)
- \((a_0),(a_4),(a_2),(a_6),(a_1),(a_5),(a_3),(a_7)\)
我们设 \(r_i\) 为最后一次递归结束后 \(i\) 位置的元素是谁，通过观察我们可以发现 \(r_i\) 满足如下递推式：
其中 \(lg=\log_2n\)。
有了 \(r\)，我们就可以不用递归，快速计算最后一层的 \(a\)，然后像递归一样向上归并即可。
void init(int lg) {
for (int i = 0; i < 1 << lg; i++) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lg - 1));
}
}
void fft(int n, C *a, int type = 1) {
for (int i = 0; i < n; i++) {
if (i < r[i]) swap(a[i], a[r[i]]);
}
static C w, t, x, y;
for (int d = 2; d <= n; d *= 2) {
w = C(cosl(2 * pi / d), type * sinl(2 * pi / d));
for (int i = 0; i < n; i += d) {
t = C(1, 0);
for (int j = i; j < i + d / 2; j++, t = t * w) {
x = a[j], y = t * a[j + d / 2];
a[j] = x + y, a[j + d / 2] = x - y;
}
}
}
}
三次变两次优化：
由于 \((a+bi)^2=a^2-b^2+2abi\)，所以我们可以将 \(f\) 放在实部，\(g\) 放在虚部，那么平方后的多项式虚部除以二就是卷积后的结果。但是这样的精度误差比较大。
快速数论变换 / NTT：
由于 FFT 需要进行大量的复数运算，常数很大，并且需要用到浮点类型，精度也不高，所以我们引入快速数论变换 / NTT。NTT 其实就是在 \(\mathbb{F}_p\) 下的 FFT（\(p\) 是质数），但是复数如何取模呢？只需考虑本原单位根如何取模。想一想前面用到了本原单位根的什么性质？设 \(\omega_n\) 为 \(n\) 次单位根（\(n\) 为偶数）。
- \(\omega_n^n=1\)；
- \(1,\omega_n,\omega_n^2,\cdots,\omega_n^{n-1}\) 互不相同；
- \(\omega_n^{\frac n2}=-1\);
- \(\omega_{\frac n2}^{\frac k2}=\omega_n^k\)，其中 \(k\) 是偶数。
这让我们想到原根，设 \(g\) 是质数 \(p\) 的原根，我们令 \(g^{\frac{p-1}{n}}\) 为 \(n\) 次单位根，需要证明满足上述性质。
证明：
第一点：\((g^{\frac{p-1}{n}})^n=g^{p-1}\equiv 1\pmod p\)。
第二点：设 \(0\le k,l<n\) 且 \(k\neq l\)，假设存在 \(k,l\) 使 \(g^{\frac{k(p-1)}{n}}\equiv g^{\frac{l(p-1)}{n}}\pmod p\)。
则由原根的性质得：\(k\frac{p-1}{n}\equiv l\frac{p-1}{n}\pmod {p-1}\)。
然后就有 \(k\equiv l\pmod{\dfrac{p-1}{(p-1,\dfrac{p-1}{n})}}\)，即 \(k\equiv l\pmod n\)，又有 \(0\le k,l<n\)，则 \(k=l\)，矛盾！
第三点：\(((g^{\frac{p-1}{n}})^{\frac n2})^2=g^{p-1}\equiv 1\pmod p\)，由二次剩余的知识知 \((g^{\frac{p-1}{n}})^{\frac n2}=\pm 1\)，但若 \((g^{\frac{p-1}{n}})^{\frac n2}=1\)，则 \(g^{\frac{p-1}{2}}=1\)，与 \(g\) 是原根矛盾！所以 \((g^{\frac{p-1}{n}})^{\frac n2}=-1\)。
第四点：\((g^{\frac{p-1}{\frac n2}})^{\frac k2}=(g^{\frac{p-1}{n}})^k\)。\(\Box\)
这就说明了我们可以用原根替代本原单位根，我们可以取 \(g^{\frac{p-1}{n}}\) 为 \(n\) 次单位根，这需要保证 \(n\mid (p-1)\)，其中 \(n\) 是 \(2\) 的幂，我们可以让 \(p\) 取 \(998244353=119\times 2^{23}+1\)，就可以满足上述条件，而且记住 \(998244353\) 的原根为 \(3\)。这样我们就可以做 NTT 了，代码如下：
using namespace std;
const int kMaxN = 4e6 + 5, kM = 998244353;
int n, m, a[kMaxN], b[kMaxN], r[kMaxN];
long long fpow(long long a, long long b, long long p, long long r = 1) {
for (a = (a % p + p) % p; b; b & 1 && (r = r * a % p), a = a * a % p, b >>= 1) {
}
return r % p;
}
void init(int lg) {
for (int i = 0; i < 1 << lg; i++) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lg - 1));
}
}
void ntt(int n, int *a, int g) {
int lg = log2(n);
for (int i = 0; i < n; i++) {
if (i < r[i]) swap(a[i], a[r[i]]);
}
for (int d = 2, w; d <= n; d *= 2) {
w = fpow(g, n / d, kM);
for (int i = 0; i < n; i += d) {
for (int j = i, t = 1, x, y; j < i + d / 2; j++, t = 1LL * t * w % kM) {
x = a[j], y = 1LL * t * a[j + d / 2] % kM;
a[j] = (x + y) % kM, a[j + d / 2] = (x - y + kM) % kM;
}
}
}
}
void idft(int n, int *a, int g, int invn = 0) {
ntt(n, a, fpow(g, kM - 2, kM)), invn = fpow(n, kM - 2, kM);
for (int i = 0; i < n; i++) {
a[i] = 1LL * a[i] * invn % kM;
}
}
int main() {
cin >> n >> m, n++, m++;
for (int i = 0; i < n; i++) {
cin >> a[i];
}
for (int i = 0; i < m; i++) {
cin >> b[i];
}
int N = pow(2, ceil(log2(n + m))), g = fpow(3, (kM - 1) / N, kM);
init(__lg(N)), ntt(N, a, g), ntt(N, b, g);
for (int i = 0; i < N; i++) {
a[i] = 1LL * a[i] * b[i] % kM;
}
idft(N, a, g);
for (int i = 0; i < n + m - 1; i++) {
cout << a[i] << ' ';
}
return 0;
}
这是 NTT 中可能会用到的素数和其原根，有需要的同学可以看一看：https://blog.miskcoo.com/2014/07/fft-prime-table。
FFT 和 NTT 都需要理解透彻，这样可以帮助你卡常。
多项式牛顿迭代：
多项式牛顿迭代过程：
多项式牛顿迭代可以用来处理下面的问题。
已知函数 \(G(x)\)，求一个多项式 \(F(x)\bmod x^n\)，满足方程 \(G(F(x))\equiv 0\pmod {x^n}\)。
当 \(n=1\)，需要单独求出满足 \(G(F(x))\equiv 0\pmod x\) 的 \(F(x)\)。
接下来假设我们已经求出了 \(F_0(x)\) 使得 \(G(F_0(x))\equiv 0\pmod {x^{\lceil\frac n2\rceil}}\)，考虑 \(G(F(x))\) 在 \(F_0(x)\) 上的泰勒展开：
由于 \(F_0(x)\equiv F(x)\pmod {x^{\lceil\frac n2\rceil}}\)，所以有 \((F(x)-F_0(x))^2\equiv 0\pmod {x^{2\lceil\frac n2\rceil}}\)，也就有 \((F(x)-F_0(x))^2\equiv 0\pmod {x^n}\)，更高次同理。
所以就有 \(G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\pmod {x^n}\)，由于 \(G(F(X))\equiv 0\pmod {x^n}\)，转化得：\(F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n}\)。
有了这个式子，这道题就做完了。
下面我们来看几道例题：
多项式求逆：
给定一个多项式 \(F(x)\)，请求出一个多项式 \(G(x)\)，满足 \(F(x)\ast G(x)\equiv 1\pmod {x^n}\)。
设 \(H(G(x))=\dfrac{1}{G(x)}-F(x)\)，我们需要求出其零点多项式。
当 \(n=1\) 时，\(G(x)\equiv \dfrac{1}{F(x)}\pmod x\)，也即 \(F(x)\) 常数项的逆元。
假设已经求出满足 \(H(G_0(x))\equiv 0\pmod {x^{\lceil\frac n2\rceil}}\) 的 \(G_0(x)\)。因为 \(G'(F(x))=-\dfrac{1}{G(x)^2}\)（可以将 \(F(x)\) 看作常数项），所以 \(G(x)\equiv G_0(x)-\dfrac{\dfrac{1}{G_0(x)}-F(x)}{-\dfrac{1}{G_0(x)^2}}\pmod {x^n}\)，也即 \(G(x)\equiv G_0(x)(2-G_0(x)F(x))\pmod {x^n}\)，在每一层用多项式乘法即可。
其时间复杂度为 \(T(n)=T(\dfrac n2)+\mathcal{O}(n\log n)=\mathcal{O}(n\log n)\)。
代码：https://www.luogu.com.cn/record/198798739。
多项式 exp：
给出 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod{\:x^n}\) 下的多项式 \(B(x)\)，满足 \(B(x) \equiv \text e^{A(x)}\)。这里保证 \([x^0]A(x)=0\)。
设 \(F(B(x))=\ln B(x)-A(x)\)，则 \(F'(B(x))=\dfrac{1}{B(x)}\)，我们需要求出其零点多项式。
当 \(n=1\) 时，\(B(x)=1\)。
假设已经求出满足 \(F(B_0(x))\equiv 0\pmod {x^{\lceil\frac n2\rceil}}\) 的 \(B_0(x)\)。那么 \(B(x)\equiv B_0(x)-\dfrac{\ln B_0(x)-A(x)}{\dfrac{1}{B_0(x)}}\pmod x^n\)，即 \(B(x)\equiv B_0(x)(1-\ln B_0(x)+A(x))\pmod x^n\)，在每一层用多项式乘法和多项式 ln 即可。
其时间复杂度为 \(T(n)=T(\dfrac n2)+\mathcal{O}(n\log n)=\mathcal{O}(n\log n)\)。
这里简单说一下为什么要保证 \([x^0]A(x)=0\)，其实应该是 \([x^0]A(x)\equiv 0\pmod p\)，其中 \(p\) 为模数。考虑 \(B(x)\) 的常数项，令 \(k=[x^0]A(x)\)，则 \([x^0]B(x)=\displaystyle\sum_{i=0}^{\infty}\dfrac{k^i}{i!}\)，然而 \(p!\) 在模 \(p\) 意义下是没有逆元的，且如果不考虑 \(p\) 的倍数，这个级数也并不收敛。（因为有无限个大于 \(0\) 的数，且没有小于 \(0\) 的数）
这也说明了我们需要注意多项式多项式变换的常数项（后面的求导再积分也是）。
代码：https://www.luogu.com.cn/record/198806655。
多项式牛顿迭代还可以求多项式开根，原理都是一样的，这里就不再赘述。
多项式初等函数：
多项式求逆：
使用多项式牛顿迭代可解决。
多项式开方：
使用多项式牛顿迭代可解决。
多项式 ln：
由 \(\ln\) 的定义知，\([x^0]f(x)=1\)，否则没有定义，所以 \([x^0]\ln f(x)=0\)。
考虑对 \(\ln f(x)\) 求导后再积分：
多项式求导和积分很好处理，再用多项式求逆可解决。
注意 \(\ln f(x)\) 的常数项为 \(0\)。
多项式 exp：
使用多项式牛顿迭代可解决。
多项式三角函数：
考虑使用欧拉公式 \(e^{ix}=\cos x+i\sin x\)。
那么就有：
可以直接使用多项式 exp 求解（\(i\) 其实就是 \(\omega_4\)）。
而其他的三角函数可以直接用 \(\sin,\cos\) 表示，都可以用多项式求逆得到。
多项式反三角函数：
这里同样需要保证 \([x^0]f(x)=0\)，
将 \(\arcsin,\arccos,\arctan\) 求导再积分可得：
这里需要注意一下 \([x^0]\arcsin f(x)=0,[x^0]\arccos f(x)=\dfrac\pi 2,[x^0]\arctan f(x)=0\)，所以 \(\arccos f(x)\) 需要保证 \([x^0]\arccos f(x)=1\) 才有定义。
多项式快速幂：
\(f(x)^k=\exp(k\ln f(x))\)，但是这里需要保证 \([x^0]f(x)=1\)。
如果 \([x^0]f(x)\) 不等于 \(1\) 且不等于 \(0\)，我们可以将整个多项式除以常数项最后再乘上。
如果 \([x^0]f(x)\) 等于 \(0\)，我们可以将 \(x^k\) 提出最后再乘上，其中 \(k\) 是满足 \([x^k]f(x)\neq 0\) 的最小的 \(k\)。
时间复杂度：\(\mathcal{O}(n\log n)\)（可是它又时常数太大，比 \(\mathcal{O}(n\log n\log k)\) 的普通快速幂还慢）
多项式的简单应用：
这里需要注意到 \(A\) 和 \(B\) 的卷积为 \((A\ast B)[k]=\displaystyle\sum_{i=0}^k A_iB_{k-i}\)。
P5488 差分与前缀和：
题目链接：https://www.luogu.com.cn/problem/P5488。
给定一个长度为 \(n\) 的序列 \(a\)，求出其 \(k\) 阶差分或前缀和。
设 \(f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}\)，容易发现差分一次相当于乘 \((1-x)\)，则 \(k\) 阶差分就是将序列乘以 \((1-x)^k\)，\(k\) 阶前缀和就是将序列乘以 \((1-x)^{-k}\)，考虑其泰勒展开式，容易发现 \(k\) 可以对模数取模，用多项式快速幂容易维护。
另一种简单常数小的方法是将 \((1-x)^k\) 以及 \((1-x)^{-k}\) 用广义二项式定理展开，由于之需求多项式 \(\bmod x^n\) 的值，可以舍去 \(n\) 以后的多项式，组合数系数可以递推求，两边多项式乘法即可。
P3338 [ZJOI2014] 力：
题目链接：https://www.luogu.com.cn/problem/P3338。
给出 \(n\) 个数 \(q_1,q_2, \dots q_n\)，定义 \(E_j~=~\displaystyle\sum_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}~-~\sum_{i = j + 1}^{n}\dfrac{q_i}{(i-j)^2}\)，求 \(E_i\)。
设：
则 \(E\) 就是 \(p\ast q\)。
