目錄
寫在前面
- 這是繼數論和組合計數類數學相關與多項式類數學相關后的第三篇數學方面內容總結。主要記錄自己近期學習的一些數學方法。內容比較雜,同時也起到對前文的一些補充作用。
- 若文章中出現錯誤,煩請告知。
- 感謝您的造訪。
一類反演問題
莫比烏斯反演
這一個內容在數論和組合計數類數學相關里面提到了。
快速莫比烏斯變換(反演)與子集卷積
莫比烏斯變換(反演)
問題提出
若 \(h,f,g\) 為下標為集合的函數,我們定義
\[h=f*g\]
表示
\[h(S) = \sum_{L \subseteq S}^{} \sum_{R \subseteq S}^{} [L \cup R = S] f(L) \times g(R)\]
容易發現,對於這個問題,我們可以用 \(O\left((2^n)^2\right)\) 的枚舉 \(L,R\) 來計算。
然而這樣復雜度較高,我們考慮類比多項式卷積的過程,可以求出 \(f,g\) 的點值,直接相乘得到 \(h\) 的點值然后再插回去。
值得注意的是為了便於表述以及規范表達,快速莫比烏斯變換就相當於點值,快速莫比烏斯反演就相當於插值。
算法原理
- 我們定義 \(f\) 的莫比烏斯變換為 \(F\) ,其中 \(F(S)=\sum_{X\subseteq S}f(X)\) ;由這個定義,我們可以推出 \(F\) 莫比烏斯反演 \(f\) 為 \(f(S) = \sum_{X \subseteq S} (-1)^{|S| - |X|}F(X)\) 。對於莫比烏斯反演的證明,可以帶入莫比烏斯變換的式子或容斥來證。
- 我們對於一個函數 \(f,g,h\) ,記它的點值式為 \(F,G,H\) 。我們將問題提出中的卷積式兩邊同時做莫比烏斯變換,得到 \[\begin{aligned}H(S) &= \sum_{L \subseteq S} \sum_{R \subseteq S} [L \cup R \subseteq S] f(L) \times g(R)\\&=\sum_{L \subseteq S} \sum_{R \subseteq S}f(L) \times g(R)\\&=\left(\sum_{L \subseteq S}f(L)\right)\times\left(\sum_{R \subseteq S}g(R)\right)\\&=F(S)\times G(S)\end{aligned}\]
至此算法原理及過程已經完全結束。似乎我們可以用 \(O(3^n)\) 枚舉子集來變換和反演,實際上我們可以讓復雜度更優。
算法實現
- 設 \(\hat{f_{S}}^{(i)}\) 表示 \(\sum_{T \subseteq S} [(S - T) \subseteq \{1, 2, ..., i\}] f_{T}\)
- 易得初始狀態:\(\hat{f_{S}}^{(0)} = f_{S}\)
- 對於每一個不包含 \(\{i\}\) 的集合 \(S\),可知 \(\hat{f_{S}}^{(i)} = \hat{f_{S}}^{(i - 1)}\)(因為 \(S\) 並沒有 \(i\) 這位),\(\hat{f}_{S \cup \{i\}}^{(i)} = \hat{f}_{S}^{(i - 1)} + \hat{f}_{S \cup \{i\}}^{(i - 1)}\)(前者的 \(T\) 沒有包含 \(\{i\}\),而后者的 \(T\) 必須包含了 \(\{i\}\))。
- 顯然,遞推了 \(n\) 輪之后,\(\hat{f}_{S}^n\) 就是所求的變換了。
用高維前綴和可以做到 \(O(n\times 2^n)\) 的遞推,求出點值和插值。
void FMT(int *A, int o) {// o 為識別因子
for (int i = 1; i < ST; i <<= 1)//ST-1 表示全集
for (int j = 0; j < ST; j++)
if (i&j) (A[j] += A[j^i]*o) %= mod;
}
例題
子集卷積
\(\text{FWT}\) :“你剛才說的那個玩意我也能做啊,要你何用?”
\(\text{FMT}\) :“……”
問題提出
若 \(h,f,g\) 為下標為集合的函數,我們定義
\[h=f*g\]
表示
\[h(S) = \sum_{X \subseteq S} f(X) \times g(S-X)\]
算法實現
回顧剛剛的集合並卷積,子集卷積的條件比集合並卷積更苛刻,即 \(L\) 和 \(R\) 的集合應該不相交。
我們可以在卷積時多加一維,維護集合的大小,如 \(f_{i,S}\) 表示集合中有 \(i\) 個元素,集合表示為 \(S\) 。顯然,當 \(i\) 和 \(S\) 的真實元素個數符合時才是對的。記數組 cnt[S]
表示集合 \(S\) 的模。初始時,我們只把 \(f_{cnt[S],S}\) 的值賦成原來的 \(f(S)\) ( \(g\) 同理),然后每一維做一遍 \(\text{FMT}\) ,點值相乘時這么寫:\(h_{i, S} = \sum_{j = 0}^{i} f_{j,S} \times g_{i - j, S}\) 。最后掃一遍把不符合實際情況的狀態賦成 \(0\) 即可。
for (int i = 0; i <= n; i++) FMT(g[i], 1);
for (int i = 0; i <= n; i++) FMT(f[i], 1);
for (int i = 1; i <= n; i++) {
for (int j = 0; j <= i; j++)
for (int k = 0; k < ST; k++)
(h[i][k] += 1ll*f[j][k]*g[i-j][k]%mod) %= mod;
FMT(h[i], -1);
for (int k = 0; k < ST; k++) if (cnt[k] != i) h[i][k] = 0;
if (i != n) FMT(h[i], 1);
}
例題
二項式反演
內容
對於函數 \(f,g\) , \(\forall p\in\mathbb{N}\) 若 \(\forall n\geq p\) ,滿足
\[f(n)=\sum_{k=p}^{n}\binom{n}{k}g(k)\]
那么 \(\forall n\geq p\)
\[g(n)=\sum_{k=p}^{n}(-1)^{n-k}\binom{n}{k}f(k)\]
證明
為了方便表達,我們取 \(p=0\) ,實質和取 \(p\in\mathbb{N}\) 的證明方法是一樣的。
\[\begin{aligned}g(n)&=\sum_{k=0}^{n}(-1)^{n-k}\binom{n}{k}f(k)\\&=\sum_{k=0}^{n}(-1)^{n-k}\binom{n}{k}\sum_{i=0}^k{k\choose i}g(i)\\&=\sum_{k=0}^n\sum_{i=0}^k(-1)^{n-k}{n\choose k}{k\choose i}g(i)\\&=\sum_{i=0}^n\left(\sum_{k=i}^n(-1)^{n-k}{n\choose k}{k\choose i}\right)g(i)\\&=\sum_{i=0}^n\left(\sum_{k=i}^n(-1)^{n-k}{n\choose i}{n-i\choose k-i}\right)g(i)\\&=\sum_{i=0}^n\left({n\choose i}\sum_{k=i}^n(-1)^{n-k}{n-i\choose n-k}\right)g(i)\\&=\sum_{i=0}^n\left({n\choose i}(1-1)^{n-i}\right)g(i)\\&=g(n)\end{aligned}\]
故成立。
應用舉例
推導錯排公式
我們記 \(f(n)\) 為 \(n\) 個數字任意放的方案數, \(g(n)\) 為 \(n\) 個數沒有一個放在自己位置上的方案數。
枚舉不在自己位置上的個數,容易得到
\[f(n)=\sum_{i=0}^n{n\choose i}g(i)\]
那么
\[\begin{aligned}g(n)&=\sum_{i=0}^n(-1)^{n-i}{n\choose i}f(i)\\&=\sum_{i=0}^n(-1)^i{n\choose i}f(n-i)\end{aligned}\]
注意到 \(f(x)=x!\) ,那么
\[\begin{aligned}g(n)&=\sum_{i=0}^n(-1)^i\frac{n!}{i!(n-i)!}(n-i)!\\&=n!\sum_{i=0}^n(-1)^i\frac{1}{i!}\end{aligned}\]
棋盤染色
有個 \(1\times n\) 的格子, \(m\) 種顏色( \(m\geq 2\) ),要求相鄰格子的顏色不相同且每種顏色都要用到,求染色方案數。
我們記 \(f(n)\) 為至多用到 \(n\) 種顏色的方案數, \(g(n)\) 為 \(n\) 為恰用到 \(n\) 種顏色的方案數。
那么
\[\begin{aligned}f(m)&=\sum_{i=2}^m{m\choose i}g(i)\\\Rightarrow g(m)&=\sum_{i=2}^m(-1)^i{m\choose i}f(n-i)\end{aligned}\]
注意到 \(f(x)=x\times(x-1)^{n-1}\) 。那么就可以帶入直接算了。
另一形式
\[a_k=\sum\limits_{i=k}^n{i\choose k}b_i\Rightarrow b_k=\sum\limits_{i=k}^n(-1)^{i-k}{i\choose k}a_i\]
證明:
\[\begin{aligned} &\sum\limits_{i=k}^n(-1)^{i-k}{i\choose k}a_i\\=& \sum\limits_{i=k}^n(-1)^{i-k}{i\choose k}\sum\limits_{j=k}^n{j\choose i}b_i\\=& \sum\limits_{i=k}^n\sum\limits_{j=k}^n(-1)^{i-k}{i\choose k}{j\choose i}b_i\\=& \sum\limits_{i=k}^n\sum\limits_{j=k}^n(-1)^{i-k}{j\choose i}{i\choose k}b_i\\=& \sum\limits_{j=k}^n\sum\limits_{i=k}^j(-1)^{i-k}{j\choose k}{j-k\choose i-k}b_i\\=& \sum\limits_{j=k}^n{j\choose k}\sum\limits_{i=k}^j(-1)^{i-k}{j-k\choose i-k}b_i\\=&\sum\limits_{j=k}^n{j\choose k}(1-1)^{j-k}b_i\\=&b_k\end{aligned}\]
例題
斯特林反演
第一類斯特林數
\(\begin{bmatrix}n\\m\end{bmatrix}\) 表示將 \(n\) 個元素排成 \(m\) 個輪換的方法數。
遞推公式:\(\begin{bmatrix}n\\m\end{bmatrix}=\begin{bmatrix}n-1\\m-1\end{bmatrix}+(n-1)\begin{bmatrix}n-1\\m\end{bmatrix}\)
含義是考慮第 \(n\) 個元素的放法:要么新開一個輪換,要么就放在前 \(n-1\) 個元素的左邊。
第二類斯特林數
\(\begin{Bmatrix}n\\m\end{Bmatrix}\) 表示將 \(n\) 個元素划分成 \(m\) 個非空子集的方法數。
遞推公式:\(\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix}\)
含義是考慮第 \(n\) 個元素的放法:要么新開一個組,要么就放在前 \(m\) 組內。
通項公式(容斥式): \(\begin{Bmatrix}n\\m\end{Bmatrix}=\frac{1}{m!}\sum\limits_{k=0}^{m}(-1)^k\binom{m}{k}(m-k)^n\)
有關通項公式的證明及運用可以參考多項式類數學相關這篇文章。
例題
反演公式
\[f(x) = \sum_{i=0}^x \begin{Bmatrix}x\\i\end{Bmatrix} g(i) \Leftrightarrow g(x) = \sum_{i=0}^x (-1) ^ {x - i}\begin{bmatrix}x\\i\end{bmatrix} f(i)\]
\[f(x) = \sum_{i=0}^x \begin{bmatrix}x\\i\end{bmatrix} g(i) \Leftrightarrow g(x) = \sum_{i=0}^x (-1) ^ {x - i}\begin{Bmatrix}x\\i\end{Bmatrix} f(i)\]
例題
- 給出 \(n\) 個點的一張簡單圖,問有多少個邊的子集,滿足保留子集中的邊后,該圖連通。(蒯自Sdchr)
- 大概就是枚舉連通塊的個數,然后塊內隨便連,然后容斥就好。
- 考慮如何求容斥系數 \(f(i)\) 。設實際上是 \(x\) 個連通塊的方案,它應該被計算 \([x=1]\) 次,實際上在所有更仔細的分塊中被統計,所以
- \[[x = 1] = \sum_{i=1}^x\begin{Bmatrix}x\\i\end{Bmatrix}f(i)\]
- 由斯特林反演
- \[\begin{aligned}f(x) &= \sum_{i=1}^x(-1)^{x-i}\begin{bmatrix}x\\i\end{bmatrix}[i=1]\\&=(-1)^{x-1}\begin{bmatrix}x\\1\end{bmatrix}\\&=(-1)^{x-1}(x-1)!\end{aligned}\]
- [BZOJ 4671]異或圖
最值反演( \(\text{min-max}\) 容斥)
公式
記 \(\max(S)\) 為集合 \(S\) 中的最大值, \(\min(S)\) 為集合 \(S\) 中的最小值, \(|S|\) 為集合 \(S\) 的元素數量,那么以下兩個等式成立
\[\max(S)=\sum_{T\subseteq S}(-1)^{|T|+1}\min(T)\]
\[\min(S)=\sum_{T\subseteq S}(-1)^{|T|+1}\max(T)\]
證明
這里只證明第一個等式好了,后邊的可以自行推出。
其實只需要證明一件事,就是除了 \(\min(T)=\max(S)\) 的那個值,其他的 \(\min\) 值都被消掉了就可以了(這里說明一下,我們假定集合中的元素兩兩相異)
先來說明 \(\max(S)\) 的系數為什么是 \(1\) ,假設中 \(S\) 最大的元素是 \(a\) ,那么我們會發現只有 \(\min(\{a\})=\max(S)\) 所以 \(\max(S)\) 的系數必須是 \(1\) 。
然后再說明為什么別的 \(\min\) 都被消掉了,假設某個元素 \(b\) 的排名是 \(k\) ,那么 \(\min(T)=b\) 當且僅當我們選出的集合是后 \(n-k\) 個的元素構成的集合的子集然后並上 \(\{b\}\) 得到的,我們會發現顯然這樣的集合有 \(2^{n-k}\) 種,而顯然這其中恰有 \(2^{n-k-1}\) 中是有奇數個元素的,恰有 \(2^{n-k-1}\) 種是有偶數個元素的,兩兩相消自然就成 \(0\) 了,當然上述等式在 \(k=n\) 的時候不成立,但是此時剩下的剛好是最大值,所以證明完畢。
拉格朗日插值法
簡介
給定 \(n + 1\) 個橫坐標不相同的點,可以唯一確定一個 \(n\) 次的多項式。最直觀的求多項式的做法就是列方程求解。但是這樣需要 \(O(n^3)\) 的時間來計算。而拉格朗日插值法則通過構造的方法,得到了一個經過 \(n + 1\) 個點的 \(n\) 次多項式。 具體的過程是這樣的,假設現在我們得到了 \(n + 1\) 個點:
\[ (x_0, y_0), \;(x_1, y_1),\; \dots,\;(x_n, y_n) \]
設拉格朗日基本多項式為
\[ \ell_j(x) = \prod_{i = 0, i \neq j}^n {x - x_i \over x_j - x_i} \]
這個基本多項式構造十分巧妙,因為注意到 \(\ell_j(x_j) = 1\) ,並且 \(\ell_j(x_i) = 0, \;\forall i \neq j\) 。那么,接着構造出這個 \(n\) 次多項式
\[ P(x) = \sum_{i = 0}^n y_i\ell_i(x) \]
根據基本多項式的性質,我們可以知道 \(P(x_i) = y_i\) ,也就是經過了這 \(n + 1\) 個點。 通過簡單的多項式乘法和多項式除法就可以在 \(O(n^2)\) 的時間求出這個多項式的系數表達。
求解
已知 \(n\) 次多項式 \(f(n)\) 上的 \(n+1\) 個點 \((x_i,y_i),i\in[0,n]\) ,求 \(f(xi)\)
int lagrange(int n, int *x, int *y, int xi) {
int ans = 0;
for (int i = 0; i <= n; i++) {
int s1 = 1, s2 = 1;
for (int j = 0; j <= n; j++)
if (i != j) {
s1 = 1ll*s1*(xi-x[j])%mod;
s2 = 1ll*s2*(x[i]-x[j])%mod;
}
ans = (1ll*ans+1ll*y[i]*s1%mod*quick_pow(s2, mod-2)%mod)%mod;
}
return (ans+mod)%mod;
}
如果 \(x\) 的取值是連續一段的話,我們可以做到 \(O(n)\) 求解。假設 \(\forall i<j,x_i<x_j\) (具體公式推導的話,如果你有興趣可以參看之后的內容。因為比較顯然,這里不再講解。)
int lagrange(int n, int *x, int *y, int xi) {
int ans = 0;
s1[0] = (xi-x[0])%mod, s2[n+1] = 1;
for (int i = 1; i <= n; i++) s1[i] = 1ll*s1[i-1]*(xi-x[i])%mod;
for (int i = n; i >= 0; i--) s2[i] = 1ll*s2[i+1]*(xi-x[i])%mod;
ifac[0] = ifac[1] = 1;
for (int i = 2; i <= n; i++) ifac[i] = -1ll*mod/i*ifac[mod%i]%mod;
for (int i = 2; i <= n; i++) ifac[i] = 1ll*ifac[i]*ifac[i-1]%mod;
for (int i = 0; i <= n; i++)
(ans += 1ll*y[i]*(i == 0 ? 1 : s1[i-1])%mod*s2[i+1]%mod
*ifac[i]%mod*(((n-i)&1) ? -1 : 1)*ifac[n-i]%mod) %= mod;
return (ans+mod)%mod;
}
例題
自然數的冪的前綴和
問題提出
給定的 \(n\) 和 \(k\) ,求
\[\sum_{i = 1}^n i^k \]
通常 \(n\) 比較大,而 \(k\) 只有幾千或者幾萬。
問題解決
我們可以知道,對於上述式子,推導公式一定是是 \(k + 1\) 次多項式。對於證明的話,我們可以參考riteme的介紹。
考慮使用拉格朗日插值法來獲得答案多項式。
首先如果我們得知了 \(n = 0, 1, \dots, k + 1\) 處的答案 \(f(n)\) ,那么給定的 \(n\) 處的答案可以寫成
\[ \begin{aligned} f(n) & = \sum_{i = 0}^{k + 1} f(i) {(n - 0)(n - 1)\cdots[n - (i - 1)][(n - (i + 1)]\cdots[n - (k + 1)] \over (i - 0)(i - 1)\cdots[i - (i - 1)][(i - (i + 1)]\cdots[i - (k + 1)]} \\ & = \sum_{i = 0}^{k + 1} f(i) {\prod_{j = 0}^{i - 1} (n - j) \prod_{j = i + 1}^{k + 1} (n - j) \over i!(-1)^{k - i + 1}(k + 1 - i)!} \\ & = \sum_{i = 0}^{k + 1} (-1)^{k - i + 1}f(i) {\prod_{j = 0}^{i - 1} (n - j) \prod_{j = i + 1}^{k + 1} (n - j) \over i!(k + 1 - i)!} \end{aligned} \]
注意到后面的分式中,分子是一個前綴積乘以一個后綴積,而分母是兩個階乘。這些都可以在 \(O(k)\) 的時間內求出。 現在剩下的問題就是如何求出 \(f(0), f(1), \dots, f(k + 1)\) 了。由於 \(g(x) = x^k\) 是個完全積性函數,所以我們可以通過歐拉篩法求出 \(g\) 函數前面的一些值。具體的就是對於質數采取直接快速冪,合數則拆出任意一個因子來算,通常是歐拉篩法中可以順便求得的最小質因子。根據素數定理,素數大約有 \(O\left(\frac{k}{\ln k}\right)\) 個。每次快速冪需要花費 \(O(\log k)\) 的時間,因此總的時間復雜度可以估計為 \(O(k)\) ,是一個非常優秀的算法。上面的方法具有通用性,只要我們可以快速的求出某個 \(k\) 次多項式的前 \(k + 1\) 個值,那么剩下的部分可以使用拉格朗日插值法在 \(O(k)\) 的時間內完成計算。
代碼實現
int lagrange(int k, int *f, int xi) {//k+2 個點對 (i, f[i]), 0 <= i <= k+1
int ans = 0; ++k;
s1[0] = xi, s2[k+1] = 1;
for (int i = 1; i <= k; i++) s1[i] = 1ll*s1[i-1]*(xi-i)%mod;
for (int i = k; i >= 0; i--) s2[i] = 1ll*s2[i+1]*(xi-i)%mod;
for (int i = 0; i <= k; i++)
(ans += 1ll*f[i]*(i == 0 ? 1 : s1[i-1])%mod*s2[i+1]%mod
*ifac[i]%mod*(((k-i)&1) ? -1 : 1)*ifac[k-i]%mod) %= mod;
return (ans+mod)%mod;
}
void pre() {//預處理出階乘逆元、插值的 k+2 個點
f[1] = ifac[0] = ifac[1] = 1;
for (int i = 2; i <= k+1; i++) ifac[i] = -1ll*mod/i*ifac[mod%i]%mod;
for (int i = 2; i <= k+1; i++) ifac[i] = 1ll*ifac[i-1]*ifac[i]%mod;
memset(isprime, 1, sizeof(isprime));
for (int i = 2; i <= k+1; i++) {
if (isprime[i]) prime[++tot] = i, f[i] = quick_pow(i, k);
for (int j = 1; j <= tot && prime[j]*i <= k+1; j++) {
isprime[i*prime[j]] = 0;
f[i*prime[j]] = 1ll*f[i]*f[prime[j]]%mod;
if (i%prime[j] == 0) break;
}
}
for (int i = 1; i <= k+1; i++) f[i] = (f[i]+f[i-1])%mod;
}
void work() {
scanf("%d", &k); pre();
while (~scanf("%d", &n)) {
if (n <= k+1) printf("%d\n", f[n]);
else printf("%d\n", lagrange(k, f, n));
}
}
例題