POJ_2888
自己一开始学polya的时候还是太浮躁了，只是相当于背过了个公式而已，其实根本没建立在理解的基础之上。今天踏下心来又看了一遍黑书，终于能够自己根据burnside引理来推导出polya公式了。
解决这个题目首先要了解burnside引理的内容，接着就是用dp的方法去具体计算黑书上所谓的C(f)（在置换f下保持不变的着色方案数）了。最后的结果是C(f)的平均值，也就是C(f)/N，由于9973和N互质，所以可以先求N的逆元，进一步就可以得到(C(f)/N)%9973的值了。
推荐一个此题讲解得非常详细的博客：。
#include#include #include #include #define MAXD 40010#define MAXM 15#define MOD 9973using namespace std;int N, M, K, prime[MAXD], isprime[MAXD], P, g[MAXM][MAXM], d[MAXD], D;struct Matrix{ int a[MAXM][MAXM]; void init() { memset(a, 0, sizeof(a)); }}unit, mat;void prepare(){ int i, j, k = 40000; P = 0; memset(isprime, -1, sizeof(isprime)); for(i = 2; i <= k; i ++) if(isprime[i]) { prime[P ++] = i; for(j = i * i; j <= k; j += i) isprime[j] = 0; }}int euler(int n){ int i, ans, x; ans = x = n; for(i = 0; i < P && prime[i] * prime[i] <= n; i ++) if(x % prime[i] == 0) { ans = ans / prime[i] * (prime[i] - 1); while(x % prime[i] == 0) x /= prime[i]; } if(x > 1) ans = ans / x * (x - 1); return ans;}void exgcd(int a, int b, int &x, int &y){ if(b == 0) x = 1, y = 0; else exgcd(b, a % b, y, x), y -= x * (a / b);}Matrix multiply(Matrix &x, Matrix &y){ int i, j, k; Matrix z; z.init(); for(i = 0; i < M; i ++) for(k = 0; k < M; k ++) if(x.a[i][k]) { for(j = 0; j < M; j ++) if(y.a[k][j]) z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j]) % MOD; } return z;}void init(){ int i, x, y; scanf("%d%d%d", &N, &M, &K); memset(g, -1, sizeof(g)); for(i = 0; i < K; i ++) { scanf("%d%d", &x, &y); -- x, -- y; g[x][y] = g[y][x] = 0; }}void divide(int n){ int i, j; D = 0; for(i = 1; i * i <= n; i ++) if(n % i == 0) { d[D ++] = i; if(n / i != i) d[D ++] = n / i; } sort(d, d + D);}void powmod(Matrix &unit, Matrix &mat, int n){ while(n) { if(n & 1) unit = multiply(mat, unit); n >>= 1; mat = multiply(mat, mat); }}void solve(){ int i, j, x, y, n, ans = 0; divide(N); for(i = 0; i < D; i ++) { n = euler(N / d[i]) % MOD; for(x = 0; x < M; x ++) for(y = 0; y < M; y ++) mat.a[x][y] = -g[x][y]; unit = mat; powmod(unit, mat, d[i] - 1); for(j = 0; j < M; j ++) ans = (ans + n * unit.a[j][j]) % MOD; } exgcd(N, MOD, x, y); x = (x % MOD + MOD) % MOD; ans = (ans * x) % MOD; printf("%d\n", ans);}int main(){ int t; prepare(); scanf("%d", &t); while(t --) { init(); solve(); } return 0;}