BZOJ 1004 Cards（polya+dp）

Description
小春现在很清闲,面对书桌上的N张牌,他决定给每张染色,目前小春只有3种颜色:红色,蓝色,绿色.他询问Sun有多少种染色方案,Sun很快就给出了答案.进一步,小春要求染出Sr张红色,Sb张蓝色,Sg张绿色.他又询问有多少种方案,Sun想了一下,又给出了正确答案. 最后小春发明了M种不同的洗牌法,这里他又问Sun有多少种不同的染色方案.两种染色方法相同当且仅当其中一种可以通过任意的洗牌法(即可以使用多种洗牌法,而每种方法可以使用多次)洗成另一种.Sun发现这个问题有点难度,决定交给你,答案可能很大,只要求出答案除以P的余数(P为质数).
Input
第一行输入 5 个整数：Sr,Sb,Sg,m,p(m<=60,m+1 < p < 100)。n=Sr+Sb+Sg。接下来 m 行，每行描述一种洗牌法，每行有 n 个用空格隔开的整数 X1X2…Xn，恰为 1 到 n 的一个排列，表示使用这种洗牌法，第 i位变为原来的 Xi位的牌。输入数据保证任意多次洗牌都可用这 m种洗牌法中的一种代替，且对每种洗牌法，都存在一种洗牌法使得能回到原状态。
Output
不同染法除以P的余数
Sample Input
1 1 1 2 7
2 3 1
3 1 2
Sample Output
2
Solution
m种洗牌就是m种置换，由Burnside引理我们只需要求出这m种置换下的不动点个数，然后加上不洗牌的不动点个数（即n!/(sr!sb!sg!)），相加除以m+1即为答案，对于每种置换求不动点个数，因为一个循环中的卡片颜色需要一样，而每种颜色都有上限，那么相当于求一个三维01背包的方案数
Code

#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
#define INF 0x3f3f3f3f
#define maxn 66
int sr,sb,sg,m,p,a[maxn],b[maxn],res,flag[maxn],dp[maxn][maxn][maxn];
int mod_pow(int a,int b,int p)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
int main()
{
    while(~scanf("%d%d%d%d%d",&sr,&sb,&sg,&m,&p))
    {
        int ans=1,n=sr+sb+sg;
        for(int i=1;i<=n;i++)ans=ans*i%p;
        for(int i=1;i<=sr;i++)ans=ans*mod_pow(i,p-2,p)%p;
        for(int i=1;i<=sb;i++)ans=ans*mod_pow(i,p-2,p)%p;
        for(int i=1;i<=sg;i++)ans=ans*mod_pow(i,p-2,p)%p;
        for(int k=1;k<=m;k++)
        {
            memset(flag,0,sizeof(flag));
            res=0;
            for(int i=1;i<=n;i++)scanf("%d",&a[i]);
            for(int i=1;i<=n;i++)
                if(!flag[i])
                {
                    int cnt=0;
                    while(flag[i]!=i)flag[i]=1,i=flag[i],cnt++;
                    b[res++]=cnt;
                }
            memset(dp,0,sizeof(dp));
            dp[0][0][0]=1;
            for(int t=0;t<res;t++)
                for(int i=1;i<=sr;i++)
                    for(int j=1;j<=sb;j++)
                        for(int k=1;k<=sg;k++)
                        {
                            if(i>=b[t])dp[i][j][k]+=dp[i-b[t]][j][k];
                            if(j>=b[t])dp[i][j][k]+=dp[i][j-b[t]][k];
                            if(k>=b[t])dp[i][j][k]+=dp[i][j][k-b[t]];
                            dp[i][j][k]%=p;
                        }
            ans=(ans+dp[sr][sb][sg])%p;
        }
        ans=ans*mod_pow(m+1,p-2,p)%p;
        printf("%d\n",ans);
    }
    return 0;
}