题意：有n张牌，有R+G+B=n的3种颜色及其数量，要求用这三种颜色去染n张牌。n张牌有m中洗牌方式，问在不同洗牌方式下本质相同的染色方案数。

解法：这道题非常有意思，题解参考Hzwer学长的。我这里再总结一下：

看到本质相同的染色方案我们很容易会想到Burnside引理和Polya定理，但是这题不能用Polya定理，为什么？因为一般的Ployd染色的颜色个数是没有限制的，于是当循环节为l颜色为c时候，方式数就是c^l（就是因为一个循环方案要相同所以染的颜色也要相同）。但是此题颜色个数有限制，不能直接每个格子有c种选择，所以不能使用Polyd定理。

那现在我们还是得保证一个循环内颜色相同但又不用Ployd呢？我们使用Burnside引理：可以想象成这样，我们必须要有R个红色，G个绿色，B个蓝色，且每一个循环节我们可以选择它染成R/G/B。那这不就是一个01背包模型，每个循环节就是一个物品，RGB就是容量限制，那么我们就可以用01背包计算方案数即可。

细节见代码及其注释。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=25;
int n,r,g,b,m,P,a[65],d[65];

void exgcd(int a,int b,int& d,int& x,int& y) { //ax+by=gcd(a,b) 
    if (!b) { d=a;x=1;y=0; } 
    else { exgcd(b,a%b,d,y,x); y-=x*(a/b); }
}

bool vis[65];
LL dp[65][N][N][N];  //dp[l][i][j][k]代表前l个循环节组成i个Rj个Gk个B的方案数 
LL solve() {  //每次计算置换群a的方案数(相当于做一次01背包) 
    for (int i=1;i<=n;i++) vis[i]=0;
    int num=0,now=1;
    for (int i=1;i<=n;i++) {  //统计循环节 
        if (vis[i]) continue;
        d[++num]=1; now=i;  //循环节数量/大小 
        vis[now]=1;
        while (!vis[a[now]]) {
            d[num]++;
            vis[a[now]]=1;
            now=a[now];
        }
    }
    for (int l=0;l<=num;l++) for (int i=0;i<=r;i++) for (int j=0;j<=g;j++) for (int k=0;k<=b;k++) 
        dp[l][i][j][k]=0;
    dp[0][0][0][0]=1;  //初始化 
    for (int l=1;l<=num;l++)  //循环节个数相当于物品个数 
        for (int i=0;i<=r;i++)
            for (int j=0;j<=g;j++)
                for (int k=0;k<=b;k++) {
                    if (i>=d[l]) dp[l][i][j][k]=(dp[l][i][j][k]+dp[l-1][i-d[l]][j][k])%P;
                    if (j>=d[l]) dp[l][i][j][k]=(dp[l][i][j][k]+dp[l-1][i][j-d[l]][k])%P;
                    if (k>=d[l]) dp[l][i][j][k]=(dp[l][i][j][k]+dp[l-1][i][j][k-d[l]])%P;
                }
    return dp[num][r][g][b];            
}

int main()
{
    scanf("%d%d%d%d%d",&r,&g,&b,&m,&P);
    n=r+g+b;
    LL ans=0;
    for (int i=1;i<=m;i++) {
        for (int j=1;j<=n;j++) scanf("%d",&a[j]);
        ans+=solve();  //累加所有置换方案数 
    }
    for (int i=1;i<=n;i++) a[i]=i;
    ans+=solve();
    int x,y,d; exgcd(m+1,P,d,x,y); 
    x=(x%P+P)%P;  //求出m+1再模P下逆元 
    cout<<ans*x%P<<endl;
    return 0;
}