暗了一下polya和burnside定理，感觉还行（就是不会证……）
这题用的是burnside
ans=在每个置换群下不动的方案数之和除以置换数
这题有个难点在取模
关于对p（p为素数）取模（涉及到了除法），我总结了两种方法：
已知x mop p=y，要求x/z mod p=？
大体思路是利用乘法逆，将/z转换成*z的逆元即可
一、利用费马小定理
z^p-1 mod p=1
所以z的逆元=power_mod（z,p-2,p）
二、利用拓展欧几里德算法，即exgcd（z，p，x，y）
while x<0 do inc(x,p)
事实上，这也是乘法逆的思想
exgcd（z，p，x，y） 这实际上是在解同余方程zx mod p=1
x不就是z的逆元吗
但是当p不是素数时，应该怎么办呢？
在做noi2002robot的时候，我运用了一种特殊的方法，最后需要求x>>1 mod p=?
我在过程中始终对p取模，到最后感觉做不下去了
后来我灵机一动，想到了多取一位 在过程中mod 10p
这样 最后输出的时候再mod p 答案就正确了
不知道这样的方法能否应用在更广泛的地方……
ps:刚才试了一下此题
在过程中一直对m*p取模，最后ans div m 再对p取模，竟然AC了
原理在哪里？
代码：
神牛Green Clouds 的代码（费马小定理）：
1 #include <cstdio> 2 #include <algorithm> 3 #include <cstring> 4 5 using namespace std ; 6 7 const int maxn =21; 8 9 int sr , sb , sg , m , mod , f[ maxn ][ maxn ][ maxn ], a[ maxn *3], n ; 10 bool used[ maxn *3]; 11 int v[ maxn *3], cnt , ans =0; 12 13 void update(int&num ,int val ){ 14 ( num += val )%= mod ; 15 } 16 17 int power(int val ,int times ){ 18 int rec =1; 19 for(; times ; times >>=1){ 20 if( times &1)( rec *= val )%= mod ; 21 ( val *= val )%= mod ; 22 } 23 return rec ; 24 } 25 26 int main( ){ 27 scanf("%d%d%d%d%d",&sr ,&sb ,&sg ,&m ,&mod ); 28 n = sr + sb + sg ; 29 for(int w =0; w ++< m +1;){ 30 if( w < m +1)for(int i =0; i ++< n ;)scanf("%d", a + i ); 31 else for(int i =0; i ++< n ;) a[ i ]= i ; 32 memset( used ,false,sizeof( used )); 33 cnt =0; 34 for(int i =0; i ++< n ;)if(! used[ i ]){ 35 v[++ cnt ]=0; 36 for(int j = i ;! used[ j ]; j = a[ j ]){ 37 ++ v[ cnt ], used[ j ]=true; 38 } 39 } 40 memset( f ,0,sizeof( f )); 41 f[0][0][0]=1; 42 for(int h =0; h ++< cnt ;){ 43 for(int i = sr ; i >=0;-- i ){ 44 for(int j = sb ; j >=0;-- j ){ 45 for(int k = sg ; k >=0;-- k ){ 46 if( i >= v[ h ])update( f[ i ][ j ][ k ], f[ i - v[ h ]][ j ][ k ]); 47 if( j >= v[ h ])update( f[ i ][ j ][ k ], f[ i ][ j - v[ h ]][ k ]); 48 if( k >= v[ h ])update( f[ i ][ j ][ k ], f[ i ][ j ][ k - v[ h ]]); 49 } 50 } 51 } 52 } 53 update( ans , f[ sr ][ sb ][ sg ]); 54 } 55 ( ans *=power( m +1, mod -2))%= mod ; 56 printf("%d\n", ans ); 57 return 0; 58 }
exgcd：
1 var s1,s2,s3,n,m,p,ans,i,j,x,y:longint; 2 a:array[0..70,0..70] of longint; 3 f:array[0..25,0..25,0..25] of longint; 4 function mo(x:longint):longint; 5 begin 6 mo:=x mod (p); 7 end; 8 procedure init; 9 begin 10 readln(s1,s2,s3,m,p); 11 n:=s1+s2+s3; 12 for i:=1 to m do 13 for j:=1 to n do 14 read(a[i,j]); 15 inc(m); 16 for i:=1 to n do a[m,i]:=i; 17 end; 18 function dp(x:longint):longint; 19 var i,j,k,tot,y,h,num:longint; 20 v:array[0..70] of boolean; 21 d:array[0..70] of longint; 22 begin 23 tot:=0; 24 fillchar(v,sizeof(v),false); 25 for i:=1 to n do 26 if not(v[i]) then 27 begin 28 num:=1; 29 v[i]:=true; 30 y:=i; 31 while a[x,y]<>i do 32 begin 33 y:=a[x,y]; 34 v[y]:=true; 35 inc(num); 36 end; 37 inc(tot); 38 d[tot]:=num; 39 end; 40 fillchar(f,sizeof(f),0); 41 f[0,0,0]:=1; 42 for h:=1 to tot do 43 for i:=s1 downto 0 do 44 for j:=s2 downto 0 do 45 for k:=s3 downto 0 do 46 begin 47 if i>=d[h] then f[i,j,k]:=mo(f[i,j,k]+f[i-d[h],j,k]); 48 if j>=d[h] then f[i,j,k]:=mo(f[i,j,k]+f[i,j-d[h],k]); 49 if k>=d[h] then f[i,j,k]:=mo(f[i,j,k]+f[i,j,k-d[h]]); 50 end; 51 exit(f[s1,s2,s3]); 52 end; 53 procedure exgcd(a,b:longint;var x,y:longint); 54 var t:longint; 55 begin 56 if a=1 then begin x:=1;y:=0;exit;end; 57 exgcd(b,a mod b,x,y); 58 t:=x;x:=y;y:=t-(a div b)*x; 59 end; 60 procedure main; 61 begin 62 ans:=0; 63 for i:=1 to m do 64 ans:=mo(ans+dp(i)); 65 exgcd(m,p,x,y); 66 while x<0 do inc(x,p); 67 writeln((ans*x) mod p); 68 end; 69 begin 70 init; 71 main; 72 end.
m*p：
1 var s1,s2,s3,n,m,p,ans,i,j,x,y:longint; 2 a:array[0..70,0..70] of longint; 3 f:array[0..25,0..25,0..25] of longint; 4 function mo(x:longint):longint; 5 begin 6 mo:=x mod (m*p); 7 end; 8 procedure init; 9 begin 10 readln(s1,s2,s3,m,p); 11 n:=s1+s2+s3; 12 for i:=1 to m do 13 for j:=1 to n do 14 read(a[i,j]); 15 inc(m); 16 for i:=1 to n do a[m,i]:=i; 17 end; 18 function dp(x:longint):longint; 19 var i,j,k,tot,y,h,num:longint; 20 v:array[0..70] of boolean; 21 d:array[0..70] of longint; 22 begin 23 tot:=0; 24 fillchar(v,sizeof(v),false); 25 for i:=1 to n do 26 if not(v[i]) then 27 begin 28 num:=1; 29 v[i]:=true; 30 y:=i; 31 while a[x,y]<>i do 32 begin 33 y:=a[x,y]; 34 v[y]:=true; 35 inc(num); 36 end; 37 inc(tot); 38 d[tot]:=num; 39 end; 40 fillchar(f,sizeof(f),0); 41 f[0,0,0]:=1; 42 for h:=1 to tot do 43 for i:=s1 downto 0 do 44 for j:=s2 downto 0 do 45 for k:=s3 downto 0 do 46 begin 47 if i>=d[h] then f[i,j,k]:=mo(f[i,j,k]+f[i-d[h],j,k]); 48 if j>=d[h] then f[i,j,k]:=mo(f[i,j,k]+f[i,j-d[h],k]); 49 if k>=d[h] then f[i,j,k]:=mo(f[i,j,k]+f[i,j,k-d[h]]); 50 end; 51 exit(f[s1,s2,s3]); 52 end; 53 procedure main; 54 begin 55 ans:=0; 56 for i:=1 to m do 57 ans:=mo(ans+dp(i)); 58 writeln((ans div m) mod p); 59 end; 60 begin 61 init; 62 main; 63 end. 64
原文：http://www.cnblogs.com/zyfzyf/p/3799170.html
