学习笔记：FFT与拉格朗日插值

多项式的表示形式#

系数表示与点值表示#

假设 $f(x)$ 是一个 $n$ 次多项式，则 $f(x)$ 的系数表示为 $f(x)=a_n{x}^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))$，其中 $\mathrm{\forall }i\ne j,x_i\ne x_j$。

• $n+1$ 个点值可以表示一个 $n$ 次多项式。

• 在点值表示下 $n$ 次多项式的乘法复杂的为 $O(n)$。

  设 $h(x)=f(x)\cdot g(x)$，则 $\mathrm{\forall }i\in [0,n],h(x_i)=f(x_i)\cdot g(x_i)$。

复数与单位根#

复数的指数形式#

$a+bi=r{e}^{i\theta }$，其中 $r=\sqrt{a^2+b^2},\tan \theta ={\displaystyle \frac{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_{1}x+\cdots +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_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+(a_1+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1})$。

令 $B(x)=a_0+a_{2}x+\cdots +a_{n-2}{x}^{n/2-1},C(x)=a_1+a_{3}x+\cdots +a_{n-2}{x}^{n/2-1}$。

则 $A(x)=B(x^2)+xC(x^2)$。

对于 $k\in [0,n/2-1]$

$$\begin{array}{rl}A({\omega }_n^k)& =B({\omega }_n^{2k})+{\omega }_n^{k}C({\omega }_n^{2k})\\ \\ & =B({\omega }_{n/2}^k)+{\omega }_n^{k}C({\omega }_{n/2}^k)\end{array}$$

对于另一半的点值，可用 $A({\omega }_n^{k+n/2})$ 表示，即

$$\begin{array}{rl}A({\omega }_n^{k+n/2})& =B({\omega }_n^{2k+n})+{\omega }_n^{k+n/2}C({\omega }_n^{2k+n})\\ \\ & =B({\omega }_{n/2}^k)-{\omega }_n^{k}C({\omega }_{n/2}^k)\end{array}$$

具体的讲，当 $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_{1}x+\cdots +a_{n-1}{x}^{n-1}$。

设 $n\times n$ 的矩阵 $\mathrm{\Omega }$，其中 ${\mathrm{\Omega }}_{i,j}={\omega }_n^{ij}$，设向量 $a=(a_0,a_1,\dots ,a_{n-1}),b=(b_0{,}_1,\dots ,b_{n-1})$。

则 IDFT 相当于求解方程 $\mathrm{\Omega }a=b$。

$$\left(\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right)=\left(\begin{array}{c}{b}_0\\ b_1\\ ⋮\\ b_{n-1}\end{array}\right)$$

如果我们能够求出 $\mathrm{\Omega }$ 的逆，则 $a={\mathrm{\Omega }}^{-1}b$。

设

$$M=\bar{\mathrm{\Omega }}\cdot \mathrm{\Omega }=\left(\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& \cdots & ({\omega }_n^{-0}{)}^{n-1}\\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& \cdots & ({\omega }_n^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(n-1)}{)}^0& ({\omega }_n^{-(n-1)}{)}^1& \cdots & ({\omega }_n^{-(n-1)}{)}^{n-1}\end{array}\right)\left(\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right)$$

其中

$$\begin{array}{rl}{M}_{i,j}& =({\omega }_n^{-i}{)}^0\cdot ({\omega }_n^0{)}^j+({\omega }_n^{-i}{)}^1\cdot ({\omega }_n^1{)}^j+\cdots ({\omega }_n^{-i}{)}^{n-1}\cdot ({\omega }_n^{n-1}{)}^j\\ \\ & =({\omega }_n^{j-i}{)}^0+({\omega }_n^{j-i}{)}^1+\cdots ({\omega }_n^{j-i}{)}^{n-1}\end{array}$$

相当于等比数列求和

