矩阵乘法总结

这几天做了几道矩阵乘法的题，有一些小技巧，记录一下。

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\n",mid);
62     return 0;
63 }
bzoj2676

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\n",A.a[1][2*n]);
46     return 0;
47 }
bzoj4417

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\n",A.a[1][mm[bit[k]-1]]);
72     return 0;
73 }
bzoj2004

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\n",ans);
80     return 0;
81 }
bzoj3120

 

posted @ 2017-10-30 10:55  Ren_Ivan  阅读(869)  评论(0)  编辑  收藏  举报