多项式全家桶
数学推导
根据代数基本定理，\(n\) 个点值可以唯一的表示一个 \(n\) 阶多项式，快速傅里叶变换可以帮助我们快速求出一个n就多项式的点值表示。
对于一个 \(n\) 阶多项式 \(A(x)\) 我们想要得到它的点值表示，朴素的方法是用 \(n\) 个值代入计算。这样计算的时间复杂度是 \(O(n^2)\)，快速傅里叶变换的思路就是通过选取特殊的点来加速计算。
单位根
对于一个复平面上的单位圆，将它 \(n\) 等分，从实轴逆时针方向的第一份就是一个 \(n\) 次单位根，记作 \(w_n^1\)，这个单位根的 \(0\) 到 \(n-1\) 次方取遍了单位元上的 \(n\) 等分点，记作 \(w_n^k\)，它们具有以下性质：
(1)
\[\begin{align}
w_{2n}^{2k} &= w_n^k \\
可拓展为&:w_{mn}^{mk} = w_n^k
\end{align}
\]
(2)
\[w_n^{k + \frac{n}{2}} = -w_n^k
\]
(3)
\[w_n^{k + n} = w_n^k
\]
快速傅里叶变换
设 \(n~(n = 2^m)\) 阶多项式 \(A(x)\) 的系数组成的向量为 \((a_0,a_1,\dots,a_{n-1})\)
\[\begin{align}
A(x) &= a_0 + a_1x + a_2x^2 + \dots + a_{n-1}x^{n-1} \\
&= (a_0 + a_2x^2 + \dots + a_{n-2}x^{n-2}) + (a_1 + a_3x^3 + \dots + a_{n-1}x^{n-1}) \\
&= A_1(x^2) + xA_2(x^2)
\end{align}
\]
将 \(w_n^k\) 代入，设 \(0 \leq k < \frac{n}{2}\)
\[\begin{align}
A(w_n^k) &= A_1(w_n^{2k}) + w_n^kA_2(w_n^{2k}) \\
&= A_1(w_\frac{n}{2}^{k}) + w_n^kA_2(w_\frac{n}{2}^{k})
\end{align}
\]
\[\begin{align}
A(w_n^{k + \frac{n}{2}}) &= A_1(w_n^{2k + n}) + w_n^{k + \frac{n}{2}}A_2(w_n^{2k + n}) \\
&= A_1(w_\frac{n}{2}^{k}) - w_n^kA_2(w_\frac{n}{2}^{k})
\end{align}
\]
因为 \(k\) 和 \(k+\frac{n}{2}\) 取遍了 \(0\) 到 \(n-1\) 的整数，所以我们只需要计算 \(A_1(w_\frac{n}{2}^k)\) 和 \(A_2(w_\frac{n}{2}^{k})\) 即可直接得到 \(A(w_n^k)\) 和 \(A(w_n^{k+\frac{n}{2}})\)。这是一个递归计算的过程，时间复杂度为 \(O(n\log n)\)。
FFT
void FFT(complex<double> *a, int len, int flag) {
for (int i = 0; i < len; i++) {
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int i = 1; i < len; i <<= 1) {
complex<double> wn(cos(pi / i), flag * sin(pi / i));
for (int j = 0; j < len; j += (i << 1)) {
complex<double> w0(1,0);
for (int k = 0; k < i; k++, w0 *= wn) {
complex<double> x = a[j + k], y = w0 * a[j + k + i];
a[j + k] = x + y;
a[j + k + i] = x - y;
}
}
}
}
NTT
void NTT(int *a, int len, int flag) {
for (int i = 0; i < len; i++) {
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int i = 1; i < len; i <<= 1) {
int gn = qpow(flag ? G : Gi, (P - 1) / (i << 1));
for (int j = 0; j < len; j += (i << 1)) {
int g0 = 1;
for (int k = 0; k < i; k++, g0 = g0 * gn % P) {
int x = a[j + k], y = g0 * a[j + k + i] % P;
a[j + k] = (x + y + P) % P;
a[j + k + i] = (x - y + P) % P;
}
}
}
if (!flag) {
int inv = qpow(len, P - 2);
for (int i = 0; i < len; i++) a[i] = a[i] * inv % P;
}
}
多项式乘法
FFT版
void Multiply(double *a, double *b, double *c, int deg1, int deg2) { //deg为阶数
for (int i = 0; i <= deg1; i++) f[i].real(a[i]);
for (int i = 0; i <= deg2; i++) f[i].imag(b[i]);
int l = max((int)ceil(log2(deg1 + deg2 + 2)), 1), len = 1 << l;
for (int i = deg1 + 1; i < len; i++) f[i] = 0;
for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
FFT(f, len, 1);
for (int i = 0; i < len; i++) f[i] = f[i] * f[i];
FFT(f, len, -1);
for (int i = 0; i <= deg1 + deg2; i++) c[i] = f[i].imag() / 2.0 / len + 0.5;
}
NTT版
void Multiply(int *a, int *b, int *c, int deg1, int deg2) { //deg为阶数
for (int i = 0; i <= deg1; i++) f[i] = a[i];
for (int i = 0; i <= deg2; i++) g[i] = b[i];
int l = max((int)ceil(log2(deg1 + deg2 + 2)), 1), len = (1 << l);
for (int i = deg1 + 1; i < len; i++) f[i] = 0;
for (int i = deg2 + 1; i < len; i++) g[i] = 0;
for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
NTT(f, len, 1); NTT(g, len, 1);
for (int i = 0; i < len; i++) f[i] = f[i] * g[i] % P;
NTT(f, len, 0);
for (int i = 0; i <= deg1 + deg2; i++) c[i] = f[i];
}
多项式求逆
已知多项式 \(A(x)\) ，求多项式 \(B(x)\) ，使得 \(A(x)*B(x)\equiv 1~(mod~x^n)\) 。
若 \(n=1\) ，则 \(A^{-1}\) 即为常数项的逆元
假设已经求出了 \(B_0(x)\) ，使得
\[A(x)B_0(x)\equiv 1 (mod~x^{\lceil\frac{n}{2}\rceil})
\]
\[A(x)B(x)\equiv 1~(mod~x^n) \Rightarrow A(x)B(x)\equiv 1~(mod~x^{\lceil\frac{n}{2}\rceil})
\]
两式相减得
\[B(x)-B_0(x)\equiv 0~(mod~x^{\lceil\frac{n}{2}\rceil})
\]
平方后得
\[B^2(x)-2B(x)B_0(x)+B_0^2(x)\equiv 0~(mod~x^n)
\]
两边同时乘\(A(x)\) 得
\[A(x)B^2(x)-2A(x)B(x)B_0(x)+A(x)B_0^2(x)\equiv 0~(mod~x^n)
\]
\[B(x)-2B_0(x)+A(x)B_0^2(x)\equiv 0~(mod~x^n)
\]
\[B(x)\equiv 2B_0(x)-A(x)B_0^2(x)~(mod~x^n)
\]
递归求解即可
void Inverse(int *a, int *b, int deg) { //deg为长度
if (deg == 1) {
b[0] = qpow(a[0], P - 2);
b[1] = 0;
return;
}
Inverse(a, b, (deg + 1) >> 1);
int l = max((int)ceil(log2(deg << 1)), 1), len = (1 << l);
for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
for (int i = 0; i < deg; i++) a0[i] = a[i];
for (int i = deg; i < len; i++) a0[i] = 0;
NTT(a0, len, 1); NTT(b, len, 1);
for (int i = 0; i < len; i++) b[i] = (b[i] * 2 % P - a0[i] * b[i] % P * b[i] % P + P) % P;
NTT(b, len, 0);
for (int i = deg; i < len; i++) b[i] = 0;
}
多项式除法（取模）
已知 \(n\) 阶多项式 \(A(x)\) ， \(m\) 阶多项式 \(B(x)\) 求 \(n-m\) 阶多项式 \(C(x)\) 和 \(m-1\) 阶多项式 \(R(x)\) 使得
\[A(x)=B(x)C(x)+R(x)
\]
发现 \(x^nA(\frac{1}{x})=A_r(x)\) 相当于把 \(A(x)\) 的系数翻转，则
\[x^nA(\frac{1}{x})=x^nB(\frac{1}{x})C(\frac{1}{x})+x^nR(\frac{1}{x})
\]
\[x^nA(\frac{1}{x})=x^mB(\frac{1}{x})x^{n-m}C(\frac{1}{x})+x^{n-m+1}x^{m-1}R(\frac{1}{x})
\]
\[A_r(x)=B_r(x)C_r(x)+x^{n-m+1}R_r(x)
\]
对 \(x^{n-m+1}\) 取模
\[A_r(x)=B_r(x)C_r(x)~(mod~x^{n-m+1})
\]
\[C_r(x)=A_r(x)B_r^{-1}(x)~(mod~x^{n-m+1})
\]
需要用到多项式求逆，将 \(C_r(x)\) 翻转后得到 \(C(x)\) ，带回原式算出 \(R(x)\)
\[R(x)=A(x)-B(x)C(x)
\]
void Module(int *a, int *b, int *c, int *r, int deg1, int deg2) { //deg为阶数
if (deg1 < deg2) return;
for (int i = 0; i <= deg1; i++) ar[i] = a[deg1 - i];
for (int i = 0; i <= deg2; i++) br[i] = b[deg2 - i];
Inverse(br, ibr, deg1 - deg2 + 1);
Multiply(ar, ibr, c, deg1, deg1 - deg2);
reverse(c, c + deg1 - deg2 + 1);
Multiply(c, b, b, deg1 - deg2, deg2);
for (int i = 0; i <= deg1; i++) r[i] = (a[i] - b[i] + P) % P;
}
浙公网安备 33010602011771号