我的第一道需要程序建矩阵的矩阵优化DP。
题目可以将不同的p分开处理。
对于p==0 || p==1 直接是0或1
对于p>1，就要DP了。这里以p==3为例：
设dp[i][s1][s2][r]为前i列，结尾为0的有s1行（0表示女生，1表示男生），结尾为01的有s2个，结尾为011的有n-s1-s2个，有r列全是1的方案数。
状态这么复杂，看起来一点也不能用矩阵优化，但我们可以将状态(s1,s2,r)hash成整数，然后建立状态之间的转移。
收获：
这种m超过10^7的一般都要用矩阵优化，如果状态很复杂，可以将复杂的状态（但一般不多）hash成整数，然后用计算机建立状态之间的转移，然后就可以用矩阵优化了。
1 /************************************************************** 2 Problem: 3120 3 User: idy002 4 Language: C++ 5 Result: Accepted 6 Time:24504 ms 7 Memory:2616 kb 8 ****************************************************************/ 9 10 #include <cstdio> 11 #include <cstring> 12 #define M 1000000007 13 #define maxs 200 14 15 typedef long long dint; 16 17 int n, p, q; 18 dint m; 19 int comb[10][10]; 20 21 void init_comb() { 22 for( int i=0; i<=9; i++ ) 23 for( int j=0; j<=i; j++ ) { 24 if( j==0 || i==j ) 25 comb[i][j]=1; 26 else 27 comb[i][j] = (comb[i-1][j]+comb[i-1][j-1]); 28 } 29 } 30 struct Matrix { 31 int n; 32 dint v[maxs][maxs]; 33 void init( int nn ) { 34 n=nn; 35 for( int i=0; i<n; i++ ) 36 for( int j=0; j<n; j++ ) 37 v[i][j]=0; 38 } 39 void make_unit( int nn ) { 40 n=nn; 41 for( int i=0; i<n; i++ ) 42 for( int j=0; j<n; j++ ) 43 v[i][j] = i==j; 44 } 45 Matrix operator*( const Matrix & b ) const { 46 const Matrix & a = *this; 47 Matrix rt; 48 memset( &rt, 0, sizeof(rt) ); 49 rt.n = b.n; 50 for( int k=0; k<n; k++ ) { 51 for( int i=0; i<n; i++ ) { 52 if( a.v[i][k] ) { 53 for( int j=0; j<n; j++ ) { 54 if( b.v[k][j] ) { 55 rt.v[i][j] += (a.v[i][k]*b.v[k][j])%M; 56 if( rt.v[i][j]>=M ) rt.v[i][j]-=M; 57 } 58 } 59 } 60 } 61 } 62 /* 63 for( int i=0; i<n; i++ ) { 64 for( int j=0; j<n; j++ ) { 65 rt.v[i][j] = 0; 66 for( int k=0; k<n; k++ ) { 67 rt.v[i][j] += (a.v[i][k]*b.v[k][j])%M; 68 if( rt.v[i][j]>=M ) 69 rt.v[i][j]-=M; 70 } 71 } 72 } 73 */ 74 return rt; 75 } 76 }; 77 78 Matrix mpow( Matrix a, dint b ) { 79 Matrix rt; 80 for( rt.make_unit(a.n); b; b>>=1,a=a*a ) 81 if( b&1 ) rt=(rt*a); 82 return rt; 83 } 84 namespace Sec1 { 85 void sov(){ 86 if( p==0 ) printf( "0 " ); 87 else if( p==1 ) printf( "1 " ); 88 } 89 } 90 namespace Sec2 { 91 int info[maxs][3], idx[9][4], tot; 92 Matrix trans; 93 void init() { 94 for( int r=0; r<=q; r++ ) 95 for( int a=0; a<=n; a++ ) { 96 info[tot][0] = a; 97 info[tot][1] = n-a; 98 info[tot][2] = r; 99 idx[a][r] = tot; 100 tot++; 101 } 102 trans.init(tot); 103 for( int s=0; s<tot; s++ ) { 104 int a=info[s][0], r=info[s][2]; 105 for( int sa=0; sa<=a; sa++ ) { 106 int na=n-sa, nr=r+(sa==n); 107 if( nr>q ) continue; 108 int ns=idx[na][nr]; 109 trans.v[s][ns] = comb[a][sa]; 110 } 111 } 112 } 113 void sov() { 114 init(); 115 trans = mpow( trans, m ); 116 int row = idx[n][0]; 117 dint ans=0; 118 for( int i=0; i<tot; i++ ) { 119 ans += trans.v[row][i]; 120 if( ans>=M ) ans-=M; 121 } 122 printf( "%lld ", ans ); 123 } 124 } 125 namespace Sec3 { 126 int info[maxs][4], idx[9][9][4], tot; 127 Matrix trans; 128 void init() { 129 for( int r=0; r<=q; r++ ) 130 for( int a=0; a<=n; a++ ) 131 for( int b=0; b<=n-a; b++ ) { 132 info[tot][0] = a; 133 info[tot][1] = b; 134 info[tot][2] = n-a-b; 135 info[tot][3] = r; 136 idx[a][b][r] = tot; 137 tot++; 138 } 139 trans.n = tot; 140 for( int s=0; s<tot; s++ ) { 141 int a=info[s][0], b=info[s][1], c=info[s][2], r=info[s][3]; 142 for( int sa=0; sa<=a; sa++ ) 143 for( int sb=0; sb<=b; sb++ ) { 144 int na=c+a-sa+b-sb, nb=sa, nr=r+(sa+sb==n); 145 if( nr>q ) continue; 146 int ns = idx[na][nb][nr]; 147 trans.v[s][ns] = comb[a][sa]*comb[b][sb]; 148 } 149 } 150 } 151 void sov(){ 152 init(); 153 trans = mpow( trans, m ); 154 int row = idx[n][0][0]; 155 dint ans=0; 156 for( int i=0; i<tot; i++ ) { 157 ans += trans.v[row][i]; 158 if( ans>=M ) ans-=M; 159 } 160 printf( "%lld ", ans ); 161 } 162 } 163 164 int main() { 165 scanf( "%d%lld%d%d", &n, &m, &p, &q ); 166 init_comb(); 167 if( p<=1 ) Sec1::sov(); 168 else if( p==2 ) Sec2::sov(); 169 else if( p==3 ) Sec3::sov(); 170 }