前置芝士
复数指数幂
有欧拉公式
显然$e^{i\pi}=1$
单位根
单位根的定义：
在复数域下，满足$x^n=1$的$x$被称为$n$次单位根。
根据代数基本定理，$n$次单位根一共有$n$个，且所有的$n$次单位根按照幅角排列是均匀分布的，第$k(0\le k<n)$个$n$次单位根为
$n$个单位根在复平面上平分单位圆。
本原单位根
$0$到$n-1$次方的值能生成所有$n$次单位根的$n$次单位根称为$n$次本原单位根。
记$n$次本原单位根为$\omega_n=e^{i\frac{2\pi}{n}}=\cos\dfrac{2\pi}{n}+i\sin\dfrac{2\pi}{n}$
设$n=2m$，有
阶与原根
欧拉定理
若正整数$a,m$满足$(a,m)=1$，则使得$a^n\equiv 1\pmod{m}$的最小正整数$n$称为$a$模$m$的阶，记作$\delta_m(a)$。
若$\delta_m(a)=\phi(m)$，则称$a$为$m$的一个原根。
阶的性质：
- $a^0,a^1,…,a^{\delta - 1}$在模$m$意义下两两不同
- $a^{\gamma}\equiv a^{\gamma’}\pmod{m}\Leftrightarrow \gamma\equiv\gamma’\pmod{\delta}$
- $\delta |\phi(m)$
原根的存在定理：
只有模$2,4,p^a,2p^a$存在原根，其中$p$为奇质数。
原根的判定定理：
设$m>1$，$g$为正整数且$(m,g)=1$，则$g$是$m$的原根当且仅当对于任意$\phi(m)$的质因子$q_i$有$g^{\frac{\phi(m)}{q_i}}\not\equiv 1\pmod{m}$。
正整数$n$若存在原根，则其最小原根是$O(n^{\frac{1}{4}})$的，要求出$n$的所有原根可以通过先求出最小原根$g$，根据$g$来求所有原根。所有原根应满足$g’=g^k$且$(k,\phi(n))=1$。
1 | ll getrt(ll p) { |
指标
对于质数$p$，假设$g$是$p$的一个原根，则$g^0,g^1,…,g^{p-2}$在模$p$意义下是$1,2,…,p-1$的一个排列。假设对于$1\le x\le p-1$有$g^c\equiv x\pmod{p}$，则称$x$的指标为$c$，记作${\rm ind}(x)=c$。
指标的性质类似于对数，对于$\forall x,y\in[1,p-1]$：
- ${\rm ind}(x\cdot y)={\rm ind}(x)+{\rm ind}(y)$
- ${\rm ind}(x^c)=c\cdot {\rm ind}(x)$
求指标直接用BSGS。
离散快速傅里叶变换（FFT）
两个多项式做乘法，用系数表示法直接算显然需要$O(n^2)$时间复杂度，于是我们考虑用点值表示法来计算，如果我们能求出两个多项式的点值表示法，那么就可以直接相乘来得到多项式相乘后的点值表示法，注意由于两个$n$次多项式相乘后是$2n$次多项式，所以我们需要求$2n$个点值。
点值表示法如果直接求，需要$O(n^2)$计算两个多项式的点值，$O(n^3)$高斯消元解出系数，我们需要将这两步的复杂度优化。
求出一个$n$次多项式在每个$n$次单位根下的点值的过程被称为离散傅里叶变换（DFT），将点值重新插值成系数表示法的过程称为逆离散傅里叶变换（IDFT）。
下面我们设进行变换的多项式为$n-1$次多项式$A(x)=\sum\limits_{i=0}^{n-1}a_i\times x^i$，其中$n$是2的幂（若不是则补齐），$m=\dfrac{n}{2}$。
FFT
考虑对$A(x)$的下角标按奇偶分类，得到
将右式提一个x出来，得到
我们设
则
虽然规模减小了，但是$A_0,A_1$系数是不一样的，那每层往下递归两次还是$O(n^2)$的。
由前置芝士我们知道，$\omega_n^{m+k}=-\omega_n^{k}$，$(\omega_n^k)^2=\omega_m^k$，于是就来考虑能不能只求小于$m$的部分。
对于$0\le k< m$，我们有
考虑$m+k$的部分，我们有
观察(2)和(5)，只要求出$0\le k<m$的部分的$A_0,A_1$，我们就能求出$A$了，问题规模缩小到原来的$\dfrac{1}{2}$。到这里，由系数表示法求点值的问题就解决了。
复杂度$O(n\log n)$
IDFT
现在我们来看，得到了$n$个点值后如何插值得到多项式的系数。
假设已知一个$n-1$次多项式$A(x)=\sum\limits_{i=0}^{n-1}a_ix^i$进行DFT后得到了点值$\{b_k\}$，那么求解$a$的过程相当于求解以下方程组
设为$\Omega\times A=B$，则$B=\Omega^{-1}\times A$
这里直接给出结论
令$M=\Omega^{-1}\times\Omega$，写出$M_{ij}$的表达式
于是求解$a$的方程就为
这就相当于给定了$B(x)=\sum\limits_{i=0}^{n-1}b_ix^i$，求点值$B(\omega_n^{-k}),0\le k<n$
该问题的形式与前面求$A(x)$的点值完全相同，唯一的区别就是$\omega_n$的指数多了一个负号，计算的时候把符号加上就行了。
时间复杂度也是求点值的复杂度$O(n\log n)$。
模板
以上只是利用FFT计算多项式乘法的基本原理，FFT本身常数很大，优化的手段很多，例如蝴蝶变换可以非递归实现FFT，还有两个多项式都需要求点值，这就需要三次FFT，而实际上可以一次FFT求出两个多项式的点值，下面的板子是最终版本，可以直接使用。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42struct Complex {
double x, y;
Complex(double x = 0.0, double y = 0.0) : x(x), y(y) {}
Complex operator +(const Complex &s)const { return Complex(x + s.x, y + s.y); }
Complex operator -(const Complex &s)const { return Complex(x - s.x, y - s.y); }
Complex operator *(const Complex &s)const { return Complex(x * s.x - y * s.y, x * s.y + y * s.x); }
};
struct Fast_Fourier_Transfrom {
Complex a[N];
int r[N];
void FFT(int len, Complex a[], int type) {
for (int i = 0; i < len; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int mid = 1; mid < len; mid <<= 1) {
Complex wn = Complex(cos(PI / mid), type * sin(PI / mid));
for (int R = mid << 1, j = 0; j < len; j += R) {
Complex w = Complex(1, 0);
for (int k = 0; k < mid; ++k, w = w * wn) {
Complex x = a[j + k], y = w * a[j + mid + k];
a[j + k] = x + y;
a[j + k + mid] = x - y;
}
}
}
}
void Mul(int n, int m) {
int len = 1, x = 0;
//假设n>=m，将第一个多项式的系数放在实部，第二个多项式的系数放在虚部
//最终a[i].x就是指数为i的项的系数
for (; len <= n + m - 2; len <<= 1, ++x);
for (int i = 0; i < len; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (x - 1));
memset(a + n, 0, sizeof(Complex) * (len - n + 1));
FFT(len, a, 1);
for (int i = 0; i <= len; ++i)
a[i] = a[i] * a[i];
FFT(len, a, -1);
for (int i = 0; i <= n + m - 2; ++i)
a[i].x = a[i].y / len / 2;
}
}FFT;
快速数论变换（NTT）
假设质数$p$满足$p=r2^l+1$，$g$是$p$的原根，用$g_n=g^{\frac{p-1}{n}}$替换FFT中的$\omega_n$。
当项数$n\le 2^l$时，与$\omega_n$相似的性质仍然成立：
- $g_{2n}^{2k}\equiv g_n^k\pmod{p}$
- $g_{2n}^n\equiv-1\pmod{p}$
- $\sum\limits_{k=0}^{n-1}g_{n}^{ik}g_n^{-kj}=\begin{cases}n&i=j\\0&i\not=j\end{cases}\pmod{p}$，其中$0\le i,j<n$
所以将$\omega_n$直接替换成$g_n$，DFT、IDFT的过程仍然正确，这样构成的算法就是NTT。值得注意的是，用NTT求解的多项式系数是模$p$后的结果，所以如果需要精确的系数，就需要选用较大的模数。
常见模数：
- $65537=2^{16}+1,g=3$
- $167772161=5\times2^{25}+1$
- $469762049=7\times2^{26}+1$
- $998244353=119\times 2^{23}+1,g=3$
- $1004535809=479\times2^{21}+1,g=3$
- $4179340454199820289=29\times 2^{57}+1,g=3$
模板
仿照FFT直接写出NTT的模板。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45struct Number_Theoretic_Transform {
ll G = 3, Gi = 332748118, mod = 998244353;
//特别注意，封装mod常数非常大，如果不是为了用MTT还是放外面好
ll a[N], b[N];
int r[N];
ll FastPow(ll a, ll b) {
ll ans = 1;
for (; b; b >>= 1, (a *= a) %= mod)
if (b & 1)(ans *= a) %= mod;
return ans;
}
ll ADD(ll a, ll b) { a += b; return a >= mod ? a - mod : a; }
ll DEL(ll a, ll b) { a -= b; return a < 0 ? a + mod : a; }
void NTT(int len, ll a[], int type) {
for (int i = 0; i < len; ++i)
if (i < r[i])
swap(a[i], a[r[i]]);
for (int mid = 1; mid < len; mid <<= 1) {
ll gn = FastPow(type ? G : Gi, (mod - 1) / (mid << 1));
for (int R = mid << 1, j = 0; j < len; j += R) {
ll g = 1;
for (int k = 0; k < mid; ++k, g = g * gn % mod) {
ll x = a[j + k], y = g * a[j + mid + k] % mod;
a[j + k] = ADD(x, y);
a[j + mid + k] = DEL(x, y);
}
}
}
}
void Mul(int n, int m) {
int len = 1, x = 0;
for (; len <= n + m - 2; len <<= 1, ++x);
for (int i = 0; i < len; ++i)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (x - 1));
memset(a + n, 0, sizeof(ll) * (len - n + 1));
memset(b + m, 0, sizeof(ll) * (len - m + 1));
NTT(len, a, 1); NTT(len, b, 1);
for (int i = 0; i < len; ++i)
a[i] = a[i] * b[i] % mod;
NTT(len, a, 0);
ll inv = FastPow(len, mod - 2);
for (int i = 0; i <= n + m - 2; ++i)
a[i] = a[i] * inv % mod;
}
}NTT;
任意模数NTT（MTT）
在NTT的基础上我们可以推广出求解模数为任意值的多项式运算方法。
其基本思想为选取若干个不同的NTT模数，通过NTT分别求出运算结果，最后利用exCRT将这些结果合并，由exCRT的性质我们知道最终合并的结果是在模NTT模数之积意义下的，所以需要保证NTT模数之积是大于运算结果的，这样就能求出运算结果的实际值，当然任意模意义下的值也能求了。
这里以取三个模数做多项式乘法为例。
设某一项系数待求答案为$x$，三个模数为$A,B,C$，则有
合并前两个方程得到
这就意味着$x\equiv x_0+k_0A\pmod{AB}$，令$x_3=x_1+k_0A$，继续合并得
于是最终得到
1 | struct Multiple_Theoretic_Transform { |
其他运算也是类似的先分别求答案再合并。
拉格朗日插值
拉格朗日插值定理
$n$个点值$(x_i,y_i),(1\le i\le n)$，满足$x_i\not=x_j,(i\not=j)$，它们唯一确定一个$n-1$次多项式$f(x)$
如果已知$x_i=i$的连续$n$个点值，就可以$O(n)$地求出任意一个点值$f(x_0)$。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25ll solve(int x) {
//求出f(x)，已知x=0~n的点值
//处理分母，这里逆元可以递推，偷懒就快速幂
muln[1] = mod - 1, mulp[1] = 1;
for (int i = 2; i <= n + 1; ++i) {
muln[i] = (mod - muln[i - 1] * ksm(i, mod - 2) % mod) % mod;
mulp[i] = mulp[i - 1] * ksm(i,mod - 2) % mod;
}
//处理分子
L[0] = x, R[n] = x - n;
for (int i = 1; i < n; ++i) {
L[i] = L[i - 1] * ((x - i + mod) % mod) % mod;
R[n - i] = R[n - i + 1] * ((x - n + i + mod) % mod) % mod;
}
ll ans = 0;
for (int i = 0; i <= n; ++i) {
if (i == 0)
(ans += f[i] * R[1] % mod * muln[n] % mod) %= mod;
else if (i == n)
(ans += f[i] * L[n - 1] % mod * mulp[n] % mod) %= mod;
else
(ans += f[i] * L[i - 1] % mod * R[i + 1] % mod * muln[n - i] % mod * mulp[i] % mod) %= mod;
}
return ans;
}
牛顿迭代
给定多项式$g(x)$，求满足$g(f(x))=0$的形式幂级数$f(x)$。
牛顿迭代用于求解$f(x)$前$n$项。
- $n=1$时，解$g(a_0)=0$得到$a_0$
- 如果已知前$n$项$f(x)\equiv f_0(x)=a_0+a_1x+…+a_{n-1}x^{n-1}\pmod{x^n}$
- 则$f(x)\equiv f_0-\dfrac{g(f_0(x))}{g’(f_0(x))}\pmod{x^{2n}}$
多项式（形式幂级数）求逆
设$h(x)$是给定的形式幂级数，求它的逆$f(x)$，则$g(f(x))=\dfrac{1}{f(x)}-h(x)=0$，得到的迭代形式为
1 | void Inv(int n, ll f[], ll h[]) { |
多项式（形式幂级数）开方
设$h(x)$是给定的形式幂级数，求它的逆$f(x)$，则$g(f(x))=f^2(x)-h(x)=0$，得到的迭代形式为
1 | ll invf[N]; |
多项式（形式幂级数）对数
考虑$(\ln(f(x)))’=f’(x)\dfrac{1}{f(x)}$，求导很明显可以$O(n)$，求积分也是$O(n)$，求逆$O(n\log n)$，因此总复杂度是$O(n\log n)$。1
2
3
4
5
6
7
8
9
10void Ln(int n, ll f[], ll h[]) {
//求ln(h)，答案放在f中，n为项数
Inv(n, f, h);
for (int i = 0; i < n - 1; ++i)b[i] = h[i + 1] * (i + 1) % mod;
for (int i = 0; i < n; ++i)a[i] = f[i];
Mul(n, n - 1);
f[0] = 0;
for (int i = 1; i < n; ++i)
f[i] = a[i - 1] * FastPow(i, mod - 2) % mod;
}
多项式（形式幂级数）指数
设$h(x)$是给定的形式幂级数，满足$[x^0]h(x)=0$，设$f(x)=\exp(h(x))$，则有$g(f(x))=\ln(f(x))-h(x)=0$，得到迭代形式为
1 | ll lnf[N]; |
多项式（形式幂级数）快速幂
求$f^k(x)\pmod{x^n}$，等同于求$\exp(k\ln f)\pmod{x^n}$，当然这里需要保证$[x^0]f(x)=1$。1
2
3
4
5
6void Pow(int n, ll f[], ll h[], ll k) {
//求h^k，答案放在f中，n为项数
Ln(n, f, h);
for (int i = 0; i < n; ++i)h[i] = f[i] * k % mod;
Exp(n, f, h);
}