[BZOJ3202][SDOI2013]项链

bzoj
luogu

sol

数论多合一。
可以考虑分两步做。
1、求有多少种不同的珠子，设其数量为(m)。
2、用(m)种颜色对(n)个珠子的项链染色，要求相邻不同色，且旋转相同视作同一种方案，求方案数。

求有多少种不同的珠子

(A=sum_{i=1}^asum_{j=1}^asum_{k=1}^a[gcd(i,j,k)=1])？
然而算重了。
如果(i,j,k)都不相等的话就多算了(6)次，如果有两个相等就多算了(3)次，如果三个都相等......那肯定是三个(1)啦，并没有多算。
这个问题怎么解决呢？我们再算(B=sum_{i=1}^asum_{j=1}^a[gcd(i,j)=1])，那么答案就是(frac 16(A+3B+2))（把每种情况补齐到(6)次，再一起除以(6)）。
求(A)跟(B)是很简单的莫比乌斯反演，具体来说就是：

[A=sum_{i=1}^alfloorfrac ai floor^3mu(i)\B=sum_{i=1}^alfloorfrac ai floor^2mu(i) ]

求染色方案

这里要用到(Burnside)定理。
这里的置换只有一种，就是旋转。设旋转(i)个位置，那么就会有产生(gcd(i,n))个等价类。
所以答案就是$$frac{sum_{i=1}^nf_{gcd(i,n)}}{n}$$其中(f_i)表示(i)个点的环染色的方案数。也可以认为是给一个长为(n)的序列染色，然后保证首尾不同色的方案数。
考虑对上式做一些简化。因为(gcd(i,n))一定是(n)的约数，所以可以枚举某个约数(d)，然后答案就变成了$$frac{sum_{d|n}f_dsum_{i=1}^n[gcd(i,n)=d]}{n}=frac{sum_{d|n}f_dvarphi(frac nd)}{n}$$

求(f_n)

考虑一下(f_n)的递推。

[f_n=(m-2)f_{n-1}+(m-1)f_{n-2},f_1=0,f_2=m(m-1) ]

讲道理的话，前一种转移表示在末尾填一个与首尾都不相同的颜色，这样有(m-2)种选择；后一种转移表示先填一个和首位相同的颜色，再填一个不同色，这样有(1 imes(m-1))种选择。
然后就可以矩乘了。

求(f_n)的通项公式

有人跟我讲这题要求通项公式我就没去写矩乘。常数优秀一点说不定跑得过？
下面讲一下求通项公式

[f_n-(m-2)f_{n-1}-(m-1)f_{n-2}=0 ]

其特征方程$$x^2-(m-2)x-(m-1)=0$$
解得其特征根$$x_1=-1,x_2=m-1$$
于是设其通项公式为$$f_i=a(-1)^i+b(m-1)^i$$
将(f_1=0,f_2=m(m-1))代入解得$$a=1-m,b=1$$
所以$$f_i=(1-m)(-1)^i+(m-1)^i$$
复杂度还是(O(log n))但是常数小很多。

求(varphi(n))

最好不要枚举因数然后暴力算(varphi(n))。
一种可行的方案是：先将(n)质因数分解，然后(dfs)搜出每一个因数，在搜的过程中自然可以求出它的(varphi)。
这样的话总的复杂度就是(O(a+T(sqrt a+d(n)log n)))。

你以为这样就做完了吗

(nle10^{14})
这样可能会产生一个问题：(n)是(P=10^9+7)的倍数。这样算出来的分子要除以(n)，不好意思不存在逆元。
怎么办呢？
我们把答案对(P^2)取模。如果(P|n)，那么分子既然是(n)的倍数就也是(P)的倍数。直接除以(P)，接着还需要除以(frac nP)，而(frac nP)在模(P)意义下是存在逆元的（(nle P^2)），于是就可以乘逆元了。

code

#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long
const ll mod = 1000000007ll;
const ll MOD = 1000000007ll*1000000007ll;
const ll Inv6 = 833333345000000041ll;
const int N = 10000005;
int T,zhi[N],pri[N/10],tot,mu[N],cnt,q[30];ll n,m,p[30],ans;
void add(ll &x,ll y){x+=y;if(x>=MOD)x-=MOD;}
ll mul(ll x,ll y){
	return (x*y-(ll)(((long double)x*y+0.5)/(long double)MOD)*MOD+MOD)%MOD;
}
ll sqr(ll x){return mul(x,x);}
ll cub(ll x){return mul(mul(x,x),x);}
ll fastpow(ll a,ll b){
	ll res=1;
	while (b) {if (b&1) res=mul(res,a);a=mul(a,a);b>>=1;}
	return res;
}
ll fpow(ll a,ll b){
	ll res=1;
	while (b) {if (b&1) res=res*a%mod;a=a*a%mod;b>>=1;}
	return res;
}
void init(){
	mu[1]=1;
	for (int i=2;i<N;++i){
		if (!zhi[i]) pri[++tot]=i,mu[i]=-1;
		for (int j=1;i*pri[j]<N;++j){
			zhi[i*pri[j]]=1;
			if (i%pri[j]) mu[i*pri[j]]=-mu[i];
			else break;
		}
	}
	for (int i=1;i<N;++i) mu[i]+=mu[i-1];
}
ll cal(ll n){
	ll A=0,B=0;
	for (ll i=1,j;i<=n;i=j+1){
		j=n/(n/i);
		add(A,mul(cub(n/i),mu[j]-mu[i-1]+MOD)%MOD);
		add(B,mul(sqr(n/i),mu[j]-mu[i-1]+MOD)%MOD);
	}
	add(A,mul(B,3));add(A,2);return mul(A,Inv6);
}
ll F(ll n){
	ll res=n&1?1-m:m-1;add(res,MOD);
	add(res,fastpow(m-1+MOD,n));return res;
}
void dfs(int i,ll cur,ll phi){
	if (i==cnt+1) {add(ans,mul(phi,F(n/cur)));return;}
	dfs(i+1,cur,phi);
	cur*=p[i];phi*=p[i]-1;dfs(i+1,cur,phi);
	for (int j=2;j<=q[i];++j)
		cur*=p[i],phi*=p[i],dfs(i+1,cur,phi);
}
int main(){
	init();
	scanf("%d",&T);while (T--){
		scanf("%lld%lld",&n,&m);m=cal(m);
		ll x=n;cnt=0;
		for (int i=1;i<=tot&&(ll)pri[i]*pri[i]<=x;++i)
			if (x%pri[i]==0){
				p[++cnt]=pri[i];q[cnt]=0;
				while (x%pri[i]==0) x/=pri[i],++q[cnt];
			}
		if (x>1) p[++cnt]=x,q[cnt]=1;
		ans=0;dfs(1,1,1);
		if (n%mod==0){
			ans/=mod;
			ans=ans*fpow(n/mod,mod-2)%mod;
		}else{
			ans%=mod;
			ans=ans*fpow(n%mod,mod-2)%mod;
		}
		printf("%lld
",ans);
	}
	return 0;
}