【BZOJ3120】Line
Description
Wayne喜欢排队……不对，是Wayne所在学校的校长喜欢看大家排队，尤其是在操场上站方阵。
某日课间操时，校长童心大发想了一个极具观赏性的列队方案，如下：
1. 方阵排成N行，每行恰好M个学生。
2. 由于校长喜欢女孩子，所以在一行上不能有连续P个男生。
3. 由于校长喜欢女孩子，所以在校长看来，一列全是男生是不好的，全男生的列数不能超过Q。
Wayne因为感冒了所以不用参加列队，不过他看着大家排队排得不亦乐乎，于是他想知道，在男女生数目无限制的情况下，有多少种列队方案？
两种方案被视作不同，表明存在至少一个二元组(i,j)而两种方案中第i行第j列的同学性别不同。另外，因为答案可能很大，所以请把答案模10^9 + 7。
Input
输入仅一行4个正整数，依次是N,M,P,Q。
Output
输出仅一行，表示答案。
Sample Input
Sample Output
HINT
【数据规模和约定】
对于5%的数据，满足P = 1。
对于另外10%的数据，满足N * M <= 20。
对于另外15%的数据，满足N <= 2，M <= 10^6。
对于另外10%的数据，满足N <= 2。
对于另外20%的数据，满足N <= 4，P <= 2，Q <= 2。
对于100%的数据，满足1 <= N <= 8，1 <= M <= 10^18，1 <= P <= 3，0 <= Q <= 3。
题解：刚看到题，第一反应是状压+矩乘，不过算了算矩阵的大小有点难以接受。不过读了读题，发现行与行之间是互不影响的，所以我们只需要知道每一列男生女生的个数即可。
用f[i][a][b][c]表示第i列有a行有连续0个男生，b行有连续1个男生，n-a-b行有连续2个男生，已经有c行全是男生的方案数，然后把后面那3维压一压变成矩阵就可以用矩乘来优化了。
注意：a+b>n的情况显然不合法，把它扔掉可以压缩矩阵大小；矩乘时特判如果该位置为0，则不进行计算，居然会快很多。
#include <cstdio>
#include <iostream>
#include <cstring>
using namespace std;
typedef long long ll;
const int P=1000000007;
int n,p,q,tot;
ll m,ans;
struct M
{
int v[200][200];
M () {memset(v,0,sizeof(v));}
int * operator [] (int a) {return v[a];}
inline M operator * (const M &b) const
{
M c;
register int i,j,k;
for(i=1;i<=tot;i++) for(k=1;k<=tot;k++) if(v[i][k])
for(j=1;j<=tot;j++) if(b.v[k][j])
c[i][j]=(c[i][j]+(ll)v[i][k]*b.v[k][j])%P;
return c;
}
}S,T;
ll C[10][10];
int _[10][10][5],__[10][5];
//a0 b1 c2
inline void pm(ll y)
{
while(y)
{
if(y&1) S=S*T;
T=T*T,y>>=1;
}
}
int main()
{
scanf("%d%lld%d%d",&n,&m,&p,&q);
int i,j,a,b,a1,b1,x,y;
for(i=0;i<=n;i++)
{
C[i][0]=1;
for(j=1;j<=i;j++) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
}
if(p==1)
{
printf("1");
return 0;
}
if(p==2)
{
for(a=0;a<=n;a++) for(b=0;b<=q;b++) __[a][b]=++tot;
for(a=0;a<=n;a++) for(a1=0;a1<=a;a1++)
{
x=n-a1;
if(a1<n) for(i=0;i<=q;i++) T[__[a][i]][__[x][i]]=C[a][a1];
else for(i=0;i<q;i++) T[__[a][i]][__[x][i+1]]=C[a][a1];
}
S[1][__[n][0]]=1;
pm(m);
for(i=1;i<=tot;i++) ans=(ans+S[1][i])%P;
printf("%lld",ans);
return 0;
}
for(a=0;a<=n;a++) for(b=0;a+b<=n;b++) for(i=0;i<=q;i++) _[a][b][i]=++tot;
for(a=0;a<=n;a++) for(b=0;a+b<=n;b++)
{
for(a1=0;a1<=a;a1++) for(b1=0;b1<=b;b1++)
{
x=n-a1-b1,y=a1;
if(a1+b1<n) for(i=0;i<=q;i++) T[_[a][b][i]][_[x][y][i]]=C[a][a1]*C[b][b1]%P;
else for(i=0;i<q;i++) T[_[a][b][i]][_[x][y][i+1]]=C[a][a1]*C[b][b1]%P;
}
}
S[1][_[n][0][0]]=1;
pm(m);
for(i=1;i<=tot;i++) ans=(ans+S[1][i])%P;
printf("%lld",ans);
return 0;
}//8 1000000000000000000 3 3