[poj2154]color

题面：

Beads of N colors are connected together into a circular necklace of N beads (N<=1000000000). Your job is to calculate how many different kinds of the necklace can be produced. You should know that the necklace might not use up all the N colors, and the repetitions that are produced by rotation around the center of the circular necklace are all neglected.

You only need to output the answer module a given number P.

INPUT：

The first line of the input is an integer X (X <= 3500) representing the number of test cases. The following X lines each contains two numbers N and P (1 <= N <= 1000000000, 1 <= P <= 30000), representing a test case.

大意：给你X个项链，每个项链有n个珠子，用n种颜色染（颜色可以不全用），求本质不同的项链方案，当且仅当一个项链任意旋转(只有旋转没有翻转)都与其他项链不同是

一开始没有仔细看题面，原来没有翻转，只有旋转

SOLUTION：如果没有(N<=10^9)，这道题就是一道裸题，直接根据polya定理计算即可，但是这道题N很大，不能直接计算，我们需要优化：

(ans=sumlimits_{i=0}^{n-1}n^{gcd(i,n)})

设(d=gcd(i,n)),(d)一定是(n)的约数，我们可以先O((sqrt{n}))处理出(n)的约数，计算每个d贡献多少次，也就是有多少个i，(gcd(i,n)=d ext{转化为}gcd(frac{i}{d},frac{n}{d})=1)，即(varphi(frac{n}{d})个),我们可以用欧拉函数

注意我们在用polya/burnside定理时最后要除以置换集的大小，以为没有翻转，我们只用除以n（不用2n），但是是在(%p)意义下，p不一定是质数所以可能没有逆元，我们直接将公式改成(varphi(n/d) imes n^{d-1})

code:

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
using namespace std;
int t,n,p,tot,cnt;
int prime[1000000],v[1000000];
inline void pre(){
 for(int i=2;i<=100000;++i){
  if(!v[i]){
   prime[++tot]=i;
   v[i]=i;
  }
  for(int j=1;j<=tot;++j){
   if(prime[j]*i>100000||prime[j]>v[i])break;
   v[i*prime[j]]=prime[j];
  }
 }
}
inline int poww(int a,int b){
     int ans=1;
     while(b){
      if(b&1)ans=ans*a%p;
       a=a*a%p;
       b>>=1;
     }
     return ans%p;
}
inline int euler(int x){
 int ans=x;
 for(int i=1;i<=tot&&prime[i]*prime[i]<=x;++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;
} 
int main(){
 pre();
 cin>>t;
 for(int i=1;i<=t;++i){
       cin>>n>>p;
       int c=0;
       for(int j=1;j<=sqrt(n);++j){
         if(n%j==0){
          c=(c+euler(n/j)%p*poww(n%p,j-1)%p)%p;
          if(j*j!=n)c=(c+euler(j)%p*poww(n%p,n/j-1)%p)%p;
         }
       }
       cout<<c%p<<endl;
 }
 return 0;
}