根据polya定理，答案应该是
[frac{1}{n}sum_{i=1}^{n}n^{gcd(i,n)}
]
但是这个显然不能直接求，因为n是1e9级别的，所以推一波式子：
[frac{1}{n}sum_{i=1}^{n}n^{gcd(i,n)}
]
[=frac{1}{n}sum_{d|n}n^dsum_{i=1}^{n}[gcd(i,n)==d]
]
[=frac{1}{n}sum_{d|n}n^dsum_{i=1}^{frac{d}{n}}[gcd(i,frac{d}{n})==1]
]
[=frac{1}{n}sum_{d|n}n^dvarphi (frac{d}{n})
]
[=sum_{d|n}n^{d-1}varphi (frac{d}{n})
]
这样就可以求了，但是注意时间还是很紧，所以开long long会T，求phi不预处理质数也会T
#include<iostream>
#include<cstdio>
using namespace std;
const int N=100005;
int T,n,mod,ans,p[N],tot;
bool v[N];
int ksm(int a,int b)
{
int r=1;
a%=mod;
while(b)
{
if(b&1)
r=r*a%mod;
a=a*a%mod;
b>>=1;
}
return r;
}
int phi(int n)
{
int r=n;
for(int i=0;p[i]*p[i]<=n;i++)
if(n%p[i]==0)
{
r=r-r/p[i];
while(n%p[i]==0)
n/=p[i];
}
if(n>1)
r=r-r/n;
return r%mod;
}
int main()
{
for(int i=2;i<=100000;i++)
if(!v[i])
{
p[tot++]=i;
for(int j=i+i;j<100000;j+=i)
v[j]=1;
}
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&n,&mod);
ans=0;
for(int i=1;i*i<=n;i++)
if(n%i==0)
{
ans=(ans+ksm(n,i-1)*phi(n/i))%mod;
if(i*i!=n)
ans=(ans+ksm(n,n/i-1)*phi(i))%mod;
}
printf("%d
",ans);
}
return 0;
}