POJ2154 以欧拉函数优化的BurnSide定理（有部分详细证明）

世界真的很大

这是burnside定理刷的第二道题，还是很有营养的，证明了不同的数论定理可以结合起来用，我先来解释一下题意吧：给你一个项链，长度为n，同时还有n种颜色，考虑旋转且不考虑翻转，求有多少种染色方案？

看过我上一篇题解（poj2409）的观众朋友应该一看就会想到burnside定理，用k枚举置换数累加n^gcd(n,k),再用公式直接求解？完全正确！！而且这道题不考虑翻转，所以直接/n即置换数就行了，很简单吧？当然不会，n的范围是10^9，，，，枚举置换的话，会超时吧？怎么办呢。。

来仔细看一下定理，累加n^gcd(n,k)，我们累加的其实只是n的gcd（n，k）次方而已，这里用d表示gcd（n，k）。d一定是n的因数，没得说，那我们其实可以直接枚举这个d，对吧，原先枚举k的目的其实也只是为了这个d罢了。枚举n的因数d，就不用枚举到n了，只用枚举到根号n，因为只要枚举了因数d，那么n/d也同样是n的因数没错，那n^gcd(n,k)就可以直接用n^d来表示，并且每一个d都有对应的n/d，时间复杂度快了不少，10^9开根号，只有10^5还要少，完全没有问题。但是啊，我仅仅枚举d可以吗，会不会有多个k和n的gcd都是同一个d呐，那不是少加了？那的确是少加了，又怎么办？

在回头看一下公式的一部分吧：d=gcd（n，k），对吧，试试这样：d/d=gcd（n/d，k/d），即1=gcd（n/d，k/d），对吧，我们想要知道在d相同时，k的个数，那不就是说k/d的个数吗，再看，k/d与n/d是互质的，k严格小于n（因为n是总置换数），那就是求与n/d互质且小于n/d的数的数量，即n的欧拉函数值；

代码：

for( i=1;i*i<n;i++)
if( n%i==0 )
ans=(ans+phi(n/i)*quickmub(n,i-1)+phi(i)*quickmub(n,n/i-1))%p;

这里的phi就是欧拉函数值。

但注意如果n正好能开方为整数的话，那么还有种情况：

if(i*i == n)
ans=(ans+phi(i)*quickmub(n,i-1))%p;

因为此时d == n/d，所以只要加一次就行了。

至于为什么是n的d-1次方呐，（i代指d）因为在这里颜色数n等于置换数n，所以直接把1/n乘到n^d里，就是n^(d-1)了，嗯，就是这样。所以说最后我们就不用再除以总置换数了，因为每次累加时已经乘了。。。

还有还有，这个欧拉函数我们怎么球呢？

我开始时使用的是欧拉筛，先线性处理欧拉函数之后o（1）查询，但看到数据范围到10^9，数组保存不了那么大。。所以果断在线处理，嗯。

但首先我们还是需要用线性筛筛一遍质数，其实也不会很多，10^9里质数不会超过40万个，数组可以存，代码如下：

void init(int n)
{
isnot[1]=true;
for(int i=2;i<=n;i++)
{
if(!isnot[i])
{
primes[ptot++]=i;
for(int j=i+i;j<=n;j+=i)
{
isnot[j]=true;
}
}
}
}

在线求欧拉函数其实需要用到欧拉函数的一个性质：

phi（n）=n*（1-1/p）。。。。。其中p是n的所有不同种类的质因数，那其实就用线性筛处理好的质数来分解n就行了，代码：

int phi(int x)
{
int temp=x;
for(int i=0;i<ptot&&primes[i]*primes[i]<=x;i++)
{
if(x%primes[i]==0)
{
temp=temp/primes[i]*(primes[i]-1);
while(x%primes[i]==0)
{
x/=primes[i];
}
}
}
if(x>1) temp=temp/x*(x-1);
return temp%p;
}

这道题就是这样了，注意的地方不多，算是2409的很不错的拓展了。嗯

完整代码：

#include<stdio.h>
#include<cstring>
using namespace std;
int p;
int primes[200015],ptot=0;
bool isnot[200015];

void init(int n)
{
isnot[1]=true;
for(int i=2;i<=n;i++)
{
if(!isnot[i])
{
primes[ptot++]=i;
for(int j=i+i;j<=n;j+=i)
{
isnot[j]=true;
}
}
}
}

int phi(int x)
{
int temp=x;
for(int i=0;i<ptot&&primes[i]*primes[i]<=x;i++)
{
if(x%primes[i]==0)
{
temp=temp/primes[i]*(primes[i]-1);
while(x%primes[i]==0)
{
x/=primes[i];
}
}
}
if(x>1) temp=temp/x*(x-1);
return temp%p;
}
int quickmub(int a,int b)
{
int ret;a=a%p;
for(ret=1;b;b>>=1,a=(a*a)%p) if(b&1) ret=(ret*a)%p;
return ret;
}

int main()
{
int t;
scanf("%d",&t);
init(200010);
for(int e=1;e<=t;e++)
{
int n;
scanf("%d%d",&n,&p);
int ans=0;
int i;
for( i=1;i*i<n;i++)
if( n%i==0 )
ans=(ans+phi(n/i)*quickmub(n,i-1)+phi(i)*quickmub(n,n/i-1))%p;
if(i*i == n)
ans=(ans+phi(i)*quickmub(n,i-1))%p;
printf("%d\n",ans%p);

}
return 0;
}