burnside & polya 第二弹 poj四题

 

  一个白天就耗在这么几道题中，表示效率太低了=_=。

   今天天气闷死了， 热还比较可以忍受，可是闷就实在无法忍受，下午的时候感觉骨头都要散掉了=_=! 被最后一道题虐了一下午=_=！

   poj 1286

  给一个长度为n的环染色， 只有三种颜色，考虑旋转同构和翻转同构， n<= 24

  其实置换方式只有2*n种，且明显构成置换群， 原本想直接枚举置换，用polya定理，裸找出循环节数即可，后来经cwx的提醒，其实可以手工算出置换的循环节数。

  旋转x个的置换，其循环节数位gcd（n，x）， n为奇数时，沿对称轴翻转的置换，循环节数为n/2+1， n为偶数是，n/2的翻转置换，其循环节数位n/2， 另外n/2的循环节数位n/2+1.

  那么代码会变得相当简单！ 看来polya的题“人类智慧”是很重要的，即使可以用计算机裸做，手算一下也会有不一样的惊喜。

  code：

# include <cstdlib># include <cstdio># include <cmath>using namespace std;typedef long long int64;int n, i;int64 ans;int64 mul(int x, int y){int64 ask = 1;while (y--) ask = ask * x;return ask;}int gcd(int x, int y) {return !x? y: gcd(y%x, x);}int main(){freopen("1286.in", "r", stdin);freopen("1286.out", "w", stdout);for (;scanf("%d", &n), n!= -1;){if (n == 0) printf("0\n");else if (n & 1){ans = mul(3, n/2+1) * n;for (i = 1; i <= n; i++)  ans += mul(3, gcd(i, n));ans /= 2*n;printf("%I64d\n", ans);}else {ans = mul(3, n/2) * n/2;ans += mul(3, n/2+1) * n/2;for (i = 1; i <= n; i++)  ans += mul(3, gcd(i, n));ans /= 2*n;printf("%I64d\n", ans);}}return 0;}

  poj2409 

  和上题如出一辙，只是颜色数是改变的，方法和上面是一模一样的。

   code：

  

# include <cstdlib># include <cstdio># include <cmath>using namespace std;typedef long long int64;int c, s, i;int64 ans;int64 mul(int x, int y){int64 ask = 1;while (y--) ask = ask * x;return ask;}int gcd(int x, int y) {return !x? y: gcd(y%x, x);}int main(){freopen("2409.in", "r", stdin);freopen("2409.out", "w", stdout);for (;scanf("%d%d", &c, &s), s!= 0;){if (s & 1){ans = mul(c, s/2+1) * s;for (i = 1; i <= s; i++)  ans += mul(c, gcd(i, s));ans /= 2*s;printf("%I64d\n", ans);}else {ans = mul(c, s/2) * s/2;ans += mul(c, s/2+1) * s/2;for (i = 1; i <= s; i++)  ans += mul(c, gcd(i, s));ans /= 2*s;printf("%I64d\n", ans);}}return 0;}

poj2154：

  长度为n的环，染n种颜色。不过这次只考虑旋转（不考虑沿对称轴翻转）引起的同构！ 

  n = 10^9 而且好心的给了3500组测试数据。

   

  和上面的式子是一样的，我们需要求sum（i^gcd（n，i）） for i in [1， n]

  即使使用快速幂都避免不了tle的命运。

  其实所有的gcd（n，i）的种类是不多的，只有n的约数级别， 一个数的约数是什么级别的？ 我直接的平均是O（logn），应该不会超过O（sqrt（n））

   而如何求出gcd（n，i） = k 有多少个i满足条件？ 很容易想到是phi（n/i）。

   可以通过弄出n的所有质因数，用一遍dfs得到所有n的约数的phi。

   这样求循环节数的复杂度就降低到了O（sqrt（n））， 然后乘上O（logn）的快速幂，勉勉强强可以过。

 // polya中存在除法操作，但是由于mod的数不是质数，不过可以通过分子的式子中的指数全-1即可。

