1 題目相關

1.1傳送門

傳送門

1.2 題目大意

給你一個$m$，現在有一個數組$f$，已知這$2^m$個數在$0$時刻的值$f[0][i](i\in [0,2^m))$，已知另一個數組$b[i](i\in [0,m])$
記$cnt(x)$為$x$在二進制下1的數量，$\oplus$為異或運算
有如下遞推式$$f[i][j]=\sum _{k=0}^{2^m?1}b[cnt[j\oplus k]]?f[i?1][k]\mathrm{m}\mathrm{o}\mathrm{d}p$$
求$f[t][i](i\in [0,2^m))$

1.3 數據范圍

原題數據范圍：
$1\le m\le 20，1\le t\le 10^{18}，2\le p\le 10^9$
$1\le f[0][i],b[i]\le 10^9$
$m,t,p,f,b\in \mathbb{Z}$
時間限制$6s$
為了解釋方便，定義$n=2^m$

2 算法

2.1 暴力

拿到一道題，首先當然先考慮暴力怎么做，當然是很方便的，首先，我們可以直接考慮模擬，總復雜度$\Theta (n^{2}tm)$

2.2 稍有優化的暴力

容易發現，那個$cnt[]$是可以預處理的，所以復雜度變為$\Theta (n^{2}t)$

2.3 換個思路

你發現，$t$實在是太大了，實在為本題的一大瓶頸
考慮矩陣乘法，對于下一個時刻本質上是進行一次矩陣乘法，總復雜度$\Theta (n^{3}logt)$

2.4 你會FWT

前面的那么多算法的復雜度并不優越，一般的出題人對于前面的算法的給分一般都只有一檔，實在不是很可觀
但是，你會FWT，你看過我的博客，(不會FWT的請點擊快速沃爾什變換(FWT))
請仔細觀察如下式子
$$f[i][j]=\sum _{k=0}^{2^m?1}b[cnt[j\oplus k]]?f[i?1][k]\mathrm{m}\mathrm{o}\mathrm{d}p$$
把它做一個轉換
$$f[i][j]=\sum _{a\oplus b=j}b[cnt[a]]?f[i?1][b]\mathrm{m}\mathrm{o}\mathrm{d}p$$
考慮構造一個數組$g[i]=b[cnt[i]]$，下一個時刻就相當于對$g$做一次異或卷積，好吧，目前的復雜度是$\Theta (tnlogn)$

2.5 預處理

由于每次卷積都是卷同一個東西，所以預處理$g^1,g^2,g^3,g^4???$
即快速冪，復雜度$\Theta (nlognlogt)$

2.6 減少FWT、IFWT次數

我們回到2.4，發現復雜度瓶頸在$FWT$和$IFWT$，然而并不需要每次$FWT()$數組對應系數乘后IDFT回來，直接在里面乘就好了，復雜度$\Theta (nlogn+nt)$

2.7 快速冪優化

發現2.6里的乘系數每次都是乘同一個數，快速冪一下，復雜度$\Theta (nlogn+nlogt)$

2.8 Tip

發現這題有一個很坑的地方，就是p可能是2的倍數，而FWT的時候要除以n，然后就會掛
例如：$n=2$,$p=6$,現在有個數當前算出來值為$8$，直接除以$2$結果是$4$，先模$p$結果為$8$，再除以$2$為結果為$1$，就掛了
然后就需要將模數擴大$n$倍
$n?p$是$10^{19}$兩個數直接乘會爆long long，慢速乘又會帶來$\Theta (logn)$的復雜度，不是很優越，所以需要一些奇技淫巧優化這個乘法

inline LL msc(LL a,LL b) {LL v=(a*b-(LL)((long double)a/mod*b+1e-8)*mod);return v<0?v+mod:v; }

代碼

到這里這個算法已經可以通過此題了，貼出AC代碼

