BZOJ 1004: [HNOI2008]Cards

/*
bzoj 1004: [HNOI2008]Cards http://www.lydsy.com/JudgeOnline/problem.php?id=1004 burside+dp
置换群 k背包dp求解不动点 burside定理 求逆元
*/
#include <cstdio>
#include <algorithm>
using namespace std;
const int Nmax=105;
const int Mmax=65;
int sr,sb,sg,n,m,mod;
int mods[Mmax][Nmax];
int book[Nmax];
int cnt;
int times[Nmax];
int f[Nmax][Nmax][Nmax];

int get_num(int x)
{
cnt=0;
for(int i=1;i<=n;i++)//初始化循环节
{
book[i]=0;
times[i]=0;
}
for(int i=0;i<=sr;i++)//初始化背包
for(int j=0;j<=sg;j++)
for(int k=0;k<=sb;k++)
f[i][j][k]=0;
for(int i=1;i<=n;i++)//求cnt个循环节的长度
{
if(book[i])
continue;
int j=i;
cnt++;
while(!book[j])
{
book[j]=1;
j=mods[x][j];
times[cnt]++;
}
}
f[0][0][0]=1;//背包
for(int i=1;i<=cnt;i++)
for(int r=sr;r>=0;r--)
for(int g=sg;g>=0;g--)
for(int b=sb;b>=0;b--)
{
if(r>=times[i])
f[r][g][b]=(f[r][g][b]+f[r-times[i]][g][b])%mod;
if(g>=times[i])
f[r][g][b]=(f[r][g][b]+f[r][g-times[i]][b])%mod;
if(b>=times[i])
f[r][g][b]=(f[r][g][b]+f[r][g][b-times[i]])%mod;
}
return f[sr][sg][sb];
}

int ex_gcd(int a,int b,int &x,int &y)//solve x,y in a*x+b*y=ex_gcd(a,b,x,y)=gcd(a,b);
{
if(b==0)
{
x=1;
y=0;
return a;
}
int ans=ex_gcd(b,a%b,x,y);
int tmp=x;
x=y;
y=tmp-a/b*y;
return ans;
//x = x0 + (b/gcd)*t
//y = y0 – (a/gcd)*t

}

int get(int a,int m,int c)//get x in a*x=c(mod m)
{
//we can solve x,y in a*x+b*y=c <=> c%gcd(a,b)==0
int x,y;
int gcd=ex_gcd(a,m,x,y);
if(c%gcd!=0)
return -1;//error
x*=c/gcd;
m=abs(m);
int ans=x%m;
while(ans<0)
ans+=m;
return ans;
}

int pow(int base,int n)
{
int ans=1;
while(n>0)
{
if(n&1)
ans=(ans*base)%mod;
base=(base*base)%mod;
n>>=1;
}
return ans;
}

int main()
{
freopen("bzoj1004.in","r",stdin);
scanf("%d%d%d%d%d",&sr,&sb,&sg,&m,&mod);
n=sr+sb+sg;
for(int i=1;i<=m;i++)
for(int j=1;j<=n;j++)
scanf("%d",&mods[i][j]);
m++;
for(int i=1;i<=n;i++)//添加不变置换
mods[m][i]=i;
int num=0;
for(int i=1;i<=m;i++)
num=(num+get_num(i))%mod;
// int x=get(m,mod,1);
int x=pow(m,mod-2);
// printf("x:%d,mod:%d,num:%d\n",x,mod,num);
int ans=(num*x)%mod;
while(ans<0)
ans-=mod;
printf("%d\n",ans);
return 0;
}