部落格園
我的部落格
快速沃爾什變換解決的卷積問題
快速沃爾什變換(FWT)是解決這樣一類卷積問題:
其中,\(\odot\) 是位運算的一種。舉個例子,給定數列 \(a,b\),求:
FWT 的思想
看到 FWT 的名字,我們可以聯想到之前學過的 FFT(很可惜,我沒有寫過 FFT 的筆記,所以沒有連結),先看看 FFT 的原理:
- 把 \(a,b\) 變換為 \(A,B\),\(O(n\log n)\);
- 透過 \(C_i=A_iB_i\) 計算,\(O(n)\);
- 把 \(C\) 變換回 \(c\),\(O(n\log n)\)。
綜上,時間複雜度是 \(O(n\log n)\) 的。
在 FFT 中,我們構造了 \(A,B\) 為 \(a,b\) 的點值表示法,這麼做滿足 \(C_i=A_iB_i\) 且容易變換。
其實 FWT 的思想也是一樣的,主要也是需要構造 \(A,B\),使得其滿足 \(C_i=A_iB_i\) 且可以快速變換。下面我們舉 \(\cup\)(按位或)、\(\cap\)(按位與)和 \(\oplus\)(按位異或)為例。
因為數列長度是 \(2\) 的冪會更好處理,所以下文認為數列長度為 \(2^n\)。
按位或
我們可以構造 \(A_i=\sum_{i\cup j=i} a_i\)。看看為什麼需要這麼構造。
首先,它滿足 \(C_i=A_iB_i\):
其次,它可以快速變換。舉順變換的例子。類比 FFT 的步驟,我們採用分治的方法來處理它。假設目前考慮到第 \(i\) 位,其中 \(A_0\) 和 \(A_1\) 是 \(i-1\) 位分治的結果:
其中,\(A_0\) 是數列 \(A\) 的左半部分,\(A_1\) 是 \(A\) 的右半部分。\(\text{merge}\) 函式就是把兩個數列像拼接字串一樣拼接起來。\(+\) 則是將兩個數列對應相加。
這麼做為什麼是正確的呢?容易發現,\(A_0\) 恰好是當前處理到的二進位制位為 \(0\) 的子數列,\(A_1\) 則是當前處理到的二進位制位為 \(1\) 的子數列。若當前位為 \(0\),則只能取二進位制位為 \(0\) 的子數列 \(A_0\) 才能使得 \(i\cup j=i\)。而若當前位為 \(1\),則兩種序列都能取。
考慮逆變換,則是將加上的 \(A_0\) 減回去:
下面我們給出程式碼實現。容易發現順變換和逆變換可以合併為一個函式,順變換時 \(\text{type}=1\),逆變換時 \(\text{type}=-1\)。
void Or(ll *a, ll type) { // 迭代實現,常數更小
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j + k] += a[i + j] * type) %= P;
}
}
}
}
按位與
同理構造 \(A_i=\sum_{i\cap j=i} a_i\)。\(C_i=A_iB_i\) 的正確性不證了。
容易發現,\(A_0\) 恰好是當前處理到的二進位制位為 \(0\) 的子數列,\(A_1\) 則是當前處理到的二進位制位為 \(1\) 的子數列。若當前位為 \(1\),則只能取二進位制位為 \(1\) 的子數列 \(A_0\) 才能使得 \(i\cap j=i\)。而若當前位為 \(0\),則兩種序列都能取。
下面我們給出程式碼實現。順變換時 \(\text{type}=1\),逆變換時 \(\text{type}=-1\)。
void And(ll *a, ll type) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j] += a[i + j + k] * type) %= P;
}
}
}
}
按位異或
發現異或有點難搞,但這怎麼會難倒沃爾什大佬呢?我們引入一個新的運算子 \(\circ\)。定義 \(x\circ y=\text{popcnt}(x\cap y)\bmod 2\),其中 \(\text{popcnt}\) 表示二進位制下 \(1\) 的個數,並重申一下 \(\cap\) 表示按位與。
不用慌,我們也不需要你真正實現一個 \(\text{popcnt}\),它僅僅只是作為一個理解的輔助罷了。
我們發現它滿足 \((x\circ y)\oplus (x\circ z)=x\circ(y\oplus z)\)。(重申一下 \(\oplus\) 表示按位異或)
感性證明:發現這個新的運算子 \(\circ\) 其實就是 \(x\) 與 \(y\) 相同位數的奇偶性。若 \((x\circ y)\oplus (x\circ z)=0\),則 \(x\) 與 \(y\)、\(x\) 與 \(z\) 相同位數個數奇偶性相同,所以 \(y\oplus z\) 和 \(x\) 相同位數個數奇偶性也是相同的 ;若 \((x\circ y)\oplus (x\circ z)=1\),則 \(x\) 與 \(y\)、\(x\) 與 \(z\) 相同位數個數奇偶性不同,所以 \(y\oplus z\) 和 \(x\) 相同位數個數奇偶性也是不同的。
設 \(A_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j\)。我們來證一下 \(C_i=A_iB_i\) 的正確性:
來看看怎麼快速計算 \(A,B\) 的值,依舊是分治:
對於 \(i\) 在當前位為 \(0\) 的子數列 \(A_0\),進行 \(\circ\) 運算時發現它和 \(0\) 計算或和 \(1\) 計算結果都不會變(因為 \(0\cap 0=0,0\cap1=0\)),所以 \(A_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j\) 中的 \(\sum_{i\circ j=1}a_j=0\)。
對於 \(i\) 在當前位為 \(1\) 的子數列 \(A_1\),進行 \(\circ\) 運算時發現它和 \(0\) 計算結果是 \(0\),和 \(1\) 計算結果是 \(1\)(因為 \(1\cap 0=0,1\cap1=1\))。
綜上,有:
也就是:
逆變換易得:
給出程式碼,順變換時 \(\text{type}=1\),逆變換時 \(\text{type}=\frac{1}{2}\)。
void Xor(ll *a, ll type) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j] += a[i + j + k]) %= P;
(a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
(a[i + j] *= type) %= P;
(a[i + j + k] *= type) %= P;
}
}
}
}
現在大家能去切前兩道模板例題,並挑戰一下後面的幾道題目了。
從另一個角度看待 FWT
我們設 \(c(i,j)\) 是 \(a_j\) 對 \(A_i\) 的貢獻係數。我們可以重新描述 FWT 變換的過程:
因為有:
所以我們可以透過簡單的證明得到:\(c(i,j)c(i,k)=c(i,j\odot k)\)。其中 \(\odot\) 是任意一種位運算。
同時,\(c\) 函式還有一個重要的性質,它可以按位處理。
舉個例子,我們變換的時候:
這麼做是比較劣的,我們將其拆分:
考慮前面的式子和後面的式子 \(i,j\) 的區別,發現只有最高位不同。
所以我們將 \(i,j\) 去除最高位的值為 \(i',j'\),並記 \(i_0\) 為 \(i\) 的最高位。有:
如果 \(i_0=0\),則有:
\(i_0=1\) 則有:
也就是說,我們只需要:
四個數就可以完成變換了。我們稱這個矩陣為位矩陣。
如果我們要進行逆變換,則需要上面的位矩陣的逆矩陣。
若逆矩陣為 \(c^{-1}\),可以透過類似操作得到原數:
逆矩陣不一定存在,比如如果有一排 \(0\) 或者一列 \(0\) 那麼這個矩陣就沒有逆,我們在構造時需要格外小心。
按位或
我們可以構造:
這樣滿足 \(c(i,j)c(i,k)=c(i,j\cup k)\)。我們發現,這和我們前面推出的 \(A=\text{merge}(A_0, A_0+A_1)\) 一模一樣!同理,下面也是一個滿足這個條件的矩陣,但我們一般使用上面這個:
雖然下面這個矩陣也滿足 \(c(i,j)c(i,k)=c(i,j\cup k)\),但這個矩陣存在一排 \(0\),不存在逆,所以不合法:
如果我們要進行逆變換,則需要對矩陣求逆,以最上面這個矩陣為例,得:
然後按照順變換的方法,把逆變換矩陣代入即可。
按位與
我們可以構造:
這樣滿足 \(c(i,j)c(i,k)=c(i,j\cap k)\)。
逆矩陣:
按位異或
我們可以構造:
這樣滿足 \(c(i,j)c(i,k)=c(i,j\oplus k)\)。
逆矩陣:
FWT 的性質
FWT 是線性變換。
若 \(FWT(X)\) 是 \(X\) 的 FWT 變換,則有:
以及:
這樣就可以實現快速卷積,參考第四道例題。
K 維 FWT
max 運算
我們重新看看我們的 \(\cup\) 運算,發現他實際上就是二進位制下的取 \(\max\)。我們將其擴充到 \(K\) 進位制,有:
若 \(j=k\),那麼上式又是:
也就是說,每一行的 \(1\) 必定只能在 \(0\) 的前面,如果在後面則不合法了。手玩一下可以發現一組合法構造:
求逆可得:
min 運算
我們重新看看我們的 \(\cap\) 運算,發現他實際上就是二進位制下的取 \(\min\)。我們將其擴充到 \(K\) 進位制,有:
若 \(j=k\),那麼上式又是:
也就是說,每一行的 \(1\) 必定只能在 \(0\) 的後面,如果在後面則不合法了。手玩一下可以發現一組合法構造:
求逆可得:
前兩者用得比較少,用得比較多的是:
不進位加法
我們重新看看我們的 \(\oplus\) 運算,發現他實際上就是二進位制下的不進位加法。我們將其擴充到 \(K\) 進位制,有:
我們構造 \(c(i,j)=\omega_{K}^j\),就可以滿足要求了:
但是每一行都一樣矩陣也沒有逆,所以我們可以構造 \(c(i,j)=\omega_{K}^{(i-1)j}\) 即可。
有下面這個矩陣:
這不就是我們熟悉的範德蒙德矩陣嗎?
現在我們也知道矩陣的逆了:
如果我們題目給出的模數是存在單位根的,我們就可以簡單實現,可以參考第六道例題。
但是單位根在模意義下可能不存在,所以我們考慮擴域,就是人為地定義一個 \(x\),滿足 \(x^K=1\),然後直接把 \(x\) 代入計算,這樣每個數都是一個關於 \(x\) 的 \(k-1\) 次多項式。我們只需要在 \(\pmod {x^K-1}\) 下計算即可。那麼矩陣可以這麼表示:
但是這麼做可能會存在零因子,也就是一個數有多種表示方法,我們無法確定一個數的真實值。
我們考慮不 \(\pmod {x^K-1}\) 了,我們 \(\bmod\) 分圓多項式 \(\Phi_{K}(x)\),他滿足 \(x\) 的階為 \(k\),且在 \(Q\) 上不可約。所以我們定義上面的計算是在 \(\pmod {\Phi_{K}(x)}\) 下進行即可。
另一方面,如何求分圓多項式,這一點可以在因式分解這道題的題解區裡瞭解。這裡給出分圓多項式的表:
還有一個問題是,\(\bmod \Phi_{K}(x)\) 常數大(因為 \(\Phi\) 本身就是一個多項式)。但是因為 \(\Phi_{K}(x)\mid x^k-1\),我們只需要在計算時 \(\bmod x^k -1\),最後再 \(\bmod \Phi_{K}(x)\) 即可。
具體實現參考第七道例題。
例題
「洛谷 P4717」 【模板】快速莫比烏斯/沃爾什變換 (FMT/FWT)
求 \(\cup\)、\(\cap\)、\(\oplus\) 的三種卷積。
\(n\le17\)
這題也就是模板題了,下文直接給出程式碼:
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define P 998244353
const ll N = 1 << 18;
ll n;
ll A[N], B[N];
ll a[N], b[N];
void init() {
for(ll i = 0; i < n; i ++) a[i] = A[i], b[i] = B[i];
}
void Or(ll *a, ll type) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j + k] += a[i + j] * type) %= P;
}
}
}
}
void And(ll *a, ll type) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j] += a[i + j + k] * type) %= P;
}
}
}
}
void Xor(ll *a, ll type) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j] += a[i + j + k]) %= P;
(a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
(a[i + j] *= type) %= P;
(a[i + j + k] *= type) %= P;
}
}
}
}
void calc() {
for(ll i = 0; i < n; i ++) (a[i] *= b[i]) %= P;
}
void print() {
for(ll i = 0; i < n; i ++) printf("%lld ", (a[i] % P + P) % P);
printf("\n");
}
int main() {
scanf("%lld", &n);
n = 1 << n;
for(ll i = 0; i < n; i ++) scanf("%lld", &A[i]);
for(ll i = 0; i < n; i ++) scanf("%lld", &B[i]);
init(); Or(a, 1); Or(b, 1); calc(); Or(a, P - 1); print();
init(); And(a, 1); And(b, 1); calc(); And(a, P - 1); print();
init(); Xor(a, 1); Xor(b, 1); calc(); Xor(a, 499122177); print();
}
「洛谷 P6097」 【模板】子集卷積
求:
\[c_k=\sum_{\substack{{i \cap j=0}\\{i\cup j=k}}} a_i b_j \]\(n\le20\)
首先,下半部分是我們喜聞樂見的 FWT 常見形式,而上半部分我們可以看成是 \(i\) 與 \(j\) 不交。有:
所以我們可以構造:
可以列舉 \(\text{popcnt}\) 的值,分開考慮。
那麼求 \(C\) 的時候有 \(C_{i,j}=\sum_{j=0}^n A_{i,k}B_{i,j-k}\)。
然後就可以做了。
#include <bits/stdc++.h>
using namespace std;
#define popcnt(x) __builtin_popcountll(x)
#define ll long long
const ll M = 20, N = 1 << M, P = 1e9 + 9;
ll n, m;
ll a[M + 1][N], b[M + 1][N], c[M + 1][N];
void Or(ll *a, ll type) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j + k] += a[i + j] * type) %= P;
}
}
}
}
int main() {
scanf("%lld" ,&m);
n = 1 << m;
for(ll i = 0; i < n; i ++) {
scanf("%lld", &a[popcnt(i)][i]);
}
for(ll i = 0; i < n; i ++) {
scanf("%lld", &b[popcnt(i)][i]);
}
for(ll i = 0; i <= m; i ++) {
Or(a[i], 1);
Or(b[i], 1);
}
for(ll i = 0; i <= m; i ++) {
for(ll j = 0; j <= i; j ++) {
for(ll k = 0; k < n; k ++) {
(c[i][k] += a[j][k] * b[i - j][k]) %= P;
}
}
}
for(ll i = 0; i <= m; i ++) {
Or(c[i], -1);
}
for(ll i = 0; i < n; i ++) {
printf("%lld ", (c[popcnt(i)][i] % P + P) % P);
}
}
「牛客 881D」Parity of Tuples
給定 \(n\times m\) 的矩陣 \(a\),定義 \(\text{cnt}(x)\) 為矩陣中有多少行對於 \(x\) 是合法的,合法的定義為這一行中每一個數 \(a_{i,j}\cap x\) 的二進位制值中都有奇數個 \(1\)。
你需要求出對於所有的 \(x\),\(\text{cnt}\) 的取值。
\(n\le10^5,m\le10,x\le2^{20}\)
再次重申,\(\cap\) 是按位與的意思。
首先我們用數學公式定義一下 \(\text{cnt}\)(因為公式複雜,所以加了 \(\tt large\)):
說明一下正確性。如果 \(\text{popcnt}(a_{i,j}\cap x)\) 是奇數的話,那麼 \((-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 的結果就是 \(-1\)。最後 \(1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 就是 \(2\),最後會被 \(\frac{1}{2^m}\) 除去;如果 \(\text{popcnt}(a_{i,j}\cap x)\) 是偶數的話,那麼 \((-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 的結果就是 \(1\)。最後 \(1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 就是 \(0\),那麼整行的結果都是 \(0\)。
然後我們發現它是可以展開的:
然後我們有一個性質:
也就是 \(\sum_{i=1}^n\text{popcnt}(a_i\cap x)\) 的奇偶性與 \(\text{popcnt}((\oplus_{i=1}^na_i)\cap x)\) 的相同。這點在上面的新的運算子 \(\circ\) 的性質中有類似的體現。
容易得到:
我們發現一加一減的可以容斥,我們容斥計算 \(f_i\) 表示 \(n\) 行的所有式子中 \((-1)^i\) 前面的係數和。
// num 處理到第幾列
// x 當前的指數
// mu 當前的係數(+1 or -1)
void dfs(ll *a, ll num, ll x, ll mu) {
if(num > m) {
f[x] += mu;
return;
}
dfs(a, num + 1, x, mu); // 不加入第 num 列,係數不變
dfs(a, num + 1, x ^ a[num], -mu);
}
這樣我們就可以進一步化簡:
我們突然發現後面這個 \((-1)^i\) 取值只有兩種,當 \(x\cap i\) 是奇數時取值為 \(-1\),否則為 \(1\)。
好了,現在我們的問題轉換為了求出:
這不就是 FWT 中的異或變換嗎:
綜上,我們發現這題就是推式子容斥之後得到 FWT 的形式。
原題需要將輸出加密:
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll P = 1e9 + 7;
#define N 100010
#define M 20
#define K 21
ll n, m, k;
ll f[1 << K];
ll a[N][M];
// num 處理到第幾列
// x 當前的指數
// mu 當前的係數(+1 or -1)
void dfs(ll *a, ll num, ll x, ll mu) {
if(num > m) {
f[x] += mu;
return;
}
dfs(a, num + 1, x, mu); // 不加入第 num 列,係數不變
dfs(a, num + 1, x ^ a[num], -mu);
}
void Xor(ll *a, ll type) {
for(ll x = 2; x <= (1 << k); x <<= 1) {
ll z = x >> 1;
for(ll i = 0; i < (1 << k); i += x) {
for(ll j = 0; j < z; j ++) {
(a[i + j] += a[i + j + z]) %= P;
(a[i + j + z] = a[i + j] - 2 * a[i + j + z]) %= P;
(a[i + j] *= type) %= P;
(a[i + j + z] *= type) %= P;
}
}
}
}
ll qpow(ll x, ll y) {
if(y == 0) return 1;
if(y % 2 == 1) return x * qpow(x, y - 1) % P;
ll tmp = qpow(x, y / 2);
return tmp * tmp % P;
}
int main() {
while(scanf("%lld %lld %lld", &n, &m, &k) != EOF) {
for(ll i = 0; i < (1 << k); i ++) f[i] = 0;
for(ll i = 1; i <= n; i ++) {
for(ll j = 1; j <= m; j ++) {
scanf("%lld", &a[i][j]);
}
dfs(a[i], 1, 0, 1);
}
Xor(f, 1);
ll pw = 1, inv = qpow(1 << m, P - 2), ans = 0;
for(ll i = 0; i < (1 << k); i ++) {
ans ^= f[i] * pw % P * inv % P;
(pw *= 3) %= P;
}
printf("%lld\n", ans);
}
}
「AT ABC212H」 Nim Counting
給定兩個數 \(N,K\),以及一個長度為 \(K\) 的整數陣列 \((A_1,A_2,\cdots, A_K)\)。
兩個人玩 Nim 遊戲。
現在透過以下方式生成一個遊戲:
任意選擇一個 \(1\le M\le N\),\(M\) 表示石子堆數。
對於每一堆,其石子數是 \(A\) 中任意一個數。
對於 \(\sum_{i=1}^N K^i\) 種遊戲,求先手獲勝的遊戲數,答案對 \(998244353\) 取模。
\(n\le2\times10^5,K\le2^{16},a_i\le2^{16}\)
根據玩 Nim 遊戲的經驗,可以發現先手獲勝當且僅當 \(\bigoplus_{i=0}^n A_i\neq 0\)。
所以我們定義 dp 式子 \(f_{i,j}\) 表示有 \(i\) 個石堆,且石堆異或和為 \(j\) 的獲勝方案數,有:
答案就是 \(\sum_{i=1}^n\sum_{j\neq0} f_{i,j}\)。
直接轉移是樸素的,發現上面的式子剛好是 FWT 異或操作,也就是:
我們定義 \(a\) 是一個全是 \(1\) 的陣列即可。
同時,我們發現其實不需要真的進行 \(n\) 次卷積,其實只需要將 FWT 變換過之後的結果 \(A\),求出 \(A+A^2+A^3+\cdots+A^n\) 即可。
上面的可以透過等比數列求和公式計算。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define P 998244353
const ll K = 1 << 20;
ll n, k, ans;
ll f[K];
void FWT(ll *a, ll type) {
for(ll x = 2; x <= K; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < K; i += x) {
for(ll j = 0; j < k; j ++) {
(a[i + j] += a[i + j + k]) %= P;
(a[i + j + k] = a[i + j] - 2 * a[i + j + k]) %= P;
(a[i + j] *= type) %= P;
(a[i + j + k] *= type) %= P;
}
}
}
}
ll qpow(ll x, ll y) {
if(y == 0) return 1;
if(y % 2 == 1) return x * qpow(x, y - 1) % P;
ll tmp = qpow(x, y / 2);
return tmp * tmp % P;
}
int main() {
scanf("%lld %lld", &n, &k);
for(ll i = 1; i <= k; i ++) {
ll x;
scanf("%lld", &x);
f[x] ++;
}
FWT(f, 1);
for(ll i = 0; i < K; i ++) {
if(f[i] == 1) f[i] = n;
else {
f[i] = f[i] * (qpow(f[i], n) - 1) % P * qpow(f[i] - 1, P - 2) % P;
}
}
FWT(f, 499122177);
for(ll i = 1; i < K; i ++) {
(ans += f[i]) %= P;
}
printf("%lld", (ans % P + P) % P);
}
「AT ARC100E」 Or Plus Max
給你一個長度為 \(2^n\) 的序列 \(a\),每個\(1\le K\le 2^n-1\),找出最大的 \(a_i+a_j\)(\(i \cup j \le K\),\(0 \le i < j < 2^n\))並輸出。
\(n\le18\)
就是要求 \(\max_{i\cup j=k} a_i+a_j\)。
我們維護 \(f_{i}\) 表示 \(\max_{i\cup j=i} a_i\),\(g_i=\text{max2}_{i\cup j=i} a_i\),\(\text{max2}\) 表示次大值。
然後就像 FWT 的或變換一樣了。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll N = 1 << 18;
ll n;
struct node {
ll mx1, mx2;
node(ll a = 0, ll b = 0):mx1(a), mx2(b) {}
friend node operator +(const node &x, const node &y) {
if(x.mx1 > y.mx1) {
return node(x.mx1, max(x.mx2, y.mx1));
}
return node(y.mx1, max(y.mx2, x.mx1));
}
} a[N];
void FWT(node *a) {
for(ll x = 2; x <= n; x <<= 1) {
ll k = x >> 1;
for(ll i = 0; i < n; i += x) {
for(ll j = 0; j < k; j ++) {
a[i + j + k] = a[i + j] + a[i + j + k];
}
}
}
}
int main() {
scanf("%lld", &n);
n = 1 << n;
for(ll i = 0; i < n; i ++) {
scanf("%lld", &a[i].mx1);
}
FWT(a);
for(ll i = 0; i < n; i ++) {
a[i].mx1 = a[i].mx1 + a[i].mx2;
}
ll ans = 0;
for(ll i = 1; i < n; i ++) {
ans = max(ans, a[i].mx1);
printf("%lld\n", ans);
}
}
「HDU 6618」 Good Numbers
定義一個正整數 \(n\) 是好數當且僅當 \(n\) 在 8 進製表示下所有的數碼出現的次數為 3 的倍數(出現 0 次亦可)。
有多少個 \(k\) 位的 8 進位制數(不含前導 0),滿足這個數是好的,且是 \(p\) 的倍數。對 \(10^9+9\) 取模。
例如:當 \(k=3,p=2\) 時,好數有 \(222,444,666\) 三個。
\(k\le10^{18},p<8\)
考慮狀壓 dp,設 \(f_{i,s,j}\) 表示第 \(i\) 位,\(8\) 種數出現次數對 \(3\) 取模的狀壓情況,以及數對 \(p\) 取模的結果為 \(j\)。
答案就是 \(f_{k,0,0}\)。
直接暴力列舉位數轉移是樸素的,瓶頸在於 \(k\),考慮最佳化掉 \(k\)。
發現我們可以使用像快速冪一樣的方法,也就是倍增 dp。
轉移公式就是:
其中 \(t\) 是轉移的位數,而 \(\oplus\) 在這裡是不進位三進位制加法。
發現這樣多了瓶頸——我們需要列舉 \(s_1\) 和 \(s_2\)。
但是我們發現這不就是 FWT 中異或的形式嗎:\(c_{i\oplus j}\gets a_ib_j\)。考慮三進位制 FWT 加速。下面給出 FWT 的程式碼,w1
是原根的一次方,w2
是原根的二次方:
void FWT(ll *a, ll type) {
for (ll x = 3; x <= N; x *= 3) {
ll k = x / 3;
for (ll i = 0; i < N; i += x) {
for (ll j = 0; j < k; j ++) {
for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k];
if (type == 1) {
tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
} else {
tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P;
}
for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l];
}
}
}
}
因為 \(1e9+9\) 存在原根 \(2\),然後就樸素實現了,注意位矩陣:
程式碼:
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll P = 1e9 + 9;
ll qpow(ll x, ll y) {
if(y == 0) return 1;
if(y % 2 == 1) return x * qpow(x, y - 1) % P;
ll tmp = qpow(x, y / 2);
return tmp * tmp % P;
}
const ll G = 2;
const ll w1 = qpow(G, (P - 1) / 3);
const ll w2 = qpow(G, (P - 1) / 3 * 2);
const ll inv3 = qpow(3, P - 2);
const ll N = 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3;
ll n, p;
ll tmp[8][N], res[8][N], one[8][N];
ll a[8][N], b[8][N];
ll pw3[8];
ll tmp1[3], tmp2[3];
void FWT(ll *a, ll type) {
for (ll x = 3; x <= N; x *= 3) {
ll k = x / 3;
for (ll i = 0; i < N; i += x) {
for (ll j = 0; j < k; j ++) {
for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k];
if (type == 1) {
tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
} else {
tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P;
}
for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l];
}
}
}
}
ll base = 1;
void fun(ll x) {
if(x == 1) {
memset(res, 0, sizeof res);
memset(tmp, 0, sizeof tmp);
memset(one, 0, sizeof one);
for(ll i = 1; i < 8; i ++) res[i % p][pw3[i]] = 1;
for(ll i = 0; i < 8; i ++) tmp[i % p][pw3[i]] = 1;
for(ll i = 0; i < 8; i ++) one[i % p][pw3[i]] = 1;
for (int i = 0; i < p; i ++) {
FWT(tmp[i], 1);
FWT(res[i], 1);
FWT(one[i], 1);
}
base = 8 % p;
return;
}
if(x % 2 == 1) {
fun(x - 1);
memset(a, 0, sizeof a);
memset(b, 0, sizeof b);
for (ll i = 0; i < p; i ++) {
for (ll j = 0; j < p; j ++) {
ll k = (i * 8 + j) % p;
for (ll x = 0; x < N; x ++)
(a[k][x] += tmp[i][x] * one[j][x]) %= P,
(b[k][x] += res[i][x] * one[j][x]) %= P;
}
}
memcpy(tmp, a, sizeof a);
memcpy(res, b, sizeof b);
(base *= 8) %= P;
return;
}
fun(x / 2);
memset(a, 0, sizeof a);
memset(b, 0, sizeof b);
for (ll i = 0; i < p; i ++) {
for (ll j = 0; j < p; j ++) {
ll k = (i * base + j) % p;
for (ll x = 0; x < N; x ++)
(a[k][x] += tmp[i][x] * tmp[j][x]) %= P,
(b[k][x] += res[i][x] * tmp[j][x]) %= P;
}
}
memcpy(tmp, a, sizeof a);
memcpy(res, b, sizeof b);
(base *= base) %= p;
}
int main() {
pw3[0] = 1;
for(ll i = 1; i <= 8; i ++) {
pw3[i] = pw3[i - 1] * 3;
}
while(scanf("%lld %lld", &n, &p) != EOF) {
fun(n);
FWT(res[0], -1);
printf("%lld\n", res[0][0]);
}
}
「CF 1103E」Radix sum
給定一個長度為 \(n\) 的序列 \(a_1,a_2,...,a_n\),對於每一個 \(p \in [0,n-1]\),求滿足下列條件的整數序列 \(i_1,i_2,...,i_n\) 的方案數,對 \(2^{58}\) 取模:
- \(\forall j \in [1,n] , i_j \in [1,n]\);
- \(\sum\limits_{j=1}^n a_{i_j} = p\),這裡的加法定義為十進位制不進位加法。
\(n\le10^5,a_i\le10^5\)
我們可以想到 dp:設計狀態 \(f_{i,s}\) 表示考慮到第 \(i\) 個數,當前加法狀態為 \(s\)。因為 FWT 變換時線性的,可以先變換為 FWT 點值表示法,然後變成自己的 \(n\) 次冪,最後再變換回來。
上面是平凡的,但是題目給出了模數 \(2^{58}\)。發現沒有單位根,所以考慮擴域。
這裡的分圓多項式 \(\Phi_{10}(x)=x^4-x^3+x^2-x+1\)。
然而我們發現 IFWT 時,需要除去進位制 \(10\),然而我們發現 \(10\) 在 \(2^{58}\) 下沒有逆元。實際上我們發現 \(5\) 在 \(2^{58}\) 下是有逆元的:\(57646075230342349\),我們只需要再除去一個 \(2\) 就可以了。設已經除以了 \(5\) 的答案為 \(x\),真正的答案為 \(y\),也就是 \(2^5y\equiv x\pmod{2^{64}}\),顯然,我們有 \(y\equiv \frac{x}{2^5}\pmod{2^{64-5}}\),也就是 \(y\equiv \frac{x}{2^5}\pmod{2^{59}}\),所以直接將最後的答案除以 \(2^5\) 即可。雖然出題人不知道為什麼要模 \(2^{58}\),但再取下模即可。
然後就是平凡實現了:
#include <bits/stdc++.h>
using namespace std;
#define ll unsigned long long
const ll P = 1ull << 58, N = 1e5 + 10;
const ll m = 5, K = 10;
ll inv5;
ll n;
ll pw[m + 1];
ll qpow(ll x, ll y) {
if(y == 0) return 1;
if(y % 2 == 1) return x * qpow(x, y - 1);
ll tmp = qpow(x, y / 2);
return tmp * tmp;
}
struct poly {
ll a[30];
poly() {memset(a, 0, sizeof a);}
ll operator [](ll x) const {return a[x];}
ll& operator [](ll x) {return a[x];}
friend poly operator *(const poly &x, const poly &y) {
poly z;
for(ll i = 0; i < K; i ++) {
for(ll j = 0; j < K; j ++) {
z[(i + j) % K] += x[i] * y[j];
}
}
return z;
}
friend poly operator *(const poly &x, const ll &y) {
poly z;
for(ll i = 0; i < K; i ++) {
z[i] += x[i] * y;
}
return z;
}
friend poly operator +(const poly &x, const poly &y) {
poly z;
for(ll i = 0; i < K; i ++) {
z[i] += x[i] + y[i];
}
return z;
}
poly w(ll x) {
poly res;
for(ll i = 0; i < K; i ++) {
res[(i + x) % K] += a[i];
}
return res;
}
} T, f[N], one;
poly qpow(poly x, ll y) {
if(y == 0) return one;
if(y % 2 == 1) return x * qpow(x, y - 1);
poly tmp = qpow(x, y / 2);
return tmp * tmp;
}
poly tmp1[30], tmp2[30];
void FWT(poly *a, ll type) {
for(ll x = K; x <= pw[m]; x *= K) {
ll k = x / K;
for(ll i = 0; i < pw[m]; i += x) {
for(ll j = 0; j < k; j ++) {
for(ll l = 0; l < K; l ++) tmp1[l] = a[i + j + l * k], tmp2[l] = poly();
if(type == 1) {
for(ll l = 0; l < K; l ++) {
for(ll v = 0; v < K; v ++) {
tmp2[l] = tmp2[l] + tmp1[v].w(l * v % K);
}
}
for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l];
} else {
for(ll l = 0; l < K; l ++) {
for(ll v = 0; v < K; v ++) {
tmp2[l] = tmp2[l] + tmp1[v].w((K - (l * v % K)) % K);
}
}
for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l] * inv5;
}
}
}
}
}
ll mod(poly x){
ll n = 4;
for(ll i = K - 1; i >= n; i --){
ll u = x[i];
for(ll j = 1; j <= n; j ++) x[i - j] -= u * T[n - j];
}
ll u = x[0];
u >>= m;
return u % P;
}
int main() {
pw[0] = 1;
for(ll i = 1; i <= m; i ++) pw[i] = pw[i - 1] * K;
T[0] = 1, T[1] = -1, T[2] = 1, T[3] = -1, T[4] = 1; // 分圓多項式phi10
one[0] = 1;
inv5 = 57646075230342349ull;
scanf("%llu", &n);
for(ll i = 1; i <= n; i ++) {
ll x;
scanf("%llu", &x);
f[x][0] ++;
}
FWT(f, 1);
for(ll i = 0; i < pw[m]; i ++) f[i] = qpow(f[i], n);
FWT(f, -1);
for(ll i = 0; i < n; i ++) cout<<mod(f[i])<<'\n';
}