#include<cstdio> #include<cctype> #include<algorithm> #define rg register typedef long long LL; template <typename T> inline T max(const T a,const T b){return a>b?a:b;} template <typename T> inline T min(const T a,const T b){return a<b?a:b;} template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;} template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;} template <typename T> inline T abs(const T a){return a>0?a:-a;} template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;} template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);} template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;} template <typename T> inline T square(const T x){return x*x;}; template <typename T> inline void read(T&x) {char cu=getchar();x=0;bool fla=0;while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}while(isdigit(cu))x=x*10+cu-'0',cu=getchar();if(fla)x=-x; } template <typename T> inline void printe(const T x) {if(x>=10)printe(x/10);putchar(x%10+'0'); } template <typename T> inline void print(const T x) {if(x<0)putchar('-'),printe(-x);else printe(x); } LL mod=998244353;const int lenth=1048576; inline LL msc(LL a,LL b) {LL v=(a*b-(LL)((long double)a/mod*b+1e-8)*mod);return v<0?v+mod:v; } inline LL pow(LL a,LL b) {LL res=1;for(;b;a=msc(a,a),b>>=1)if(b&1)res=msc(a,res);return res; } LL n,m,t,a[lenth],b[lenth],c[21];int cnt[lenth]; inline void FWT(LL*A,const int fla) {for(rg int i=1;i<n;i<<=1)for(rg int j=0;j<n;j+=(i<<1))for(rg int k=0;k<i;k++){const LL x=A[j+k],y=A[j+k+i];A[j+k]=(x+y)%mod;A[j+k+i]=(x+mod-y)%mod;}if(fla==-1)for(rg int i=0;i<n;i++)A[i]/=n; } int main() {read(m),read(t),read(mod),n=1<<m,mod*=n;for(rg int i=0;i<n;i++)read(a[i]);for(rg int i=1;i<n;i++)cnt[i]=cnt[i^(i&-i)]+1;for(rg int i=0;i<=m;i++)read(c[i]);for(rg int i=0;i<n;i++)b[i]=c[cnt[i]];FWT(a,1),FWT(b,1);for(rg int i=0;i<n;i++)a[i]=msc(a[i],pow(b[i],t));FWT(a,-1);for(rg int i=0;i<n;i++)print(a[i]),putchar('\n');return 0; }

3 優化

然而，這道題可以進行優化

3.1 加強版數據范圍（毒瘤）

$1\le t\le 2^{10000}$

3.2 思考

對于之前的那個算法，顯然是無法接受這個數據范圍的，所以我們現在需要換個思路

新定義:
定義兩個位置$i$、$j$的距離$dis(i,j)=cnt[i\oplus j]$

對于一個位置$x$，每次能向它轉移的位置$y$有$2^m$個，所以暴力復雜度非常不優。但是你會發現，轉移的位置雖然多，但是系數的種類卻不多，對$x$位置產生貢獻的系數按與$x$位置的距離分只有$m+1$種，我們就可以觀察是否可以在這里進行突破
當然可以，不然我說這個干什么
容易發現，這$m+1$種系數就是輸入的那個$b[]$，而在$t$時刻，仍然是有一個$b_t[]$的(在之后，這個下標代表時刻)，只是數值上發生了變化
那么如果我們能夠快速的求出$t$時刻的b數組，那么我們剩下要做的就是一遍FWT，復雜度只剩$\Theta (nlogn)$了

3.3 分析

要計算$b_t[]$，我們要從$b_{t?1}[]$推過來
轉移顯然是$\Theta (m^2)$的，那么轉移的系數是什么呢？
考慮$b_{t?1}[j]$對$b_t[i]$的貢獻

• step1
  把數值字母化，對于一組$x,y,z$滿足$dis(x,y)=j,dis(x,z)=i$
  $b_{t?1}[j]$相當于$t?1$時刻$x$對$y$的貢獻
  $b_t[i]$相當于$t$時刻$x$對$z$的貢獻
• step2
  已知$dis(x,y)=j$即$z$中對$dis$貢獻為$1$的位有$j$位，為$0$的有$m?j$位
  我們枚舉$y$變到$z$的過程中，貢獻$1$變$0$的位數$s_1$和貢獻$0$變$1$的位數$s_1$
  容易發現$j?s_1+s_0=i$，并且$dis(y,z)=s_0+s_1$，貢獻系數為$$b[s_0+s_1]?C(j,s_1)?C(m?j,s_0)$$
  枚舉$s_0$，那么$s_1=j?i+s_0$，把式子算出來即可

