BZOJ1004 [HNOI2008]Cards(Polya计数)

#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<string>
#include<algorithm>
#include<map>
#include<queue>
#include<vector>
#include<cmath>
#include<utility>
using namespace std;
typedef long long LL;
const int N = 68, INF = 0x3F3F3F3F;
#define MS(a, num) memset(a, num, sizeof(a))
#define PB(A) push_back(A)
#define FOR(i, n) for(int i = 0; i < n; i++)
int r, g, b, m, p, n;
int dis[N][N];
bool vis[N];
LL dp[N][N][N];

LL Ext_gcd(LL a,LL b,LL &x,LL &y){//扩展欧几里得
   if(b==0) { x=1, y=0; return a; }
   LL ret= Ext_gcd(b,a%b,y,x);
   y-= a/b*x;
   return ret;
}
LL Inv(LL a,int m){   ///求逆元
   LL d,x,y,t= (LL)m;
   d= Ext_gcd(a,t,x,y);
   if(d==1) return (x%t+t)%t;
   return -1;
}

LL solve(){
    LL ans = 0;
    for(int x = 0; x < m; x++){
        memset(vis, 0, sizeof(vis));
        vector<int> v;
        for(int i = 1; i <= n; i++){
            if(!vis[i]){
                int cnt = 0;
                int tp = i;
                while(!vis[tp]){
                    cnt++;
                    vis[tp] = 1;
                    tp = dis[x][tp];
                }
                v.push_back(cnt);
            }
        }
        memset(dp, 0, sizeof(dp));
        dp[0][0][0] = 1;
        for(int t = 0; t < v.size(); t++){
            for(int i = r; i >= 0; i--){
                for(int j = b; j >= 0; j--){
                    for(int k = g; k >= 0; k--){
                        if(i == 0 && j == 0 && k == 0){
                            continue;
                        }
                        dp[i][j][k] = 0;
                        if(i >= v[t]){
                            dp[i][j][k] = (dp[i][j][k] + dp[i - v[t]][j][k]) % p;
                        }
                        if(j >= v[t]){
                            dp[i][j][k] = (dp[i][j][k] + dp[i][j - v[t]][k]) % p;
                        }
                        if(k >= v[t]){
                            dp[i][j][k] = (dp[i][j][k] + dp[i][j][k - v[t]]) % p;
                        }
                    }
                }
            }
        }
        ans = (ans + dp[r][b][g]) % p;
    }
    ans = ans * Inv(m, p) % p;
    return ans;
}
int main(){
    while(~scanf("%d %d %d %d %d", &r, &b, &g, &m, &p)){
        n = r + b + g;
        bool f = 1;
        for(int i = 0; i < m; i++){
            int cnt = 0;
            for(int j = 1; j <= n; j++){
                scanf("%d", &dis[i][j]);
                if(dis[i][j] == j){
                    cnt++;
                }
            }
            if(cnt == n){
                f = 0;
            }
        }
        if(f){
            for(int i = 1; i <= n; i++){
                dis[m][i] = i;
            }
            m++;
        }
        printf("%lld\n", solve());
    }
    return 0;
}