【BZOJ3202】项链（SDOI2013）-数论+组合数学综合

测试地址：项链
题目大意： 用$n(\le 10^{14})$个珠子串成一个环形项链，每个珠子是正三棱柱，每一面上写一个$1$ ~ $a(\le 10^7)$中的数字，当且仅当写的三个数字$\gcd$为$1$时，这个珠子是合法的珠子。两个珠子旋转或翻转后相同就看做相同。用$n$个合法的珠子串成一个环形项链，要求相邻两个珠子不同，求本质不同的方案数，两个项链经旋转后相同就看做相同。
做法： 本题是一道数论+组合数学集大成题，需要用到的比较重要的知识有：Burnside引理，莫比乌斯反演，质因数分解及大数因数枚举，欧拉函数，数列通项公式推导。
很明显，我们应该先算出珠子的种类数$m$。用Burnside引理分析，置换群中应该有$6$个置换，即置换到$1$ ~ $3$的全排列，现在分析每个置换的不动点数目（即经过置换后不变的方案数）。
首先，有一个置换是不动的，那么不动点数应该是${\sum }_{i=1}^a{\sum }_{j=1}^a{\sum }_{k=1}^a[\gcd (i,j,k)=1]$；
接着，有两个置换中，是长度为$3$的循环，那么当且仅当三个数都相等，且$\gcd$等于$1$时满足条件，因此显然对每种置换只有$1$个不动点（ { 1 , 1 , 1 } \{1,1,1\} {1,1,1}），加起来就是$2\cdot 1$。
最后，有三个置换中，是一个长度为$2$的循环加一个点，那么那个循环中的两个点应该为同一个数，两个相同的数的$\gcd$就是本身，所以不动点数实际上就是两个数$\gcd$为$1$的方案数，即${\sum }_{i=1}^a{\sum }_{j=1}^a[gcd(i,j)=1]$，加起来就是上面那个式子乘$3$。
那么合法的珠子种类数就是上面那些东西加起来再乘上$\frac{1}{6}$了。而上面的两个和式，都可以用基础的莫比乌斯反演知识推出来一个可以数论分块的式子，只需要先线性筛求出莫比乌斯函数$\mu$的前缀和即可。
现在有了$m$，就可以讨论第二个问题了：如何求本质不同的项链数。还是用Burnside引理分析，我们发现一个项链旋转$i$个单位后，置换中形成循环节长度为$\frac{n}{\gcd (n,i)}$的很多循环。具体到环上，我们发现这要求的就是，用$m$种颜色染每个点，相邻两个点颜色不同，且染色序列形成长为$\gcd (n,i)$的循环节的方案数。因为是循环，所以上面要求的就是用$m$中颜色染一个长为$\gcd (n,i)$的环，相邻两个点颜色不同的方案数，我们把这记作$f(\gcd (n,i))$。
求$f(n)$是一个经典的组合数学问题，我还在机房大佬的《高考数学你掌握了吗？》书里看到了这个题…证法比较神，虽然证明很经典，传播范围很广，但我在这里还是写一下，不想看的可以直接跳过。
我们来看，如果用$m$种颜色染一条长度为$n$的链，方案数显然为$m(m-1{)}^{n-1}$。现在要拿出一个神奇的结论：$m(m-1{)}^{n-1}=f(n)+f(n-1)$。为什么呢？如果令链的两端颜色不同，那么显然可以接起来成为环，方案数为$f(n)$，否则如果令链的两端颜色相同，那么可以把两端缩成一个点再接成环，方案数为$f(n-1)$。因此我们得到递推公式：
$f(n)=m(m-1{)}^{n-1}-f(n-1)$
两边同除$(m-1{)}^n$得：
$\frac{f(n)}{(m-1{)}^n}=\frac{m}{m-1}-\frac{1}{m-1}\cdot \frac{f(n-1)}{(m-1{)}^{n-1}}$
令$g(n)=\frac{f(n)}{(m-1{)}^n}$，有：
$g(n)=-\frac{1}{m-1}\cdot g(n-1)+\frac{m}{m-1}$
解这样的数列的通项公式就是高考难度的了（虽然还是让不懂的我非常难过），即如果有：
$g(n)=A\cdot g(n-1)+B$
那么令$h(n)=g(n)+x$，构造$h$为等比数列，则有：
$(g(n)+x)=q(g(n-1)+x)$
将这个式子和上面的式子对照，可得当$x=\frac{B}{A-1}$时，$h$就为等比数列。
在上面那个式子中算出$x=-1$，因此$h(n)=g(n)-1$是首项为$-1$，公比为$-\frac{1}{m-1}$的等比数列，那么算出$h(n)$的通项为：
$h(n)=(-1{)}^n\cdot \frac{1}{(m-1{)}^{n-1}}$
于是$g(n)=h(n)+1=(-1{)}^n\cdot \frac{1}{(m-1{)}^{n-1}}+1$，因此$f(n)=g(n)\cdot (m-1{)}^n=(m-1{)}^n+(-1{)}^n\cdot (m-1)$。
于是最后我们得到了$f(n)=(m-1{)}^n+(-1{)}^n\cdot (m-1)$这个重要的结论。最后我们只要想办法求出答案即可：
$ans=\frac{1}{n}{\sum }_{i=1}^{n}f(\gcd (n,i))$
$n$很大，不能直接计算。根据数论直觉，我们想到应该枚举$d=\gcd (n,i)$，那么使得$n$与$i$的$\gcd$为$d$的$i$的数量有多少个呢？两边除$d$，其实要求的就是与$\frac{n}{d}$互质的比它小的数有多少个，这显然等于欧拉函数$\phi (\frac{n}{d})$。因此上式可以化为：
$ans=\frac{1}{n}{\sum }_{d|n}\phi (\frac{n}{d})\cdot f(d)$
$n$的因数的个数就比$n$小很多了，可以通过质因数分解+搜索的方式枚举出来，实测不会爆。而因为$n$的质因数分解中至多有一个$>\sqrt{n}$的质因数，所以我们只需要处理出$\sqrt{n}$内的质数表即可，这可以在一开始求莫比乌斯函数的线性筛中同时求出。至于欧拉函数，也可以在枚举因数的同时求出。
于是经过漫长的讨论，这一题总算是解决了…吗？注意到一个问题，$n$很大，以至于它有可能是模数$10^9+7$的倍数，那么我们就不能直接对$n$求逆元了。解决的办法是，将上述算法中的模数改为$(10^9+7{)}^2$，那么如果$n$是$10^9+7$的倍数，则上面算出的答案对这个新模数取模后应该也是$10^9+7$的倍数，因此把$n$和算出的答案同时除$mod$，再对$n$求关于$mod$的逆元即可。而对于其他情况，直接将算出的答案再对$mod$取模进行计算就行了。至于会乘爆的问题，建议使用一种叫$O(1)$快速乘的奇葩操作，可以保证两个long long范围内相乘取模的结果不错，可以参看我代码中的mult()函数。这样我们就解决了这一题。
我傻逼的地方：在调用mult()函数前，如果两个数没有先取模，还是有可能乘爆…这样就炸了一个点。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1000000007;
const ll MOD=mod*mod;
const ll inv6=(5ll*MOD+1ll)/6ll;
int T;
ll maxn=0,n[15],x[15],a[15],m,mu[10000010],prime[10000010],totans;
vector<int> p[15];
bool vis[10000010]={0};