現在知道轉移的系數了，那么就可以用矩陣乘法優化，復雜度是$\Theta (m^{3}logt)$的，算法總復雜度為$\Theta (m^{3}logt+nlogn)$

3.4 代碼

貼出AC代碼

#include<cstdio> #include<cctype> #include<algorithm> #include<cstring> #define rg register typedef long long LL; template <typename T> inline T max(const T a,const T b){return a>b?a:b;} template <typename T> inline T min(const T a,const T b){return a<b?a:b;} template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;} template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;} template <typename T> inline T abs(const T a){return a>0?a:-a;} template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;} template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);} template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;} template <typename T> inline T square(const T x){return x*x;}; template <typename T> inline void read(T&x) {char cu=getchar();x=0;bool fla=0;while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}while(isdigit(cu))x=x*10+cu-'0',cu=getchar();if(fla)x=-x; } template <typename T> inline void printe(const T x) {if(x>=10)printe(x/10);putchar(x%10+'0'); } template <typename T> inline void print(const T x) {if(x<0)putchar('-'),printe(-x);else printe(x); } LL mod=998244353;const int lenth=1048576; inline LL msc(LL a,LL b) {LL v=(a*b-(LL)((long double)a/mod*b+1e-8)*mod);return v<0?v+mod:v; } const int maxn=21; int size=21; struct Martix {LL a[maxn][maxn];Martix(){};Martix(const int x){if(x==0)memset(a,0,sizeof(a));else if(x==1){memset(a,0,sizeof(a));for(rg int i=0;i<size;i++)a[i][i]=1;}}Martix operator *(const Martix&b)const{Martix x=0;for(rg int i=0;i<size;i++)for(rg int j=0;j<size;j++)for(rg int k=0;k<size;k++)x.a[i][j]=(x.a[i][j]+msc(a[i][k],b.a[k][j]))%mod;return x;} }bz,final; template <typename T,typename sum>inline T pow(T x,const sum y) {T res=1;for(rg sum i=1;i<=y;i<<=1,x=x*x)if(i&y)res=res*x;return res; } LL n,m,t,a[lenth],b[lenth],c[21];int cnt[lenth]; inline void FWT(LL*A,const int fla) {for(rg int i=1;i<n;i<<=1)for(rg int j=0;j<n;j+=(i<<1))for(rg int k=0;k<i;k++){const LL x=A[j+k],y=A[j+k+i];A[j+k]=(x+y)%mod;A[j+k+i]=(x+mod-y)%mod;}if(fla==-1)for(rg int i=0;i<n;i++)A[i]/=n; } LL C[21][21]; int main() {read(m),read(t),read(mod),n=1<<m,mod*=n;for(rg int i=0;i<n;i++)read(a[i]);for(rg int i=1;i<n;i++)cnt[i]=cnt[i^(i&-i)]+1;for(rg int i=0;i<=m;i++)read(c[i]);final.a[0][0]=1;for(rg int i=0;i<=m;i++){C[i][0]=1;for(rg int j=1;j<i;j++)C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;C[i][i]=1;}for(rg int j=0;j<=m;j++)for(rg int i=0;i<=m;i++){for(rg int s0=0;s0<=m;s0++){const int s1=j-i+s0;if(s0>m-j)continue;if(s1>j||s1<0)continue;bz.a[i][j]=(bz.a[i][j]+msc(msc(c[s0+s1],C[m-j][s0]),C[j][s1]))%mod;}}final=final*pow(bz,t);for(rg int i=0;i<n;i++)b[i]=final.a[0][cnt[i]];FWT(a,1),FWT(b,1);for(rg int i=0;i<n;i++)a[i]=msc(a[i],b[i]);FWT(a,-1);for(rg int i=0;i<n;i++)print(a[i]),putchar('\n');return 0; }

4 總結

這是一道非常不錯的題，想要AC此題并不是很難，我寫這題的博客是因為它的優化非常巧妙，也沒有用到很難的知識點，值得深思
還有要把慢速乘背下來

總結

以上是生活随笔為你收集整理的Little Pony and Elements of Harmony(CF 453 D)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。