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]));

}

巴特西

LGP3307题解

最新文章

热门文章