LGP3307题解
题意有点儿神秘，而且出题人可能有点大病（
-
项链由 \(n\) 颗珠子构成，相邻的珠子不能相同。
-
每颗珠子上有 \(3\) 个数字，这 \(3\) 个数之间没有顺序，且 \(\gcd\) 为 \(1\)。
为什么珠子是三菱柱啊
明确一下思路：先算出有多少种可能的珠子，再计算方案数。
有多少种可能的珠子莫反一下就好了，很简单。
然后考虑 polya。
设 \(f(n,m)\) 是长度为 \(n\) 的环且环上的数值域为 \([1,m]\)，不存在相邻的数相同的环的数量。（下面将 \(f(n,m)\) 写成 \(f(n)\)）
来考虑最后两个位置。
-
靠左的与第一个相同：\(f(n-2)\times(m-1)\)。
-
靠左的与第一个不同：\(f(n-1)\times(m-2)\)。
所以：\(f(n)=(m-2)\times f(n-1)+(m-1)\times f(n-2)\)。
注意到这玩意儿是 \([x^n]\frac{m(m-1)x^2}{1-(m-2)x-(m-1)x^2}=m(m-1)[x^{n-2}]\frac{1}{1-(m-2)x-(m-1)x^2}\)。
对于后面这玩意儿，如果你懒得算的话，可以设一个 \(p[k][m]=x^{kn}\pmod{1-(m-2)x-(m-1)x^2}\)，然后像快速幂那样乘起来。（其实就和矩阵的光速幂是一个玩意儿）
你也可以直接把通项推出来，得到 \(f(n)=(m-1)^n+(-1)^n(m-1)\)。
然后。。。如果 \(n\) 是 \(1e9+7\) 的倍数，你可就有大麻烦了。。。
需要将模数修改为 \((1e9+7)^2\)，然后写龟速乘。
#include<cstdio>
typedef long long ll;
const int SZ=1<<10|5,M=1e7+5,mod=1e9+7,MOD=mod-1;
int T,mx,top,a[15],pri[M],mu[M],Smu[M];ll n[15];bool zhi[M];
int len,ans,k[15];ll p[15];
int p1[SZ],p2[SZ],p3[SZ];
inline void sieve(const int&M){
mu[1]=1;
for(int i=2;i<=M;++i){
if(!zhi[i])pri[++top]=i,mu[i]=-1;
for(int x,j=1;j<=top&&(x=i*pri[j])<=M;++j){
zhi[x]=true;if(!(i%pri[j]))break;mu[x]=-mu[i];
}
}
for(int i=1;i<=M;++i)Smu[i]=Smu[i-1]+mu[i];
}
inline int pw2(const int&n){
return 1ll*n*n%mod;
}
inline int pw3(const int&n){
return 1ll*pw2(n)*n%mod;
}
inline int Getlim(const int&n){
int L,R,x,sum1(0),sum2(0);
for(L=1;L*L<=n;++L)sum1=(sum1+1ll*(mod+mu[L])*pw3(x=n/L))%mod;
for(;L<=n;L=R+1){
if(x*L>n)--x;R=n/x;
sum1=(sum1+1ll*(Smu[R]-Smu[L-1]+mod)*pw3(x))%mod;
}
for(L=1;L*L<=n;++L)sum2=(sum2+1ll*(mod+mu[L])*pw2(x=n/L))%mod;
for(;L<=n;L=R+1){
if(x*L>n)--x;R=n/x;
sum2=(sum2+1ll*(Smu[R]-Smu[L-1]+mod)*pw2(x))%mod;
}
return 166666668ll*(sum1+3ll*sum2+2)%mod;
}
inline int pow(const int&idx){
return 1ll*p1[idx&1023]*p2[idx>>10&1023]%mod*p3[idx>>20]%mod;
}
int kx;
inline void DFS(const ll&n,const int&id,const ll&d,const int&phi,const int&m){
if(id==len+1){
ans=(ans+1ll*(pow((n/d)%MOD)+((n/d)&1?mod+1-m:m-1))*phi)%mod;return;
}
DFS(n,id+1,d,phi,m);
ll fd=p[id];int fphi=p[id]%mod-1;
for(int i=1;i<=k[id];++i){
DFS(n,id+1,d*fd,1ll*phi*fphi%mod,m);
fd=fd*p[id];fphi=1ll*fphi*p[id]%mod;
}
}
inline void init(int x,const int&lim){
p1[0]=p2[0]=p3[0]=1;
p1[1]=x;
for(int i=1;i<=1024;++i){
if((i<<0*10)>lim)return;p1[i]=1ll*p1[i-1]*x%mod;
}
x=p1[1024];p2[1]=x;
for(int i=1;i<=1024;++i){
if((i<<1*10)>lim)return;p2[i]=1ll*p2[i-1]*x%mod;
}
x=p2[1024];p3[1]=x;
for(int i=1;i<=1024;++i){
if((i<<2*10)>lim)return;p3[i]=1ll*p3[i-1]*x%mod;
}
}
inline int INV(int n){
int b=mod-2,ans(1);for(;b;b>>=1,n=1ll*n*n%mod)if(b&1)ans=1ll*ans*n%mod;return ans;
}
inline int Solve(ll n,int m){
ll N=n;
m=Getlim(m);init(m-1,n>MOD?MOD:n);len=ans=0;
for(int i=2;1ll*i*i<=n;++i)if(!(n%i)){
p[++len]=i;k[len]=0;while(!(n%i))++k[len],n/=i;
}
if(n^1)p[++len]=n,k[len]=1;
DFS(N,1,1,1,m);
return 1ll*ans*INV(N%mod)%mod;
}
namespace HELL_MOD{
typedef long long ll;
const ll mod=1000000007ll*1000000007ll;
ll p1[SZ],p2[SZ],p3[SZ],p4[SZ],p5[SZ];ll ans;
inline ll Mod(const ll&n){
return n<0?n+mod:n>=mod?n-mod:n;
}
inline ll Add(const ll&a,const ll&b){
return a+b>=mod?a+b-mod:a+b;
}
inline ll times(ll a,ll b){
const unsigned long long x=a,y=b,MOD=mod;
return (x*y-(unsigned long long)((long double)x/MOD*y)*MOD+MOD)%mod;
if(a<b)a^=b^=a^=b;
ll ans(0);for(;b;b>>=1,a=Add(a,a))if(b&1)ans=Add(ans,a);return ans;
}
inline ll pw2(const ll&n){
return times(n,n);
}
inline ll pw3(const ll&n){
return times(pw2(n),n);
}
inline ll Getlim(const int&n){
int L,R,x;ll sum1(0),sum2(0);
for(L=1;L*L<=n;++L)sum1=Add(sum1,times(Mod(mu[L]),pw3(x=n/L)));
for(;L<=n;L=R+1){
if(x*L>n)--x;R=n/x;
sum1=Add(sum1,times(Add(Mod(Smu[R]),Mod(mod-Smu[L-1])),pw3(x)));
}
for(L=1;L*L<=n;++L)sum2=Add(sum2,times(Mod(mu[L]),pw2(x=n/L)));
for(;L<=n;L=R+1){
if(x*L>n)--x;R=n/x;
sum2=Add(sum2,times(Add(Mod(Smu[R]),Mod(mod-Smu[L-1])),pw2(x)));
}
return times(833333345000000041ll,(sum1+3ll*sum2+2)%mod);
}
inline ll pow(const ll&idx){
return times(times(times(times(p1[idx&1023],p2[idx>>10&1023]),p3[idx>>20&1023]),p4[idx>>30&1023]),p5[idx>>40]);
}
inline void DFS(const ll&n,const int&id,const ll&d,const ll&phi,const ll&m){
if(id==len+1){
ans=Add(ans,times(Add(pow(n/d),(n/d)&1?mod+1-m:m-1),phi));
return;
}
DFS(n,id+1,d,phi,m);
ll fd=p[id],fphi=p[id]%mod-1;
for(int i=1;i<=k[id];++i){
DFS(n,id+1,d*fd,phi*fphi,m);
fd=fd*p[id];fphi=fphi*p[id];
}
}
inline void init(ll x,const ll&lim){
p1[0]=p2[0]=p3[0]=p4[0]=p5[0]=1;
p1[1]=x;
for(int i=1;i<=1024;++i){
if((1ll*i<<0*10)>lim)return;p1[i]=times(p1[i-1],x);
}
x=p1[1024];p2[1]=x;
for(int i=1;i<=1024;++i){
if((1ll*i<<1*10)>lim)return;p2[i]=times(p2[i-1],x);
}
x=p2[1024];p3[1]=x;
for(int i=1;i<=1024;++i){
if((1ll*i<<2*10)>lim)return;p3[i]=times(p3[i-1],x);
}
x=p3[1024];p4[1]=x;
for(int i=1;i<=1024;++i){
if((1ll*i<<3*10)>lim)return;p4[i]=times(p4[i-1],x);
}
x=p4[1024];p5[1]=x;
for(int i=1;i<=1024;++i){
if((1ll*i<<4*10)>lim)return;p5[i]=times(p5[i-1],x);
}
}
inline int Solve(ll n,ll m){
ll N=n;
m=Getlim(m);init(m-1,n);len=ans=0;
for(int i=2;1ll*i*i<=n;++i)if(!(n%i)){
p[++len]=i;k[len]=0;while(!(n%i))++k[len],n/=i;
}
if(n^1)p[++len]=n,k[len]=1;DFS(N,1,1,1,m);
return 1ll*(ans/::mod)*INV(N/::mod)% ::mod;
}
}
signed main(){
scanf("%d",&T);
for(int i=1;i<=T;++i){
scanf("%lld%d",n+i,a+i);
if(a[i]>mx)mx=a[i];
}
sieve(mx);
for(int i=1;i<=T;++i)printf("%d\n",n[i]%mod?Solve(n[i],a[i]):HELL_MOD::Solve(n[i],a[i]));
}
