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;
}