poj 2154 Color (组合数学 polya计数法)

大致题意：给你n种颜色的珠子，和一个长度也是n((n<=1000000000)的项链，用这n种珠子串成这个项链，项链可以旋转，经过旋转所得的项链视为同一种项链，现在告诉你颜色总数和项链的长度n，由于答案可能超出int范围，题目还会给你一个取模的数p，求共能组成几条不同的项链（取模后）。

Time Limit: 2000MS　　Memory Limit: 65536K

样例：

Sample Input
5
1 30000
2 30000
3 30000
4 30000
5 30000

Sample Output
1
3
11
70
629

题目一看也是一道计数法的题目，应用polya定理应该能够解决，置换群中的置换元素个数和每个元素的循环节个数求法已经在前面的博文中详细给出了分析，不再赘述。这道题看似比poj 2409 Let it Bead简单，因为它所在的置换群中只有旋转置换，少了翻转置换，但如果采取上题中所给出的枚举i的方法就会超时，因为n现在的上限很大，10^9，我们分析下时间复杂度，看下面一段程序：

1 ll polya(){
2     int i,j;
3     ll ans=0;
4     for(i=0;i<s;i++){
5         ans+=(ll)pow(1.0*c,gcd(s,i));
6     }
7     return ans/s;
8 }

如果我们忽略求pow()和gcd()的时间，运行时间也得是O(n)，n是10^9，所以2s内是过不了的。

其实我们根本不用枚举i的值，因为我们知道gcd(n,i)的值一定是s的因子，所以它是有限多的，如果我们求出有多少数(A0~Ax)使得gcd(n,Ai)相等，那就可以大大降低复杂度。枚举gcd(n,i)变成了枚举因子，相同gcd下的i的个数也不难求，因为n/gcd(n,i)就跟i是互质的了，所以这个个数也就是n/gcd(n,i)内和n/gcd(n,i)互质的正整数个数（保证i和n不再有其他因子），欧拉函数在这就派上用场了。我们再分析下枚举gcd(n,i)的时间复杂度，（复杂度不是很好分析，可能说的不对），n分解质因数最多20多个质数因子(2^31是int范围)，最坏情况应该是每个质因子都只有一个，利用求因子公式，那么n的因子大约有(2^20)多个，化成十进制是10^7，这应该是最坏的情况了，在2s内应该还是能跑完的。

找因子的方法是dfs搜索，下面dfs的代码是参考某位大神的，因为写的实在太精髓了，强势沿用。求pow就用exp_mod模版就行，分解质因数，求欧拉函数都有模版，分解质因数和求欧拉函数可以先打素数表，这样求会快一些。至此，这个题也就能A掉了。

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cstring>
 4 #include<iostream>
 5 #include<cmath>
 6 using namespace std;
 7 bool is[50010];
 8 int prime[6005];
 9 int k,n,mod,ans;
10 void getPrime(){　　　//素数筛
11     k=0;
12     memset(is,0,sizeof(is));
13     for(int i=2;i<=50005;i++)if(!is[i]){ 
14         prime[++k]=i; 
15         for(int j=2;j*i<=50005;j++){  
16             is[i*j]=1;  
17         }  
18     }  
19 }
20 int d[30],c[30],top;
21 int exp_mod(int a,int b,int n){　　//快速幂模版
22     int ans=1;
23     a=a%n;
24     while(b){
25         if(b&1){
26             ans=(a*ans)%n;
27         }
28         b>>=1;
29         a=(a*a)%n;
30     }
31     return ans;
32 }
33 void divide(int n) {　　//分解质因数
34     top = 0;
35     for (int i = 1; i <= k; i++) {
36         int t = 0;
37         while (n % prime[i] == 0) t++, n /= prime[i];
38         if (t >= 1) {
39             d[++top] = prime[i];
40             c[top] = t;
41         }
42         if (n == 1) break;
43         if (n / prime[i] < prime[i]) {
44             d[++top] = n;
45             c[top] = 1;
46             break;
47         }
48     }
49 }
50 int euler(int n) {　　//求欧拉函数
51     int ret = n;
52     for (int i = 1; i <=k; i++) {
53         if (n % prime[i] == 0) {
54             ret -= ret / prime[i];
55             while (n % prime[i] == 0) n /= prime[i];
56         }
57         if (n == 1) break;
58         if (n / prime[i] < prime[i]) {
59             ret -= ret / n;
60             break;
61         }
62     }
63     return ret % mod;
64 }
65 void dfs(int step,int now){　　//枚举gcd(n,i)的dfs，本题的关键部分
66     int i,t;
67     if(step>=top+1){
68         ans=(ans+euler(n/now)*exp_mod(n%mod,now-1,mod)%mod)%mod;　　//这里由于最后要除以n索性直接求n^now-1，当然求n^now，最后在乘逆元也可以
69         return;　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 
70     }
71     for(i=0,t=1;i<=c[step];t*=d[step],i++){　　//这里注意细节gcd(n,1)和gcd(n,n)的情况是一样的求其中一个就行
72         dfs(step+1,now*t);
73     }
74 }
75 int main(){
76     int cas;
77     getPrime();
78     scanf("%d",&cas);
79     while(cas--){
80         scanf("%d%d",&n,&mod);
81         divide(n);
82         ans=0;
83         dfs(1,1);
84         printf("%d\n",ans);
85     }
86     return 0;
87 }