好题！！！没话说……
用到的知识面很多，这题的难点在于公式的推导。
原始推导过程见：http://hi.baidu.com/spellbreaker/item/d8bb3bda5af30be6795daa93
这个过程有点小问题，改进后的见：http://blog.csdn.net/acm_cxlove/article/details/7868589
下面是自己敲的代码：
/*************************************************************************
> File Name: xh.cpp
> Author: XINHUA
> Mail: 525799145@qq.com
> Created Time: 2013/7/19 星期五 15:11:25 新华
************************************************************************/ #include<iostream>
#include<stdio.h>
#include<algorithm>
#include<iomanip>
#include<cmath>
#include<string>
#include<vector>
using namespace std;
int prime[],m,n;
bool f[]={};
__int64 mod;
struct ma
{
__int64 a[][];
void init()
{
a[][]=a[][]=;
a[][]=a[][]=;
}
}e;
//使用二分模拟乘法计算（a和b的范围太大）
__int64 mulmod(__int64 a,__int64 b)
{
a%=mod;b%=mod;
if(a<) a+=mod;
if(b<) b+=mod;
__int64 ans=;
while(b)
{
if(b&)
{
ans+=a;
if(ans>=mod) ans-=mod;
}
a=a<<;
if(a>=mod) a-=mod;
b=b>>;
}
return ans%mod;
}
ma operator*(ma ab,ma b)
{
ma ans;
int i,j,k;
for(i=;i<;i++)
{
for(j=;j<;j++)
{
ans.a[i][j]=;
for(k=;k<;k++)
ans.a[i][j]=(ans.a[i][j]+mulmod(ab.a[i][k],b.a[k][j]))%mod;
}
}
return ans;
}
ma operator^(ma ab,int b)
{
ma ans;
ans.init();
while(b)
{
if(b&) ans=ans*ab;
b>>=;
ab=ab*ab;
}
return ans;
} void init_prime()
{
int i,j;
m=;
for(i=;i<=;i++)
{
if(f[i]==) prime[m++]=i;
for(j=;j<m&&i*prime[j]<=;j++)
{
f[i*prime[j]]=;
if(i%prime[j]==) break;
}
}
}
__int64 euler(__int64 n)
{
__int64 ans=;
for(int i=;i<m&&prime[i]*prime[i]<=n;i++)
{
if(n%prime[i]==)
{
ans*=prime[i]-;
n/=prime[i];
while(n%prime[i]==)
{
n/=prime[i];
ans=(ans*prime[i])%mod;
}
}
}
if(n>) ans*=n-;
return ans%mod;
}
__int64 get_T(int k)
{
if(k==) return ;
else if(k==) return ;
ma temp=e^(k-);
__int64 f=*temp.a[][]+temp.a[][];
__int64 g=*(f-(*temp.a[][]+temp.a[][])-);
return (f+g)%mod;
}
int main()
{
init_prime();
e.a[][]=;e.a[][]=;e.a[][]=-;e.a[][]=;
while(scanf("%d%I64d",&n,&mod)!=EOF)
{
mod=(__int64)n*mod;
__int64 sum=;
for(int i=;i*i<=n;i++)
{
if(n%i==)
{
sum=(sum+mulmod(euler(i),get_T(n/i)))%mod;
if(i*i!=n)
sum=(sum+mulmod(euler(n/i),get_T(i)))%mod;
}
}
sum/=n;
printf("%I64d\n",sum%(mod/n));
}
return ;
}