題意
求長度為 \(n\),元素為 \(1\) 到 \(m\) 之間的整數,且滿足下面兩個條件的序列的數量:
- \(i\) 的出現次數不超過 \(l_i\)
- \(a_{p_i}\neq s_i, \forall i=1,\dots,r\)
\(n\leq 20 ,m\leq 100, r\leq 1000\)
前置知識 1. 子集卷積
首先本題有一個簡單的狀壓 dp 做法:列舉 i 出現的位置集合,然後
這個做法複雜度是 \(O(m3^n)\),顯然是過不了的。
我們令 \(F_i(x) = \sum_{S}[S\cap T_i=\emptyset][|S|\leq l_i]x^S\),則答案為\([x^U]\prod_{i=1}^m F_i(x)\)。\(U\) 表示全集。
注意這裡的卷積是子集卷積而不是子集或卷積。二者的區別
即子集卷積要求沒有重複元素。
子集卷積的一般做法是,將 \(F_i(x)\) 的係數改為 \(n\) 次多項式,\(x^Sy^k\) 的係數是集合為 \(S\),集合元素數量(重複元素算多次)為 \(k\) 的方案數。則 \(F_i(x,y)=\sum_S[S\cap T_i=\emptyset][|S|\leq l_i]x^Sy^{|S|}\)。這樣最後 \(x^Uy^n\) 的係數就是最終答案了。
加一維的好處是,我們已經把無法最佳化的子集卷積轉化為可以用 FWT 最佳化的子集或卷積了。只不過子集或卷積的係數都是多項式。單次子集或卷積的加減次數為 \(O(n2^n)\),而多項式的加減複雜度是 \(O(n)\),所以子集卷積的複雜度為 \(O(n^22^n)\)。
總複雜度從 \(O(m3^n)\) 最佳化到了 \(O(mn^22^n)\),還是無法透過 \(n=20\) 的資料。
PS:子集卷積的另一道板題是 [WC2018]州區劃分
前置知識 2. 集合冪級數的對數和指數運算
係數是整數的集合冪級數是無法定義對數和指數運算的(log 和 exp 在模剩餘系下沒有定義)。但本題集合冪級數的係數都是多項式,因此可藉助多項式的對數和指數運算定義和計算。
定義一個集合冪級數 \(F\) 的對數級數為(要求常數項為 1):
這個級數有意義是因為\(k>n\)時\((F-1)^k\)就是0了。
同理可以定義指數級數為(要求常數項為 0):
同樣的,\(k>n\) 時 \(F^k\) 就是0了,所以級數是有意義的。
集合冪級數的對數和指數怎麼算呢?很簡單,FWT一遍,然後對 \(2^n\) 個多項式分別求對數或指數,再IFWT一遍,就結束了。
現在回到原題。我們需要計算 \(\prod_{i=1}^m F_i(x)\),那麼因為 log 和 exp 可以把連乘和連加互相轉化,所以即計算 \(\exp(\sum_{i=1}^m\log F_i(x))\)。那麼我們只需快速計算所有 \(\log F_i(x)\) 即可。
——我會多項式全家桶!只需\(O(2^nn\log n)\) 就能計算 \(\log F_i\) 了。
——醒醒,\(n\leq 20\)。多項式全家桶跑不過暴力遞推。
關於 log 和 exp 的暴力遞推求法(當然不是硬上定義式求),可以記下這個公式:
那麼直接暴力還是 \(O(mn^22^n)\) 的(要做 \(m2^n\) 次多項式 log)。如何進一步最佳化呢?
本題題解
我們把上面的演算法流程寫下來,看看能發現什麼。
Poly F[M][S],ans[S];
int main(){
// 算 F
for(int i=1;i<=m;++i){
FWT(F[i]);
for(int j=0;j<1<<n;++j)
Log(F[i][j]);
IFWT(F[i]);
for(int j=0;j<1<<n;++j)
ans[j]+=F[i][j];
}
FWT(ans);
for(int j=0;j<1<<n;++j)Exp(ans[j]);
IFWT(ans);
cout<<ans[(1<<n)-1][n]<<endl;
}
看瓶頸部分。我們需要對所有的 \(i,S\) 計算 \(\log([x^S]FWT(F_i))\)。
因此,本質不同的 \([x^S]FWT(F_i)\) 只有 \(O(n^2)\) 種。我們可以預處理出這些多項式的 log。這樣需要用的時候直接取出來即可。
另外,因為 FWT 是線性的,所以 log 的 IFWT 和 exp 的 FWT 可以抵消掉。這樣我們就最佳化掉了對 \(F_i\) 的所有變換。
總複雜度 \(O(n^4+mn2^n+n^22^n)\)。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=21,M=105,S=1<<20,mod=998244353;
int sub(int x,int y){return x^(x&y);}
ll inv[N],c[N][N];
valarray<ll> Log(const valarray<ll>& f){
int n=f.size()-1;
valarray<ll> g(n+1);
for(int i=1;i<=n;++i){
ll tmp=0;
for(int j=1;j<i;++j)
tmp=(tmp+j*g[j]%mod*f[i-j])%mod;
g[i]=(f[i]-inv[i]*tmp%mod+mod)%mod;
}
return g;
}
valarray<ll> Exp(const valarray<ll>& f){
int n=f.size()-1;
valarray<ll> g(n+1);
g[0]=1;
for(int i=1;i<=n;++i){
for(int j=1;j<=i;++j)
g[i]=(g[i]+j*f[j]%mod*g[i-j])%mod;
g[i]=g[i]*inv[i]%mod;
}
return g;
}
void FWT_OR(valarray<ll>* a,int n){
for(int k=1;k<n;k<<=1)
for(int i=0;i<n;i+=k<<1)
for(int j=0;j<k;++j)
(a[i+j+k]+=a[i+j])%=mod;
}
void IFWT_OR(valarray<ll>* a,int n){
for(int k=1;k<n;k<<=1)
for(int i=0;i<n;i+=k<<1)
for(int j=0;j<k;++j)
(a[i+j+k]+=mod-a[i+j])%=mod;
}
valarray<ll> g[N][N];
void init(int n){
inv[1]=1;
for(int i=2;i<=n;++i)inv[i]=inv[mod%i]*(mod-mod/i)%mod;
for(int i=0;i<=n;++i){
c[i][0]=1;
for(int j=1;j<=i;++j)
c[i][j]=c[i-1][j]+c[i-1][j-1];
}
for(int i=0;i<=n;++i)
for(int j=0;j<=n;++j){
g[i][j].resize(n+1);
for(int t=0;t<=j;++t)g[i][j][t]=c[i][t];
g[i][j]=Log(g[i][j]);
}
}
int n,m,lim[M],q,t[M],x,y;
valarray<ll> f[S];
int main(){
ios::sync_with_stdio(0);cin.tie(0);
cin>>n>>m;
init(n);
for(int i=1;i<=m;++i)cin>>lim[i];
cin>>q;
while(q--)cin>>x>>y,--x,t[y]|=1<<x;
for(int i=0;i<1<<n;++i){
f[i].resize(n+1);
for(int j=1;j<=m;++j)
f[i]=(f[i]+g[__builtin_popcount(sub(i,t[j]))][lim[j]])%mod;
f[i]=Exp(f[i]);
}
IFWT_OR(f,1<<n);
cout<<f[(1<<n)-1][n]<<'\n';
}
valarray是個好東西。用起來有點py那味了(