这几天做了几道矩阵乘法的题，有一些小技巧，记录一下。
1.bzoj2676随时记录答案
这是需要随时统计答案，首先，肯定是二分答案，如何判定呢
显然节点数应该是Q*R的，记录当前的命和连胜次数
但是期望得分还是没有办法转移，考虑新建一个节点用来记录答案，
在每个节点，我们统计的是当前这步出现这个状态的概率，
原有节点直接按照概率转移，构造矩阵，
那么同时在每个节点，他也有x的几率对答案做出贡献，直接修改矩阵就好了，最后a[1][tot]就是最终的答案
1 #include<cstdio> 2 #include<cstring> 3 #include<iostream> 4 #include<algorithm> 5 #include<cmath> 6 #define eps 1e-8 7 #define min(x,y) ((x>=y)?(y):(x)) 8 #define id(x,y) (min(x-1,n1-1)*n2+min(y,n2)) 9 using namespace std; 10 int n,n1,n2,m,s; 11 struct mart{ 12 double a[102][102]; 13 mart(){memset(a,0,sizeof a);} 14 inline mart operator * (mart b){ 15 mart c; 16 for(int j=1;j<=m;j++) 17 for(int i=1;i<=m;i++){ 18 if(a[i][j]<1e-11)continue; 19 for(int k=1;k<=m;k++){ 20 if(b.a[j][k]<1e-11)continue; 21 c.a[i][k]+=a[i][j]*b.a[j][k]; 22 } 23 } 24 return c; 25 } 26 }; 27 inline mart qp(mart a,int b){ 28 mart c; 29 for(int i=1;i<=m;i++)c.a[i][i]=1; 30 while(b){ 31 if(b&1)c=c*a; 32 a=a*a;b>>=1; 33 } 34 return c; 35 } 36 inline bool check(double x){ 37 mart A,B; 38 B.a[m][m]=1; 39 for(int i=1;i<=n1;i++) 40 for(int j=1;j<=n2;j++) 41 { 42 B.a[id(i,j)][id(i+1,j+1)]=x; 43 B.a[id(i,j)][m]=x*i; 44 if(j!=1)B.a[id(i,j)][id(1,j-1)]=1-x; 45 } 46 B=qp(B,n); 47 A.a[1][n2]=1; 48 A=A*B; 49 return s<A.a[1][m]; 50 } 51 int main(){ 52 scanf("%d%d%d%d",&n,&n1,&n2,&s); 53 m=n1*n2+1; 54 double l=0,r=1,mid; 55 if(!check(r)){puts("Impossible.");return 0;} 56 while(r-l>eps){ 57 mid=(l+r)/2.0; 58 if(check(mid))r=mid; 59 else l=mid; 60 } 61 printf("%0.6lf ",mid); 62 return 0; 63 }
2.bzoj4417前缀和
可以发现，这可转移需要用到和他奇偶性不同的答案的前缀和，
考虑怎么在矩阵中记录前缀和，
那么，我们将矩阵边长乘二，我们在前n个记录与他奇偶性不同的前缀和，后n个记录相同的前缀和
转移时首先要交换这两部分，还要将不同的前缀和更新出的答案加到当前这项的前缀和上，
要注意最后一步转移时只要统计答案而不需转移前缀和就好了
1 #include<cstdio> 2 #include<cstring> 3 #include<iostream> 4 #include<algorithm> 5 #include<cmath> 6 #define mod 30011 7 using namespace std; 8 int n,m; 9 struct mart{ 10 int a[105][105]; 11 mart(){memset(a,0,sizeof a);} 12 mart operator * (mart b){ 13 mart c; 14 for(int j=1;j<=2*n;j++) 15 for(int i=1;i<=2*n;i++){ 16 c.a[i][j]=0; 17 for(int k=1;k<=2*n;k++) 18 c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod; 19 } 20 return c; 21 } 22 }A,B,C; 23 mart qp(mart a,int b){ 24 mart c; 25 for(int i=1;i<=2*n;i++)c.a[i][i]=1; 26 while(b){ 27 if(b&1)c=c*a; 28 a=a*a; b>>=1; 29 } 30 return c; 31 } 32 int main(){ 33 scanf("%d%d",&n,&m); 34 for(int i=1;i<=n;i++){ 35 B.a[i][n+i]=B.a[n+i][i]=1; 36 for(int j=-1;j<=1;j++) 37 if(i+j>=1&&i+j<=n)B.a[n+i][n+i+j]=1; 38 } 39 A.a[1][1]=1; 40 C=qp(B,m-1); 41 A=A*C; 42 for(int i=1;i<=n;i++) 43 B.a[i][n+i]=0; 44 A=A*B; 45 printf("%d ",A.a[1][2*n]); 46 return 0; 47 }
3.bzoj2004合理状态便于转移
这题是状态压缩+矩阵乘，关键就是想怎么定义状态，怎么方便转移
我先想的是停这个车站的车在前面哪个车站停了，发现无法转移。又想了个十进制状压，233，
后来想到了接近正解的定义方法，f[i][S]表示第一辆车走到i，他后面p个车站是否有车经过，我尝试先让第一辆走，但是依旧很乱，无法转移
那么不如修改一下，f[i][S]表示最后一辆车走到i，前面p个车站是否有车停留。这样我们强行让最后一辆车走就可以了，
注意，因为我们是按顺序递推，不能跳过任何一步，如果他前面那个车站没车，他一定只能走一步，否则枚举任意一个空位转移就好了
1 #include<cstdio> 2 #include<cstring> 3 #include<iostream> 4 #include<algorithm> 5 #include<cmath> 6 #include<map> 7 #define mod 30031 8 using namespace std; 9 map<int,int> mm; 10 int n,m,k,p; 11 int bit[13]; 12 struct mart{ 13 int a[130][130]; 14 mart(){memset(a,0,sizeof a);} 15 mart operator * (mart b){ 16 mart c; 17 for(int i=1;i<=m;i++) 18 for(int j=1;j<=m;j++){ 19 if(a[i][j]) 20 for(int k=1;k<=m;k++) 21 c.a[i][k]=(c.a[i][k]+a[i][j]*b.a[j][k])%mod; 22 } 23 return c; 24 } 25 }A,B; 26 int getnum(int x){ 27 int cnt=0; 28 while(x){cnt+=x&1;x>>=1;} 29 return cnt; 30 } 31 int tmpa[130]; 32 void multyA(){ 33 memset(tmpa,0,sizeof tmpa); 34 for(int i=1;i<=m;i++) 35 for(int j=1;j<=m;j++) 36 tmpa[i]=(tmpa[i]+A.a[1][j]*B.a[j][i])%mod; 37 for(int i=1;i<=m;i++) 38 A.a[1][i]=tmpa[i]; 39 } 40 mart qp(int x){ 41 while(x){ 42 if(x&1)multyA(); 43 B=B*B;x>>=1; 44 } 45 } 46 bool vis[2050]; 47 int main(){ 48 bit[0]=1; 49 for(int i=1;i<=11;i++)bit[i]=bit[i-1]<<1; 50 scanf("%d%d%d",&n,&k,&p); 51 for(int i=0,j;i<bit[p];i++)if((getnum(i)==k)&&(i&1)){ 52 int v=(i>>1),vv; 53 if((v&1)==0){ 54 v|=1; 55 if(!vis[i]){mm[i]=++m;vis[i]=1;} 56 if(!mm[v]){mm[v]=++m;vis[v]=1;} 57 B.a[mm[i]][mm[v]]++; 58 } 59 else{ 60 for(j=1;j<=p;j++) 61 if((v&bit[j-1])==0){ 62 vv=v|bit[j-1]; 63 if(!vis[i]){mm[i]=++m;vis[i]=1;} 64 if(!mm[vv]){mm[vv]=++m;vis[vv]=1;} 65 B.a[mm[i]][mm[vv]]++; 66 } 67 } 68 } 69 A.a[1][mm[bit[k]-1]]=1; 70 qp(n-k); 71 printf("%d ",A.a[1][mm[bit[k]-1]]); 72 return 0; 73 }
4.bzoj3120去除冗余状态，合并相同状态
这题乍一看和上一题一样，但是要记录连续的三列，$O(2^{8^{3^3}})$发现不可做
会发现这样有好多状态其实是本质相同的，所以我们只需要记录分别有多少个位置是xx0，x01，011，而不需要具体知道那几位是他们，
转移时枚举这列出现多少个1，又有多少个1放在了x01的后面，组合数转移即可
1 #include<cstdio> 2 #include<cstring> 3 #include<iostream> 4 #include<algorithm> 5 #include<cmath> 6 #include<map> 7 #define LL long long 8 #define mod 1000000007 9 #define id(_,__,___) ((233*(_))+(223*(__))+(211*(___))) 10 using namespace std; 11 struct data{int x,y,z;}d[2005];//x:出现过几列都是1 y:xx0 z:x01 n-y-z:011 12 map<int,int> mm; 13 int n,p,q,tot,ans; 14 LL m; 15 int C[15][15]; 16 void init(){ 17 for(int i=0;i<=10;i++){ 18 C[i][0]=1; 19 for(int j=1;j<=i;j++) 20 C[i][j]=C[i-1][j-1]+C[i-1][j]; 21 } 22 } 23 struct mart{ 24 int a[180][180]; 25 mart(){memset(a,0,sizeof a);} 26 mart operator * (mart b){ 27 mart c; 28 for(int i=1;i<=tot;i++) 29 for(int j=1;j<=tot;j++) if(a[i][j]) 30 for(int k=1;k<=tot;k++) 31 c.a[i][k]=(c.a[i][k]+1ll*a[i][j]*b.a[j][k]%mod)%mod; 32 return c; 33 } 34 }A,B; 35 int tmpa[180]; 36 void multya(){ 37 memset(tmpa,0,sizeof tmpa); 38 for(int i=1;i<=tot;i++) 39 for(int j=1;j<=tot;j++) 40 tmpa[i]=(tmpa[i]+1ll*A.a[1][j]*B.a[j][i]%mod)%mod; 41 for(int i=1;i<=tot;i++) 42 A.a[1][i]=tmpa[i]; 43 } 44 void qp(LL b){ 45 while(b){ 46 if(b&1)multya(); 47 B=B*B;b>>=1; 48 } 49 } 50 int main(){ 51 init(); 52 scanf("%d%lld%d%d",&n,&m,&p,&q); 53 for(int t=0;t<=q;t++) 54 for(int i=0;i<=n;i++) 55 for(int j=0,k;j<=n;j++)if(i+j<=n){ 56 k=n-i-j; 57 if((i==0)&&(t==0))continue; 58 if((i==0&&j==0)&&(t<2))continue; 59 if((k)&&(p==2))continue; 60 if((k||j)&&(p==1))continue; 61 mm[id(t,i,j)]=++tot; 62 d[tot]=(data){t,i,j}; 63 } 64 for(int i=1,x1,y1,z1,w1;i<=tot;i++){ 65 x1=d[i].x,y1=d[i].y,z1=d[i].z,w1=n-y1-z1; 66 for(int j=0;j<=y1+z1;j++){ 67 if(j==n&&x1==q)continue; 68 for(int k=0;k<=j&&k<=z1;k++){//共添加j个1，其中k个添加到x01后 69 if(j-k>y1)continue; 70 if(j==n)B.a[mm[id(x1,y1,z1)]][mm[id(x1+1,n-j,j-k)]]+=C[z1][k]*C[y1][j-k]; 71 else B.a[mm[id(x1,y1,z1)]][mm[id(x1,n-j,j-k)]]+=C[z1][k]*C[y1][j-k]; 72 } 73 } 74 } 75 A.a[1][mm[id(0,n,0)]]=1; 76 qp(m); 77 for(int i=1;i<=tot;i++) 78 ans=(ans+A.a[1][i])%mod; 79 printf("%d ",ans); 80 return 0; 81 }