}
return ans;
}
//以上为矩阵快速幂乘
void Prime(){
for(int i=2;i<=sqrt(N+1.0);i++){
if(flag[i]) continue;
prime[cnt++]=i;
for(int j=2;j*i<=sqrt(N+1.0);j++)
flag[i*j]=true;
}
}
int Eular(int n){
int ret=1;
for(int i=0;i
n/=prime[i];ret*=prime[i]-1;
while(n%prime[i]==0){n/=prime[i];ret=(ret*prime[i])%MOD;}
}
}
if(n>1) ret*=n-1;
return ret%MOD;
}
//以上为打素数表，求欧拉函数
LL Get_T(int k){
if(k==1) return 1;
else if(k==2) return 5;
Matrix temp=init^(k-2);
LL g=2*(f-(3*temp.m[0][1]+temp.m[1][1])-1);
return (g+f)%MOD;
}
//计算T值
LL Polya(){
LL sum=0;
int i;
//Burnside定理，枚举循环个数
for(i=1;i*i
sum=(sum+MultMod(Eular(i),Get_T(n/i)))%MOD;
sum=(sum+MultMod(Eular(n/i),Get_T(i)))%MOD;
}
if(i*i==n) sum=(sum+MultMod(Get_T(i),Eular(i)))%MOD;
return sum/n;
}
int main(){
Prime(); www.2cto.com
//构造矩阵
init.m[0][0]=3;init.m[0][1]=1;init.m[1][0]=-1;init.m[1][1]=0;
while(scanf("%d%I64d",&n,&MOD)!=EOF){
MOD=(LL)n*MOD;
printf("%I64d\n",Polya()%(MOD/n));
}
return 0;
}
作者：ACM_cxlove