学习笔记:FFT与拉格朗日插值
多项式的表示形式
系数表示与点值表示
假设 \(f(x)\) 是一个 \(n\) 次多项式,则 \(f(x)\) 的系数表示为 \(f(x) = a_nx^n + a_{n - 1}x^{n - 1} + \cdots + a_0\)
\(f(x)\) 的点值表示为 \((x_0, f(x_0)), \ (x_1, f(x_1)), \dots, (x_n, f(x_n))\),其中 \(\forall i \neq j, \ x_i \neq x_j\)。
-
\(n + 1\) 个点值可以表示一个 \(n\) 次多项式。
-
在点值表示下 \(n\) 次多项式的乘法复杂的为 \(O(n)\)。
设 \(h(x) = f(x) \cdot g(x)\),则 \(\forall i \in [0, n], \ h(x_i) = f(x_i) \cdot g(x_i)\)。
复数与单位根
复数的指数形式
\(a + bi = re^{i\theta}\),其中 \(r = \sqrt {a^2 + b^2}, \ \tan \theta = \dfrac{b}{a}\)。
欧拉公式:\(e^{ix} = \cos x + i\sin x\)。
单位根
\(x^n = 1\) 在复数域上的根称为 \(n\) 次单位根。\(n\) 次单位根有 \(n\) 个,形式为 \(\omega_{n}^k = e^{i\frac{2k\pi}{n}}\)。
单位根在复平面上等分单位圆。
单位根的性质:
- \(\omega_n^k = \omega_{2n}^{2k}\)。
- \(\omega_{2n}^{k + n} = -\omega_{2n}^{k}\)
快速傅里叶变换(Fast Fourier Transform)
离散傅里叶变换(Discrete Fourier Transform)
将多项式 \(A(x) = a_0 + a_1x + \dots + a_{n - 1}x^{n - 1}\) 转化为其点值形式 \((\omega_n^k, \ A(\omega_n^k)), \ (k = 0, 1, \dots, n - 1)\)。
不妨令 \(n\) 为 \(2\) 的整次幂。
把原式拆成奇偶两部分,即 \(A(x) = (a_0 + a_2x^2 + \cdots + a_{n - 2}x^{n - 2}) + (a_1 + a_3x^3 + \cdots + a_{n - 1}x^{n - 1})\)。
令 \(B(x) = a_0 + a_2x + \cdots + a_{n - 2}x^{n / 2 - 1}, \ \ C(x) = a_1 + a_3x + \cdots + a_{n - 2}x^{n / 2 - 1}\)。
则 \(A(x) = B(x^2) + xC(x^2)\)。
对于 \(k \in [0, n / 2 - 1]\)
对于另一半的点值,可用 \(A(\omega_n^{k + n / 2})\) 表示,即
具体的讲,当 \(n = 8\) 时,各项存在如下关系:
\(T(n) = 2T(n / 2) + O(n)\)。
可以直观看出,以这种方式将系数表示转化为点表示的时间复杂度为 \(O(n \log n)\)。
对于具体实现
-
分治(递归,常数大)。
-
蝴蝶变换(bit-reversal permutation)(非递归,常数小)。
观察上图第一行和最后一行各项的二进制表示,互相颠倒。
因此可以预处理最后一行的系数,自下而上递推。
逆离散傅里叶变换(Inverse DFT)
将多项式的点值表示 \((\omega_n^k, \ b_k), \ (k = 0, 1, \dots, n - 1)\) 转化为其系数表示 \(A(x) = a_0 + a_1x + \dots + a_{n - 1}x^{n - 1}\)。
设 \(n \times n\) 的矩阵 \(\Omega\),其中 \(\Omega_{i, j} = \omega_n^{ij}\),设向量 \(a = (a_0, a_1, \dots, a_{n - 1}), \ b = (b_0, _1, \dots, b_{n - 1})\)。
则 IDFT 相当于求解方程 \(\Omega a = b\)。
如果我们能够求出 \(\Omega\) 的逆,则 \(a = \Omega^{-1}b\)。
设
其中
相当于等比数列求和
令 \(\overline\Omega\) 满足 \(\overline\Omega_{i, j} = \omega^{-ij}\),有 \(\overline\Omega\cdot\Omega = nI\),因此
相当于给定 \(B(x) = b_0 + b_1x + \cdots + b_{n - 1}x^{n - 1}\),求点值 \(a_k = B(\omega_n^{-k}), \ (0 \le k < n)\)。
例题
A * B Problem Plus
题意:给定 \(a, \ b\),求 \(a \times b\),\(a, \ b < 10^{50001}\)(范围可以更大)。
将 \(a\) 看作 \(a_0 + a_110^1 + a_210^2 + \cdots\),fft 后再处理进位。
SP8372 TSUM
题意:给定长度为 \(N\) 的数列 \(s\),对于任意可能存在的 \(V\),求满足 \(s_i + s_j + s_k = V, \ \ i < j < k\) 的对数,\(|s_i| \le 20000, \ N \le 40000\)。
构造多项式 \(A(x) = \sum x^{s_i}\)。
那么 \(\sum\limits_{i}\sum\limits_{j}\sum\limits_{k}[s_i + s_j + s_k = V] = [x^V]A^3(x)\)。
这样求得的数对不一定满足偏序关系,也不一定互不相同。
偏序关系很好处理,最后将答案除上 \(3!\) 即可。
现要除去存在位置相同的数对,考虑容斥。
性质 \(a_1: i = j\),性质 \(a_2: i = k\),性质 \(a_3:j = k\)。
-
减去满足三个性质中一个的方案。
钦定两个位置相等。
令 \(B(x) = \sum x^{2s_i}\),产生的方案为 \([x^V](A(x)B(x))\)。
-
加上满足三个性质中两个的方案。
\(i = j = k\)。
令 \(C(x) = \sum x^{3s_i}\),这部分的方案为 \([x^V]C(x)\)。
-
减去满足三个性质中三个的方案。
与第二部分相同。
则最终答案 \(D(x) = A^3(x) - 3B(x) + 2C(x)\)。
可以先对 \(A(x), \ B(x), \ C(x)\) 分别做 dft,计算出 \(D(x)\) 的点值,最后再对 \(D(x)\) 做 idft。
FFT 在字符串匹配中的应用
普通的单模式串匹配
问题概述:给定模式串 \(a\)(长度为 \(m\))、文本串 \(b\)(长度为 \(n\)),求所有位置 \(x\),满足 \(B\) 串以 \(x\) 结尾的连续 \(m\) 个字符与 \(a\) 串完全相同。
定义匹配函数 \(C(x, y) = [a_x - b_y]^2\),当 \(C(x, y) = 0\) 时则称 \(a\) 的第 \(x\) 个字符与 \(B\) 的第 \(y\) 个字符匹配。
定义完全匹配函数 \(P(x) = \sum\limits_{i = 0}^{m - 1}[a_i - b_{(x -m + i + 1)}]^2\)。
当且仅当 \(P(x) = 0\) 时以 \(x\) 结尾的连续 \(m\) 个字符与 $$ 匹a配。
令 \(s\) 串满足 \(s_i = a_{(m - i - 1)}\),即 \(a\) 串的翻转。
于是 \(P(x) = \sum\limits_{i = 0}^{m - 1}[s_{(m - i + 1)} - b_{(x -m + i + 1)}]^2\)。
暴力展开:\(P(x) = \sum\limits_{i = 0}^{m - 1} s_{(m - i + 1)}^2 + \sum\limits_{i = 0}^{m - 1}b_{(x -m + i + 1)}^2 - \sum\limits_{i = 0}^{m - 1}2\cdot s_{(m - i + 1)}b_{(x - m + i + 1)}\)。
第一项为常数。
第二项可以预处理 \(f(x) = \sum_{i = 0}^xb_i^2\),计算 \(f(x) - f(x - m)\)。
第三项等效于求 \(g(x) = \sum\limits_{i + j = x}s_ib_j\)。
则 \(P(x) = T + f(x) - f(x - m) - 2g(x)\)。
把字符串当作多项式,即 \(B(x) = b_0 + b_1x + b_2x^2 + \cdots\)。
则 \(g\) 为 \(S\) 与 \(B\) 的卷积。
最后检验多项式 \(P\) 有多少系数为 \(0\) 的项。
带通配符的单模式串匹配
问题概述:给定模式串 \(a\)(长度为 \(m\)),文本串 \(b\)(长度为 \(n\)),串的某些字符是 “通配符”(可以与任意字符匹配)。求所有位置 \(x\),满足 \(b\) 串以 \(x\) 结尾的连续 \(m\) 个字符与 \(a\) 串完全相同。
总结思路:
- 定义匹配函数。
- 定义完全匹配函数。
- 快速计算每一位的完全匹配函数(翻转模式串化为卷积)。
设通配符的数值为 \(0\)。
定义匹配函数 \(C(x, y) = [a_x - b_y]^2a_xb_y\)。
则完全匹配函数 \(P(x) = \sum\limits_{i = 0}^{m - 1}[a_i - b_{(x -m + i + 1)}]^2a_ib_{(x - m + i + 1)}\)。
将 \(a\) 翻转,\(P(x) = \sum\limits_{i = 0}^{m - 1}[s_{(m - i + 1)} - b_{(x -m + i + 1)}]^2\cdot s_ib_{(x - m + i + 1)}\)。
暴力展开:
把右边写成多项式乘积的形式,则
例题
P4173 残缺的字符串
模板题。
#include<bits/stdc++.h>
using namespace std;
using namespace numbers;
using ll = long long;
constexpr int N = 6e5 + 5, tot = 1 << 19;
struct Complex {
double x, y;
Complex operator + (const Complex &o) const {
return {x + o.x, y + o.y};
}
Complex operator - (const Complex &o) const {
return {x - o.x, y - o.y};
}
Complex operator * (const Complex &o) const {
return {x * o.x - y * o.y, x * o.y + y * o.x};
}
} a[N], b[N], c[N];
int rev[N];
void fft(Complex a[], int Inv) {
for(int i = 0; i < tot; ++ i) {
if(i < rev[i]) {
swap(a[i], a[rev[i]]);
}
}
for(int mid = 1; mid < tot; mid <<= 1) {
auto w1 = Complex({cos(pi / mid), Inv * sin(pi / mid)});
for(int i = 0; i < tot; i += mid * 2) {
auto wk = Complex({1, 0});
for(int j = 0; j < mid; ++ j, wk = wk * w1) {
auto x = a[i + j], y = a[i + j + mid];
a[i + j] = x + wk * y;
a[i + j + mid] = x - wk * y;
}
}
}
if(Inv == -1) {
for(int i = 0; i < tot; ++ i) {
a[i].x /= tot;
}
}
}
string s, t;
int m, n, A[N], B[N];
int main() {
cin.tie(0)->sync_with_stdio(0);
cin >> m >> n >> s >> t;
reverse(s.begin(), s.end());
for(int i = 0; i < m; ++ i) {
A[i] = s[i] == '*' ? 0 : s[i] - 'a' + 1;
}
for(int i = 0; i < n; ++ i) {
B[i] = t[i] == '*' ? 0 : t[i] - 'a' + 1;
}
for(int i = 0; i < tot; ++ i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << 18);
}
for(int i = 0; i < tot; ++ i) {
a[i] = Complex({A[i] * A[i] * A[i], 0});
b[i] = Complex({B[i], 0});
}
fft(a, 1), fft(b, 1);
for(int i = 0; i < tot; ++ i) {
c[i] = c[i] + a[i] * b[i];
}
for(int i = 0; i < tot; ++ i) {
a[i] = Complex({A[i], 0});
b[i] = Complex({B[i] * B[i] * B[i], 0});
}
fft(a, 1), fft(b, 1);
for(int i = 0; i < tot; ++ i) {
c[i] = c[i] + a[i] * b[i];
}
for(int i = 0; i < tot; ++ i) {
a[i] = Complex({A[i] * A[i], 0});
b[i] = Complex({B[i] * B[i], 0});
}
fft(a, 1), fft(b, 1);
for(int i = 0; i < tot; ++ i) {
c[i] = c[i] - Complex({2, 0}) * a[i] * b[i];
}
fft(c, -1);
vector<int> ans;
for(int i = m - 1; i < n; ++ i) {
ll v = c[i].x + 0.5;
if(v == 0) {
ans.push_back(i - m + 2);
}
}
cout << ans.size() << '\n';
for(int x : ans) cout << x << ' ';
return 0;
}
CF528D
题意:给定只含 ATGC
的两个字符串 \(s, \ t\) 和一个非负整数 \(k\),求有多少 \(i\) 满足任意 \(j \in [0, |t|)\),都存在 \(p \in [0, |s|)\) 使得 \(|i + j - p|\le k\) 且 \(s_p = t_j\)。
记 \(s\) 的长度为 \(n\),\(t\) 的长度为 \(m\)。
预处理第 \(i\) 个位置能否与字符 \(c\) 匹配,记做 \(b_{i, c}\)。
定义完全匹配函数 \(P(x)\),表示以 \(x\) 结尾的字符串与 \(t\) 匹配的个数,匹配成功当且仅当 \(P(x) = m\)。
\(P(x) = \sum\limits_{c\in t}\sum\limits_{t_i = c}b_{(x - m + i + 1, c)}\)。
单独考虑每个字符,记 \(f(x, A) = \sum_{i = 0}^{m - 1}[t_i = A][b_{x - m + i + 1, A}]\)。
按照套路将 \(t\) 翻转。
则 \(f(x, A) = \sum_{i = 0}^{m - 1}[t'_{m -i - 1} = A][b_{x - m + i + 1, A}]\)。
令 \(T(x) = \sum_{i \ge 0} [t'_i = A]x^i, \ B(x) = \sum_{i \ge 0}[b_{i, A}]x^i\),则 \(f(i, A) = [x^i]T(x)B(x)\)
拉格朗日插值
拉格朗日插值定理
\(n\) 个点值 \((x_i, y_i), \ (1\le i \le n)\),满足 \(x_i \neq x_j, \ (i\neq j)\),它们 唯一 确定一个 \(n - 1\) 次多项式 \(f(x)\) 满足
构造一个多项式 \(f(x) = f_1(x) + f_2(x) + \cdots + f_n(x)\)。
满足 \(f_i(x_j) = \begin{cases}y_j & i = j\\0 & i \neq j\end{cases}\)。
待定系数,\(f_1(x) = A(x - x_2)(x - x_3)\cdots(x - x_n)\)。
带入 \(x_1\),解出 \(A = \dfrac{y_1}{(x_1 - x_2)(x_1 - x_3)\cdots(x_1 - x_n)}\),其余 \(f_i\) 同理。
-
与 \(\text{CRT}\) 的关系?
多项式模多项式:如果 \(f(x) = q(x)g(x) + r(x)\),其中 \(r(x)\) 的次数小于 \(g(x)\),则记做 \(f(x)\equiv r(x)\pmod {g(x)}\)。
多项式互质:两个多项式的公因式只有常数,则称这两个多项式互质。
对于任意 \(i \in [1, n]\):
\[f(x) - f(x_i) =a_1(x - x_i) + a_2(x - x_i^2) + \cdots a_{n - 1}(x^{n - 1} - x_i^{n - 1}) \equiv 0 \pmod{x - x_i} \]有同余方程组
\[\begin{cases} f(x) \equiv f(x_1) \pmod {x - x_1}\\ \\ f(x) \equiv f(x_2) \pmod {x - x_2}\\ \\ \vdots\\ \\ f(x) \equiv f(x_{n - 1}) \pmod {x - x_{n - 1}}\\ \end{cases} \]显然有模数两两互质,满足中国剩余定理的使用条件。
-
暴力实现:先算 \(g(x) = \prod(x - x_i)\),再求 \(n\) 次 \(g(x) / (x - x_i)\),复杂度 \(O(n^2)\)。
-
特殊情况1:只要求一个 \(f(x_0)\),复杂度 \(O(n^2)\)。
-
特殊情况2:只要求一个 \(f(x_0)\),满足 \(x_i = i\),复杂度 \(O(n)\)(可以预处理)。
例题
CF622F
题意:求 \(\sum_{i = 1}^ni^k\pmod {10^9 + 7}, \ n \le 10^9, \ k \le 10^6\)。
\(f(n) = \sum_{i = 1}^ni^k\) 是关于 \(n\) 的 \(k + 1\) 次多项式。
证明:
对于第二类斯特林数,满足公式
那么
有组合恒等式
则
所以求出 \([0, k + 1]\) 的函数值,直接拉插即可。
Polynomia
题意:\(f(n)\) 是 \(n\) 次多项式,给定 \(f(i), \ i \in[0, n]\),询问 \(s(R) - s(L - 1)\)。
\(f(n)\) 是 \(n\) 次多项式,则其前缀和为 \(n + 1\) 次多项式。
先对 \(f(n)\) 做一遍拉插求出 \(f(n + 1)\),因而得到 \(s(n + 1)\),再对 \(s\) 做拉插。