$$M_{i,j}\{\begin{array}{ll}{\displaystyle \frac{1-(w_n^{j-i}{)}^n}{1-w_n^{j-i}}}=0& w_n^{j-i}\ne 1\text{ 即 }i\ne j\\ \\ n& i=j\end{array}$$

令 $\bar{\mathrm{\Omega }}$ 满足 ${\bar{\mathrm{\Omega }}}_{i,j}={\omega }^{-ij}$，有 $\bar{\mathrm{\Omega }}\cdot \mathrm{\Omega }=nI$，因此

$$a={\displaystyle \frac{1}{n}}\bar{\mathrm{\Omega }}b$$

相当于给定 $B(x)=b_0+b_{1}x+\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_1{10}^1+a_2{10}^2+\cdots$，fft 后再处理进位。

submission

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 _i\sum _j\sum _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$。

1. 减去满足三个性质中一个的方案。

   钦定两个位置相等。

   令 $B(x)=\sum x^{2s_i}$，产生的方案为 $[x^V](A(x)B(x))$。

2. 加上满足三个性质中两个的方案。

   $i=j=k$。

   令 $C(x)=\sum x^{3s_i}$，这部分的方案为 $[x^V]C(x)$。

3. 减去满足三个性质中三个的方案。

   与第二部分相同。

则最终答案 $D(x)=A^3(x)-3B(x)+2C(x)$。

可以先对 $A(x),B(x),C(x)$ 分别做 dft，计算出 $D(x)$ 的点值，最后再对 $D(x)$ 做 idft。

submission

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 _{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 _{i=0}^{m-1}[s_{(m-i+1)}-b_{(x-m+i+1)}{]}^2$。

暴力展开：$P(x)=\sum _{i=0}^{m-1}{s}_{(m-i+1)}^2+\sum _{i=0}^{m-1}{b}_{(x-m+i+1)}^2-\sum _{i=0}^{m-1}2\cdot s_{(m-i+1)}{b}_{(x-m+i+1)}$。

第一项为常数。

第二项可以预处理 $f(x)=\sum _{i=0}^x{b}_i^2$，计算 $f(x)-f(x-m)$。

第三项等效于求 $g(x)=\sum _{i+j=x}{s}_i{b}_j$。

则 $P(x)=T+f(x)-f(x-m)-2g(x)$。

把字符串当作多项式，即 $B(x)=b_0+b_{1}x+b_2{x}^2+\cdots$。

则 $g$ 为 $S$ 与 $B$ 的卷积。

最后检验多项式 $P$ 有多少系数为 $0$ 的项。

带通配符的单模式串匹配#

问题概述：给定模式串 $a$（长度为 $m$），文本串 $b$（长度为 $n$），串的某些字符是 “通配符”（可以与任意字符匹配）。求所有位置 $x$，满足 $b$ 串以 $x$ 结尾的连续 $m$ 个字符与 $a$ 串完全相同。

总结思路：

1. 定义匹配函数。
2. 定义完全匹配函数。
3. 快速计算每一位的完全匹配函数（翻转模式串化为卷积）。

设通配符的数值为 $0$。

定义匹配函数 $C(x,y)=[a_x-b_y{]}^2{a}_x{b}_y$。

则完全匹配函数 $P(x)=\sum _{i=0}^{m-1}[a_i-b_{(x-m+i+1)}{]}^2{a}_i{b}_{(x-m+i+1)}$。

将 $a$ 翻转，$P(x)=\sum _{i=0}^{m-1}[s_{(m-i+1)}-b_{(x-m+i+1)}{]}^2\cdot s_i{b}_{(x-m+i+1)}$。

暴力展开：

$$\begin{array}{rl}P(x)& =\sum _{i=0}^{m-1}[s_{(m-i+1)}-b_{(x-m+i+1)}{]}^2\cdot s_{(m-i+1)}{b}_{(x-m+i+1)}\\ \\ & =\sum _{i=0}^{m-1}{s}_{(m-i+1)}^3{b}_{(x-m+i+1)}+\sum _{i=0}^{m-1}{s}_{(m-i+1)}{b}_{(x-m+i+1)}^3-\sum _{i=0}^{m-1}2\cdot s_{(m-i+1)}^2{b}_{(x-m+i+1)}^2\\ \\ & =\sum _{i+j=x}{s}_i^3{b}_j+\sum _{i+j=x}{s}_i{b}_j^3-2\cdot \sum _{i+j=x}{s}_i^2{b}_j^2\end{array}$$

把右边写成多项式乘积的形式，则

$$\begin{array}{r}\\ P(x)=\sum s_i^3{x}^i\cdot \sum b_i{x}^i+\sum s_i{x}^i\cdot \sum b_i^3{x}^i-2\cdot \sum s_i^2{x}^i\cdot \sum b_i^2{x}^i\\ \end{array}$$

例题#

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 _{c\in t}\sum _{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}^{\prime }=A][b_{x-m+i+1,A}]$。

令 $T(x)=\sum _{i\ge 0}[t_i^{\prime }=A]x^i,B(x)=\sum _{i\ge 0}[b_{i,A}]x^i$，则 $f(i,A)=[x^i]T(x)B(x)$

submission

拉格朗日插值#

拉格朗日插值定理#

$n$ 个点值 $(x_i,y_i),(1\le i\le n)$，满足 $x_i\ne x_j,(i\ne j)$，它们 唯一 确定一个 $n-1$ 次多项式 $f(x)$ 满足

$$f(x)=\sum _{i=1}^n{y}_i\prod _{j\ne i}{\displaystyle \frac{x-x_j}{x_i-x_j}}$$

构造一个多项式 $f(x)=f_1(x)+f_2(x)+\cdots +f_n(x)$。

满足 $f_i(x_j)=\{\begin{array}{ll}{y}_j& i=j\\ 0& i\ne j\end{array}$。

待定系数，$f_1(x)=A(x-x_2)(x-x_3)\cdots (x-x_n)$。

带入 $x_1$，解出 $A={\displaystyle \frac{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)(\mathrm{mod}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(\mathrm{mod}x-x_i)$$

  有同余方程组

  $$\{\begin{array}{l}f(x)\equiv f(x_1)(\mathrm{mod}x-x_1)\\ \\ f(x)\equiv f(x_2)(\mathrm{mod}x-x_2)\\ \\ ⋮\\ \\ f(x)\equiv f(x_{n-1})(\mathrm{mod}x-x_{n-1})\end{array}$$

  显然有模数两两互质，满足中国剩余定理的使用条件。

• 暴力实现：先算 $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}^n{i}^k(\mathrm{mod}{10}^9+7),n\le {10}^9,k\le {10}^6$。

$f(n)=\sum _{i=1}^n{i}^k$ 是关于 $n$ 的 $k+1$ 次多项式。

证明：

对于第二类斯特林数，满足公式

$$n^m=\sum _{i=0}^m\left\{\begin{array}{c}m\\ i\end{array}\right\}{n}^{\underset{―}{i}}$$

那么

$$f(n)=\sum _{i=1}^n\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}{i}^{\underset{―}{j}}=\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}j!\cdot \sum _{i=1}^n\left(\begin{array}{c}i\\ j\end{array}\right)$$

有组合恒等式

$$\sum _{l=0}^n\left(\begin{array}{c}l\\ k\end{array}\right)=\left(\begin{array}{c}n+1\\ k+1\end{array}\right)$$

则

$$f(n)=\sum _{i=0}^k\left\{\begin{array}{c}k\\ i\end{array}\right\}i!\left(\begin{array}{c}n+1\\ i+1\end{array}\right)=\sum _{i=0}^k\left\{\begin{array}{c}k\\ i\end{array}\right\}{\displaystyle \frac{(n+1{)}^{\underset{―}{i+1}}}{i+1}}$$

所以求出 $[0,k+1]$ 的函数值，直接拉插即可。

submission

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$ 做拉插。

submission

reference#

Ebola:浅谈FFT在字符串匹配中的应用