Polya计数法
关于群、置换、置换群、Burnside定理、Polya定理的详细内容是看得组合数学课件，地址：
http://ishare.iask.sina.com.cn/download/explain.php?fileid=15008967
练习：
poj2409
题意：给一个包含s个珠子的项链，用c种颜色对其染色，问存在多少个不同的等价类。
解：项链可以进行旋转和翻转；
翻转：如果s是奇数，则存在s种置换，每种置换包含s/2+1个循环。
如果s是偶数，存在s/2种以边的中点为中心轴的翻转，每种包含s/2个循环，另外还存在s/2种以点为中心的翻转，每种包含s/2+1个循环；
旋转：旋转i个珠子得到的循环数为gcd(i, s)
证明：设初始值为x，每次旋转i个珠子。有序列：x, x + i, x + i*2, ... , x + i*t
假设x + i*t的时候回到原位置，即 x = (x + i*t)%n -> i*t%n = 0;
也就是 t*i = 0(mod n)
解这个线性同余方程。
d = gcd(n, i);
i'x - n'y = 0/d;
e = x*(0/d)%n = 0;
解的剩余系：e, e + n/d, e + n/d*2 ... e+n/d*(d-1);
发现解的剩余系中第一个有效的解为t = n/d = n/gcd(n, i);
所以旋转i个珠子得到的循环长度为t = n/gcd(i, n),
所以旋转i个珠子得到的循环数为n/t = gcd(i, n);
View Code
//#pragma comment(linker,"/STACK:327680000,327680000") #include <iostream> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #include <algorithm> #include <string> #include <set> #include <functional> #include <numeric> #include <sstream> #include <stack> #include <map> #include <queue> #define CL(arr, val) memset(arr, val, sizeof(arr)) #define REP(i, n) for((i) = 0; (i) < (n); ++(i)) #define FOR(i, l, h) for((i) = (l); (i) <= (h); ++(i)) #define FORD(i, h, l) for((i) = (h); (i) >= (l); --(i)) #define L(x) (x) << 1 #define R(x) (x) << 1 | 1 #define MID(l, r) (l + r) >> 1 #define Min(x, y) (x) < (y) ? (x) : (y) #define Max(x, y) (x) < (y) ? (y) : (x) #define E(x) (1 << (x)) #define iabs(x) (x) < 0 ? -(x) : (x) #define OUT(x) printf("%I64d\n", x) #define Read() freopen("data.in", "r", stdin) #define Write() freopen("data.out", "w", stdout); typedef long long LL; const double eps = 1e-8; const double PI = acos(-1.0); const int inf = ~0u>>2; using namespace std; LL gcd(LL a, LL b) { return b ? gcd(b, a%b) : a; } LL pow(LL a, LL b) { LL res = 1; while(b != 0) { if(b&1) res *= a; a = a*a; b >>= 1; } return res; } int main() { //freopen("data.in", "r", stdin); LL c, s, ans; while(cin >> c >> s) { if(!c && !s) break; ans = 0; for(int i = 1; i <= s; ++i) { ans += pow(c, gcd(s, i)); //printf("%lld %lld %lld %lld\n", c, gcd(s, i), pow(c, gcd(s, i)), ans); } if(s&1) { ans += pow(c, s/2 + 1)*s; } else { ans += pow(c, s/2 + 1)*(s/2); ans += pow(c, s/2)*(s/2); } cout << ans/(s*2) << endl; } return 0; }
poj 2154
这个题题意跟上一题类似，不过没有翻转操作。但是这个题的颜色和珠子数量都在10^9范围内，从1-10^9的范围内枚举会TLE。可以根据一些数学规律对起优化：
上一题的思路是枚举旋转的个数i，得到循环的长度L，现在可以反过来想一下。枚举循环的长度，看存在多少i满足n/gcd(n, i) = L (n能整除L)
设 a = n/L = gcd(n, i), i = a*t.
则当且仅当gcd(L, t) == 1的时候满足:gcd(n, i) = gcd(a*L, a*t) = a
问题转换成求满足gcd(L, t) == 1条件t的个数，前提是 1<= t <= L - 1
也就是求小于L的数中有多少和L互质的数，可以容斥原理，当然也是很明显的eular函数。
可以在[1...sqrt(n)]的范围内枚举L(n%L == 0)，则另外一个被N整除的数就是n/L。这样做比较慢，要跑1300+ ms
代码：
View Code
//#pragma comment(linker,"/STACK:327680000,327680000") #include <iostream> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #include <algorithm> #include <string> #include <set> #include <functional> #include <numeric> #include <sstream> #include <stack> #include <map> #include <queue> #define CL(arr, val) memset(arr, val, sizeof(arr)) #define REP(i, n) for((i) = 0; (i) < (n); ++(i)) #define FOR(i, l, h) for((i) = (l); (i) <= (h); ++(i)) #define FORD(i, h, l) for((i) = (h); (i) >= (l); --(i)) #define L(x) (x) << 1 #define R(x) (x) << 1 | 1 #define MID(l, r) (l + r) >> 1 #define Min(x, y) (x) < (y) ? (x) : (y) #define Max(x, y) (x) < (y) ? (y) : (x) #define E(x) (1 << (x)) #define iabs(x) (x) < 0 ? -(x) : (x) #define OUT(x) printf("%I64d\n", x) #define Read() freopen("data.in", "r", stdin) #define Write() freopen("data.out", "w", stdout); typedef long long LL; const double eps = 1e-8; const double PI = acos(-1.0); const int inf = ~0u>>2; using namespace std; const int N = 32622; int MOD, cnt; int prime[N]; bool vis[N]; void init() { LL i, j; CL(vis, true); for(i = 2; i*i <= 1e9; ++i) { if(vis[i]) { for(j = i + i; j*j <= 1e9; j += i) vis[j] = false; } } cnt = 0; for(i = 2; i*i <= 1e9; ++i) { if(vis[i]) prime[cnt++] = i; } } LL eular(LL n) { int i, res = 1; for(i = 0; i < cnt && prime[i]*prime[i] <= n; ++i) { if(n%prime[i] == 0) { n /= prime[i]; res *= prime[i] - 1; while(n%prime[i] == 0) { n /= prime[i]; res *= prime[i]; } } } if(n > 1) res *= n - 1; return res; } LL exp_mod(LL a, LL b) { LL res = 1; while(b != 0) { if(b&1) {res = (res*a)%MOD;} a = (a*a)%MOD; b >>= 1; } return res; } int main() { //Read(); init(); int T; int l, n, p, ans; scanf("%d", &T); while(T--) { scanf("%d%d", &n, &p); MOD = p; ans = 0; for(l = 1; l*l <= n; ++l) { if(n%l == 0) { ans = (ans + eular(l)*exp_mod(n, n/l - 1))%MOD; //printf("%lld %lld\n", eular(l), exp_mod(n, n/l - 1)); if(l*l == n) continue; ans = (ans + eular(n/l)*exp_mod(n, l - 1))%MOD; //printf("%lld %lld\n", eular(n/l), exp_mod(n, l - 1)); } } printf("%d\n", ans); } return 0; }
可以再进一步优化，上一份代码发现时间主要用在枚举[1...sqrt(n)]范围内的L上面了。
这里可以对n进行质因子分解，然后dfs构造出n的所有因子，时间300+ ms
View Code
//#pragma comment(linker,"/STACK:327680000,327680000") #include <iostream> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #include <algorithm> #include <string> #include <set> #include <functional> #include <numeric> #include <sstream> #include <stack> #include <map> #include <queue> #define CL(arr, val) memset(arr, val, sizeof(arr)) #define REP(i, n) for((i) = 0; (i) < (n); ++(i)) #define FOR(i, l, h) for((i) = (l); (i) <= (h); ++(i)) #define FORD(i, h, l) for((i) = (h); (i) >= (l); --(i)) #define L(x) (x) << 1 #define R(x) (x) << 1 | 1 #define MID(l, r) (l + r) >> 1 #define Min(x, y) (x) < (y) ? (x) : (y) #define Max(x, y) (x) < (y) ? (y) : (x) #define E(x) (1 << (x)) #define iabs(x) (x) < 0 ? -(x) : (x) #define OUT(x) printf("%I64d\n", x) #define Read() freopen("data.in", "r", stdin) #define Write() freopen("data.out", "w", stdout); typedef long long LL; const double eps = 1e-8; const double PI = acos(-1.0); const int inf = ~0u>>2; using namespace std; const int N = 42622; int MOD, cnt, cntf; int prime[N]; int fac[N], num[N]; bool vis[N]; int n, ans; void init() { LL i, j; CL(vis, true); for(i = 2; i*i <= LL(1e9+10); ++i) { if(vis[i]) { for(j = i + i; j*j <= LL(1e9+10); j += i) vis[j] = false; } } cnt = 0; for(i = 2; i*i <= LL(1e9+10); ++i) { if(vis[i]) prime[cnt++] = i; } } void dd(int n) { int i; CL(fac, 0); CL(num, 0); cntf = 0; for(i = 0; i < cnt && prime[i]*prime[i] <= n; ++i) { if(n%prime[i] == 0) { while(n%prime[i] == 0) {n /= prime[i]; num[cntf]++;} fac[cntf++] = prime[i]; } } if(n > 1) {num[cntf]++; fac[cntf++] = n;} } LL eular(LL n) { int i, res = 1; for(i = 0; i < cnt && prime[i]*prime[i] <= n; ++i) { if(n%prime[i] == 0) { n /= prime[i]; res *= prime[i] - 1; while(n%prime[i] == 0) { n /= prime[i]; res *= prime[i]; } } } if(n > 1) res *= n - 1; return res; } LL exp_mod(LL a, LL b) { LL res = 1; while(b != 0) { if(b&1) {res = (res*a)%MOD;} a = (a*a)%MOD; b >>= 1; } return res; } void dfs(int i, int l) { //构造L if(i == cntf) { ans = (ans + eular(l)*exp_mod(n, n/l - 1))%MOD; if(l*l != n) ans = (ans + eular(n/l)*exp_mod(n, l - 1))%MOD; return ; } int tmp = 1, j; for(j = 0; j <= num[i]; ++j) { if(LL(l*tmp)*LL(l*tmp) > n) return ; dfs(i + 1, l*tmp); tmp *= fac[i]; } } int main() { //Read(); init(); int T, p; scanf("%d", &T); while(T--) { scanf("%d%d", &n, &p); MOD = p; ans = 0; dd(n); dfs(0, 1); //cout << endl; printf("%d\n", ans); } return 0; }
POJ 2888
题意：还是给n(n < 10^9)个珠子，m( m <= 10)种颜色，但是现在要求这m种颜色里面有一部分的颜色是不能相邻的。
解：因为m比较小，可以利用颜色关系构图，因为n个珠子构成一个环，想当于从第i种颜色为起点长度为n的回路有多少。《十个利用矩阵乘法解决的经典题目》中的第8题讲的就是这个，如果i，j不能相邻，则矩阵A[i][j] = 0,否则等于1。其他的做法跟上一题一样，枚举L，ans = (ans + eular(L)*A[][]^(n/L))%MOD;
不过ans最后要除掉n，因为共有n种置换，对于除法取模，直接求逆元就可以。
根据费马小定理，当gcd(a, p) = 1时 a^(p -1) ≡ 1 (mod p) -> a*a^(p-2) ≡ 1(mod p)。
结果就是 ans = (ans*exp_mod(n, 9973-2))%9973;
View Code
//#pragma comment(linker,"/STACK:327680000,327680000") #include <iostream> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #include <algorithm> #include <string> #include <set> #include <functional> #include <numeric> #include <sstream> #include <stack> #include <map> #include <queue> #define CL(arr, val) memset(arr, val, sizeof(arr)) #define REP(i, n) for((i) = 0; (i) < (n); ++(i)) #define FOR(i, l, h) for((i) = (l); (i) <= (h); ++(i)) #define FORD(i, h, l) for((i) = (h); (i) >= (l); --(i)) #define L(x) (x) << 1 #define R(x) (x) << 1 | 1 #define MID(l, r) (l + r) >> 1 #define Min(x, y) (x) < (y) ? (x) : (y) #define Max(x, y) (x) < (y) ? (y) : (x) #define E(x) (1 << (x)) #define iabs(x) (x) < 0 ? -(x) : (x) #define OUT(x) printf("%I64d\n", x) #define Read() freopen("data.in", "r", stdin) #define Write() freopen("data.out", "w", stdout); typedef long long LL; const double eps = 1e-8; const double PI = acos(-1.0); const int inf = ~0u>>2; using namespace std; const int M = 11; const int N = 42622; const int MOD = 9973; int cnt, cntf; int prime[N]; int fac[N], num[N]; bool vis[N]; int n, m, ans; struct Mat { int mat[M][M]; }; Mat A, C; Mat operator * (Mat a, Mat b) { Mat c; memset(c.mat, 0, sizeof(c.mat)); int i, j, k; for(k = 0; k < m; ++k) { for(i = 0; i < m; ++i) { if(a.mat[i][k] <= 0) continue; //不要小看这里的剪枝，cpu运算乘法的效率并不是想像的那么理想（加法的运算效率高于乘法，比如Strassen矩阵乘法） for(j = 0; j < m; ++j) { if(b.mat[k][j] <= 0) continue; //剪枝 c.mat[i][j] += a.mat[i][k] * b.mat[k][j]; c.mat[i][j] %= MOD; } } } return c; } Mat operator ^ (Mat a, int k) { Mat c; int i, j; for(i = 0; i < m; ++i) for(j = 0; j < m; ++j) c.mat[i][j] = (i == j); //初始化为单位矩阵 for(; k; k >>= 1) { if(k&1) c = c*a; a = a*a; } return c; } void init() { LL i, j; CL(vis, true); for(i = 2; i*i <= LL(1e9+1000); ++i) { if(vis[i]) { for(j = i + i; j*j <= LL(1e9+1000); j += i) vis[j] = false; } } cnt = 0; for(i = 2; i*i <= LL(1e9+1000); ++i) { if(vis[i]) prime[cnt++] = i; } } void dd(int n) { int i; CL(fac, 0); CL(num, 0); cntf = 0; for(i = 0; i < cnt && prime[i]*prime[i] <= n; ++i) { if(n%prime[i] == 0) { while(n%prime[i] == 0) {n /= prime[i]; num[cntf]++;} fac[cntf++] = prime[i]; } } if(n > 1) {num[cntf]++; fac[cntf++] = n;} } LL eular(LL n) { int i, res = 1; for(i = 0; i < cnt && prime[i]*prime[i] <= n; ++i) { if(n%prime[i] == 0) { n /= prime[i]; res *= prime[i] - 1; while(n%prime[i] == 0) { n /= prime[i]; res *= prime[i]; } } } if(n > 1) res *= n - 1; return res; } LL exp_mod(LL a, LL b) { LL res = 1; while(b != 0) { if(b&1) {res = (res*a)%MOD;} a = (a*a)%MOD; b >>= 1; } return res; } int Count(int k) { int i, j, res = 0; for(i = 0; i < m; ++i) { for(j = 0; j < m; ++j) C.mat[i][j] = A.mat[i][j]; } C = C^k; for(i = 0; i < m; ++i) { res = (res + C.mat[i][i])%MOD; } return res; } void dfs(int i, int l) { if(i == cntf) { ans = (ans + eular(l)*Count(n/l))%MOD; if(l*l != n) ans = (ans + eular(n/l)*Count(l))%MOD; return ; } int tmp = 1, j; for(j = 0; j <= num[i]; ++j) { if(LL(l*tmp)*LL(l*tmp) > n) return ; dfs(i + 1, l*tmp); tmp *= fac[i]; } } int main() { //Read(); init(); int i, j, x, T; scanf("%d", &T); while(T--) { scanf("%d%d%d", &n, &m, &x); for(i = 0; i < m; ++i) { for(j = 0; j < m; ++j) A.mat[i][j] = 1; } while(x--) { scanf("%d%d", &i, &j); i--; j--; A.mat[j][i] = A.mat[i][j] = 0; } ans = 0; dd(n); dfs(0, 1); ans = (ans*exp_mod(n, MOD-2))%MOD; printf("%d\n", ans); } return 0; }
