这神题,真是给跪了orz
首先你需要做过bzoj1002,我们得到了一个递推式
f[n]=3f[n−1]−f[n−2]+2
,用来算n个点的方案数。
然后这题旋转同构,你需要polya,因为n很大,所以要欧拉函数优化,所以答案就是
∑i|nϕ(i)f[n/i]
。你需要
O(n√)
枚举因子d,
O(d√)
算欧拉函数,用矩阵快速幂
O(33logd)
算f。
然而最后还要除以n,n,m居然不保证互质,不一定存在逆元!怎么办呢qaq
我们有
ab mod p=a mod bpb
,当
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;
}