hdu2481 Toy

这题就是【bzoj1002 轮状病毒】和【poj2888 Magic Bracelet】的综合版。用前者找到的规律和后者求不动点的方法就可以了。
需要注意的是$m$不一定是质数，需要利用$a/b\% m=a\%(b*m)/b$在模$n*m$下计算。但是这样乘法就会爆long long了，需要快速乘。

#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();
}