bzoj3113/hdu2481 Toy（矩阵树定理+polya+矩阵快速幂）

这神题，真是给跪了orz
首先你需要做过bzoj1002，我们得到了一个递推式$f[n]=3f[n-1]-f[n-2]+2$，用来算n个点的方案数。
然后这题旋转同构，你需要polya，因为n很大，所以要欧拉函数优化，所以答案就是$\sum\limits_{i|n}\phi(i)f[n/i]$。你需要$O(\sqrt n）$枚举因子d，$O(\sqrt d)$算欧拉函数，用矩阵快速幂$O(3^3logd)$算f。

然而最后还要除以n，n,m居然不保证互质，不一定存在逆元！怎么办呢qaq
我们有$\frac{a}{b}\ mod\ p=\frac{a\ mod\ bp}{b}$，当$b|a$时成立。
因此我们就一直%nm就好了，最后可以直接除以n。
问题又来了，这样乘法就会爆ll了！那我们可以写快速乘。又带上了一个log qaq
这样就可以通过hdu啦！
但是bzoj只给了1s…
怎么办呢，我们还有神奇的O(1)快速乘！orz orz Orz

然后这题就完了呢！owo

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define ll long long
#define ld long double
#define inf 0x3f3f3f3f
#define N 100010
inline int read(){
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    return x*f;
}
int n,prime[N],tot=0,phi[N];ll mod,ans;
bool notprime[N];
/*inline ll mul(ll a,ll b){
    ll res=0,f=1;if(a<0) a=-a,f=-1;
    for(;a;a>>=1,b=(b+b)%mod) if(a&1) (res+=b)%=mod;return res*f;
}*/
inline ll mul(ll a,ll b){
    ll tmp=a*b-(ll)((ld)a/mod*b+1e-8)*mod;return tmp<0?tmp+mod:tmp;
}
struct Matrix{
    ll a[3][3];
    Matrix(bool t){memset(a,0,sizeof(a));if(t) for(int i=0;i<3;++i) a[i][i]=1;}
    friend Matrix operator*(Matrix a,Matrix b){
        Matrix res(0);
        for(int i=0;i<3;++i)
            for(int j=0;j<3;++j)
                for(int k=0;k<3;++k)
                    (res.a[i][j]+=mul(a.a[i][k],b.a[k][j]))%=mod;return res;
    }friend Matrix operator^(Matrix a,int k){
        Matrix res(1);for(;k;k>>=1,a=a*a) if(k&1) res=res*a;return res;
    }
}a(0),trans(0);
inline void getprime(){
    notprime[1]=1;phi[1]=1;
    for(int i=2;i<=100000;++i){
        if(!notprime[i]) prime[++tot]=i,phi[i]=i-1;
        for(int j=1;prime[j]*i<=100000;++j){
            notprime[prime[j]*i]=1;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
}
inline int ph(int x){
    if(x<=100000) return phi[x]%mod;
    int res=x,xx=x;
    for(int j=1;prime[j]*prime[j]<=x&&xx!=1;++j){
        if(xx%prime[j]) continue;
        res-=res/prime[j];while(xx%prime[j]==0) xx/=prime[j];
    }if(xx!=1) res-=res/xx;return res%mod;
}
inline ll calc(int x){
    if(x==1) return 1;if(x==2) return 5;
    Matrix res=a*(trans^(x-2));return res.a[0][0];
}
int main(){
//  freopen("a.in","r",stdin);
    getprime();a.a[0][0]=5;a.a[0][1]=1;a.a[0][2]=2;trans.a[0][0]=3;trans.a[1][0]=-1;
    trans.a[2][0]=1;trans.a[0][1]=1;trans.a[2][2]=1;
    while(~scanf("%d",&n)){
        mod=read();mod*=n;ans=0;int x;
        for(x=1;x*x<n;++x)
            if(n%x==0) (ans+=mul(ph(n/x),calc(x))+mul(ph(x),calc(n/x)))%=mod;
        if(x*x==n) (ans+=mul(ph(x),calc(x)))%=mod;mod/=n;
        printf("%lld\n",(ans/n+mod)%mod);
    }return 0;
}