學習筆記：FFT與拉格朗日插值

Lu_xZ發表於2024-05-09
筆記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]\)

\[\begin{aligned} 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{aligned} \]

對於另一半的點值，可用 \(A(\omega_n^{k + n / 2})\) 表示，即

\[\begin{aligned} 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{aligned} \]

具體的講，當 \(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\)。

\[\begin{pmatrix} (\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}\\ \vdots & \vdots& \ddots & \vdots \\ (\omega_n^{n- 1})^0 & (\omega_n^{n - 1})^1 & \cdots & (\omega_n^{n - 1})^{n - 1}\\ \end{pmatrix} \begin{pmatrix} a_0\\ a_1\\ \vdots\\ a_{n - 1} \end{pmatrix} = \begin{pmatrix} b_0\\ b_1\\ \vdots\\ b_{n - 1} \end{pmatrix} \]

如果我們能夠求出 \(\Omega\) 的逆，則 \(a = \Omega^{-1}b\)。

設

\[M = \overline\Omega\cdot\Omega = \begin{pmatrix} (\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}\\ \vdots & \vdots& \ddots & \vdots \\ (\omega_n^{-(n- 1)})^0 & (\omega_n^{-(n - 1)})^1 & \cdots & (\omega_n^{-(n- 1)})^{n - 1}\\ \end{pmatrix} \begin{pmatrix} (\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}\\ \vdots & \vdots& \ddots & \vdots \\ (\omega_n^{n- 1})^0 & (\omega_n^{n - 1})^1 & \cdots & (\omega_n^{n - 1})^{n - 1}\\ \end{pmatrix} \]

其中

\[\begin{aligned} 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{aligned} \]

相當於等比數列求和

\[M_{i, j} \begin{cases} \dfrac{1 - (w_n^{j - i})^n}{1 - w_n^{j - i}} = 0 \quad & w_n^{j - i} \neq 1\text{ 即 } i \neq j \\ \\ n \quad & i = j \end{cases} \]

令 \(\overline\Omega\) 滿足 \(\overline\Omega_{i, j} = \omega^{-ij}\)，有 \(\overline\Omega\cdot\Omega = nI\)，因此

\[a = \dfrac{1}{n}\overline\Omega b \]

相當於給定 \(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 後再處理進位。

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\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\)。

  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\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\) 串完全相同。

總結思路：

  1. 定義匹配函式。
  2. 定義完全匹配函式。
  3. 快速計算每一位的完全匹配函式（翻轉模式串化為卷積）。

設萬用字元的數值為 \(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)}\)。

暴力展開：

\[\begin{aligned} P(x) &= \sum\limits_{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\limits_{i = 0}^{m - 1} s_{(m - i + 1)}^3b_{(x - m + i + 1)} + \sum\limits_{i = 0}^{m - 1}s_{(m - i + 1)}b_{(x - m + i + 1)}^3 - \sum\limits_{i = 0}^{m - 1}2\cdot s_{(m - i + 1)}^2b_{(x - m + i + 1)}^2\\ \\ &= \sum_{i + j = x}s_i^3b_j + \sum_{i + j = x}s_ib_j^3 -2 \cdot \sum_{i + j = x}s_i^2b_j^2\\ \end{aligned} \]

把右邊寫成多項式乘積的形式，則

\[\begin{aligned} \\ P(x) = \sum s_i^3x^i \cdot \sum b_ix^i +\sum s_ix^i \cdot \sum b_i^3x^i -2\cdot \sum s_i^2x^i\cdot \sum b_i^2x^i\\ \\ \end{aligned} \]

例題

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)\)

submission

拉格朗日插值

拉格朗日插值定理

\(n\) 個點值 \((x_i, y_i), \ (1\le i \le n)\)，滿足 \(x_i \neq x_j, \ (i\neq j)\)，它們 唯一 確定一個 \(n - 1\) 次多項式 \(f(x)\) 滿足

\[f(x) = \sum_{i = 1}^ny_i\prod_{j \neq i}\dfrac{x - x_j}{x_i - x_j} \]

構造一個多項式 \(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\) 次多項式。

證明：

對於第二類斯特林數，滿足公式

\[n^m = \sum\limits_{i = 0}^m\begin{Bmatrix}m\\i\end{Bmatrix}n^{\underline{i}} \]

那麼

\[f(n) = \sum\limits_{i = 1}^n \sum\limits_{j = 0}^k\begin{Bmatrix}k\\j\end{Bmatrix}i^{\underline{j}} = \sum\limits_{j = 0}^k\begin{Bmatrix}k\\ j\end{Bmatrix}j!\cdot \sum\limits_{i = 1}^n\begin{pmatrix}i\\ j\end{pmatrix} \]

有組合恆等式

\[\sum\limits_{l = 0}^n\begin{pmatrix}l\\ k\end{pmatrix} =\begin{pmatrix}n + 1\\ k + 1\end{pmatrix} \]

則

\[f(n)= \sum\limits_{i = 0}^k\begin{Bmatrix}k\\ i\end{Bmatrix}i!\begin{pmatrix}n + 1\\ i + 1\end{pmatrix}= \sum\limits_{i = 0}^k\begin{Bmatrix}k\\ i\end{Bmatrix}\dfrac{(n + 1)^{\underline{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在字串匹配中的應用

相關文章