hdu2481 Toy

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LL long long
const int maxx=100000;
int n,m,tot,have[maxx+10],prm[maxx+10];
LL p;
int phi(int x)
{
int ret=x;
for (int i=1;i<=tot&&(LL)prm[i]*prm[i]<=x;i++)
if (x%prm[i]==0)
{
ret=ret/prm[i]*(prm[i]-1);
while (x%prm[i]==0) x/=prm[i];
}
if (x>1) ret=ret/x*(x-1);
return ret%p;
}
LL multi(LL base,LL x)
{
LL ret=0;
while (x)
{
if (x&1)
{
ret+=base;
if (ret>p) ret-=p;
}
base+=base;
if (base>p) base-=p;
x>>=1;
}
return ret;
}
struct mat
{
LL a[5][5];
mat operator * (const mat &mm) const
{
mat ret;
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
{
ret.a[i][j]=0;
for (int k=1;k<=3;k++)
ret.a[i][j]=(ret.a[i][j]+multi(a[i][k],mm.a[k][j]))%p;
}
return ret;
}
mat pow(int k)
{
mat ret,base;
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
{
base.a[i][j]=a[i][j];
ret.a[i][j]=(i==j);
}
for (;k;k>>=1,base=base*base)
if (k&1) ret=ret*base;
return ret;
}
}ori,trans,res;
LL calc(int n)
{
res=ori*trans.pow(n-1);
return res.a[1][1];
}
void solve()
{
LL ans=0;
p=(LL)m*n;
trans.a[2][1]=p-1;
for (int i=1;(LL)i*i<=n;i++)
if (n%i==0)
{
ans=(ans+multi(phi(n/i),calc(i)))%p;
if (i*i<n) ans=(ans+multi(phi(i),calc(n/i)))%p;
}
printf("%d\n",ans/n);
}
int main()
{
ori.a[1][1]=trans.a[1][2]=trans.a[3][1]=trans.a[3][3]=1;
ori.a[1][3]=2;
trans.a[1][1]=3;
for (int i=2;i<=maxx;i++)
{
if (!have[i]) prm[++tot]=i;
for (int j=1;(LL)i*prm[j]<=maxx;j++)
{
have[i*prm[j]]=1;
if (i%prm[j]==0) break;
}
}
while (scanf("%d%d",&n,&m)==2) solve();
}