# include <cstdlib># include <cstdio># include <cmath>using namespace std;typedef long long int64;const int sqN = 100000+10;int n, mo;int tot, phi[sqN], pri[sqN];bool step[sqN];int Div[sqN], dphi[sqN], num[sqN], have[sqN];int64 ans;void prepare_prime(){int i, j;phi[1] = 1; for (i = 2; i <= 100000; i++){if (!step[i])  step[pri[++tot] = i] = 1, phi[i] = i-1;for (j = 1; j <= tot; j++){  if (i*pri[j] > 100000) break;  step[i*pri[j]] = 1;  if (i % pri[j] == 0)  {  phi[i*pri[j]] = phi[i] * pri[j];  break;  }  else phi[i*pri[j]] = phi[i] * (pri[j]-1);}}}void get_prime(){int i, it;for (i = 1; i <= have[0]; i++) have[i] = num[i] = 0;have[0] = 0; it = n;for (i = 1; i <= tot && pri[i] <= n; i++)if (!(n % pri[i])){  have[++have[0]] = pri[i];  while (!(it % pri[i])) it/= pri[i], num[have[0]]++;}if (it != 1) have[++have[0]] = it, num[have[0]]++;}void dfs(int it, int now, int fai){int i, fn, fp;if (it == have[0]+1){Div[++Div[0]] = now; dphi[Div[0]] = fai;return;}dfs(it+1, now, fai);fn = now*have[it]; fp = fai*(have[it]-1);dfs(it+1, fn, fp);for (i = 2; i <= num[it]; i++){  fn*= have[it]; fp*= have[it];  dfs(it+1, fn, fp);}}int64 mul(int64 x, int y){int64 ask = 1;for (;y;y>>=1, x= x*x%mo)  if (y&1) ask = ask*x%mo;return ask;}int main(){int i, task;freopen("2154.in", "r", stdin);freopen("2154.out", "w", stdout);prepare_prime();scanf("%d", &task);while (task--){ans = 0;scanf("%d %d", &n, &mo);for (i = 1; i <= Div[0]; i++) Div[i] = 0; Div[0] = 0;get_prime();dfs(1, 1, 1);for (i = 1; i <= Div[0]; i++)          ans =  (ans + mul((int64)n,n/Div[i]-1)* dphi[i] % mo) % mo;printf("%d\n", (int)ans);    }return 0;}

  poj2888

  和上题的差不多，n仍然为10^9， 颜色数减到了10， 但是颜色出现了限制—— 给定（x，y）这样的二元组，表示x颜色和y颜色不能相邻！

   那么不能直接使用polya了。

   burnside是可以不管涂染颜色的限制的，但是我们得把所有的涂染方案枚出来么？这不可能。

   其实我们只需要求出burnside中每个置换的长度为1的循环数即可。

   熟悉hnoi2008 cards 的同学可以发现，我们仍然可以使用相同方法的dp来求：对于每个置换，找出它的所有循环，把循环内部涂上同一种颜色的方案数就是burnside中针对方案的置换的1长度的循环数，当然，裸dp是不现实的，不过这样的dp优化应该是矩乘乘法的经典问题了。

    那么同上，枚举n的每个因数，配合欧拉函数和矩阵乘法，运用burnside定理计算。

   这里的除法不能搞上面的猥琐了，不过mod的数是prime了！ 好爽哦，可以直接用除法的逆元了！

  code (速度好慢~~~)

# include <cstdlib># include <cstdio># include <cmath># include <cstring>using namespace std;const int mo = 9973, sqN = 100000+100, M = 15;bool step[sqN];int  have[sqN], pri[sqN], phi[sqN], divi[sqN], dphi[sqN], num[sqN];int f[M][M], g[M][M], tmp[M][M];bool can[M][M];int tot, n, m, k, ans;void prepare_prime(){  int i, j;  phi[1] = 1;  for (i = 2; i <= 100000; i++)  {    if (!step[i])       pri[++tot] = i, phi[i] = i-1;    for (j = 1; j <= tot; j++)    {      if (pri[j]*i > 100000) break;      step[pri[j]*i] = true;      if (i % pri[j] == 0)      {phi[i*pri[j]] = phi[i] * pri[j];        break;         }      else phi[i*pri[j]] = phi[i] * (pri[j] - 1);    }  }}void get_prime(){  int i;  //memset(have, 0, sizeof(have));  //memset(num, 0, sizeof(num));  //memset(divi, 0, sizeof(divi));   for (i = 1; i <= have[0]; i++) have[i] = num[i] = 0;  for (i = 1; i <= divi[0]; i++) divi[i] = 0;  have[0] = divi[0] = 0;   int it = n;   for (i = 1; i <= tot && pri[i] <= n; i++)  if (n % pri[i]==0)  {    have[++have[0]] = pri[i];    while (it % pri[i] == 0) it /= pri[i], num[have[0]]++;  }  if (it != 1) have[++have[0]] = it, num[have[0]] = 1;}void dfs(int it, int now, int fai){  int tn, tf, i;  if (it == have[0]+1)  {    divi[++divi[0]] = now;    dphi[divi[0]] = fai;    return;  }  dfs(it+1, now, fai);  tn = now * have[it]; tf = fai * (have[it]-1);  dfs(it+1, tn, tf);  for (i = 2; i <= num[it]; i++)  {    tn *= have[it]; tf *= have[it];    dfs(it+1, tn, tf);  }}void mul(int c[M][M], int a[M][M], int b[M][M]){  int i, j, k;  memset(tmp, 0, sizeof(tmp));  for (i = 1; i <= m; i++)    for (j = 1; j <= m; j++)      for (k = 1; k <= m; k++)        (tmp[i][j] += a[i][k]*b[k][j]);  //memcpy(c, tmp, sizeof(tmp));  for (i = 1; i <= m; i++)    for (j = 1; j <= m; j++)      c[i][j] = tmp[i][j]%mo;}int ask(int Q){  int i, j;  for (i = 1; i <= m; i++)    for (j = 1; j <= m; j++)      f[i][j] = g[i][j] = can[i][j];  for (Q--;Q;Q>>=1, mul(g, g, g))    if (Q & 1) mul(f, f, g);  int ask = 0;  for (i = 1; i <= m; i++)    ask = (ask + f[i][i]) % mo;  return ask;}int pow(int x, int y){  int ask = 1;  for (;y;y>>=1, x=x*x% mo)   if (y & 1) ask = ask*x%mo;  return ask;}int main(){  int i, task, x, y;  freopen("2888.in", "r", stdin);  freopen("2888.out", "w", stdout);  prepare_prime();  scanf("%d", &task);  while (task--)  {     ans = 0;     scanf("%d%d%d", &n, &m, &k);     memset(can, true, sizeof(can));     for (i = 1; i <= k; i++)       scanf("%d%d", &x, &y), can[x][y]=can[y][x]=0;     get_prime();      dfs(1, 1, 1);      for (i = 1; i <= divi[0]; i++)       ans = (ans + dphi[i]%mo * ask(n/divi[i])) % mo;       ans = ans *pow(n% mo, mo-2) % mo;     printf("%d\n", ans);   }   return 0;}