ll add(ll x,ll y,ll mod)
{
	return ((x+y)%mod+mod)%mod;
}

ll mult(ll x,ll y,ll mod)
{
	x%=mod,y%=mod;
	return (x*y-(ll)(((long double)x*y+0.5)/(long double)mod)*mod+mod)%mod;
}

ll power(ll a,ll b,ll mod)
{
	ll s=1,ss=a;
	while(b)
	{
		if (b&1) s=mult(s,ss,mod);
		ss=mult(ss,ss,mod);
		b>>=1;
	}
	return s;
}

void calc_prime(int limit)
{
	mu[1]=1;
	for(int i=2;i<=limit;i++)
	{
		if (!vis[i])
		{
			mu[i]=-1;
			prime[++prime[0]]=i;
			for(int j=1;j<=T;j++)
				if (n[j]%i==0)
				{
					p[j].push_back(i);
					while(x[j]%i==0) x[j]/=i;
				}
		}
		for(int j=1;j<=prime[0]&&i*prime[j]<=limit;j++)
		{
			vis[i*prime[j]]=1;
			if (i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=2;i<=limit;i++)
		mu[i]=(mu[i]+mu[i-1]+MOD)%MOD;
}

ll solve(ll n,bool type)
{
	ll ans=0;
	for(ll i=n;i;i=n/(n/i+1ll))
	{
		int l=n/(n/i+1ll)+1ll,r=i;
		if (type) ans=add(ans,mult((mu[r]-mu[l-1]+MOD)%MOD,mult(n/i,mult(n/i,n/i,MOD),MOD),MOD),MOD);
		else ans=add(ans,mult((mu[r]-mu[l-1]+MOD)%MOD,mult(n/i,n/i,MOD),MOD),MOD);
	}
	return ans;
}

ll f(ll n)
{
	ll ans=power((m-1ll+MOD)%MOD,n,MOD);
	ans=add(ans,mult((n%2)?(MOD-1ll):1ll,(m-1ll+MOD)%MOD,MOD),MOD);
	return ans;
}

ll dfssolve(ll now,ll nowphi,int id)
{
	return mult(nowphi,f(n[id]/now),MOD);
}

void dfs(ll now,ll nowphi,int id,int step,int num)
{
	if (step>=p[id].size()) return;
	if (step==0&&num==0)
		totans=add(totans,dfssolve(now,nowphi,id),MOD);
	if (num) totans=add(totans,dfssolve(now,nowphi,id),MOD);
	if (n[id]%(now*p[id][step])==0)
	{
		if (!num) dfs(now*p[id][step],nowphi*(p[id][step]-1ll),id,step,num+1);
		else dfs(now*p[id][step],nowphi*p[id][step],id,step,num+1);
	}
	dfs(now,nowphi,id,step+1,0);
}

int main()
{
	scanf("%d",&T);
	for(int i=1;i<=T;i++)
	{
		scanf("%lld%lld",&n[i],&a[i]);
		maxn=max(maxn,(ll)sqrt(n[i])+1ll);
		maxn=max(maxn,a[i]);
		x[i]=n[i];
	}
	
	calc_prime(maxn);
	for(int i=1;i<=T;i++)
	{
		m=mult(add(solve(a[i],1),mult(solve(a[i],0),3ll,MOD)+2ll,MOD),inv6,MOD);
		
		if (x[i]>1) p[i].push_back(x[i]);
		
		totans=0;
		dfs(1ll,1ll,i,0,0);
		if (n[i]%mod==0) printf("%lld\n",totans/mod*power(n[i]/mod,mod-2,mod)%mod);
		else printf("%lld\n",totans%mod*power(n[i],mod-2,mod)%mod);
	}
	
	return 0;
}