SDOI 2013 项链 题解

题目传送门

题目大意： 一串合法的项链满足：每颗珠子是一个无序三元组，且三个数都在 $[1,M]$ 范围内，它们的 $\gcd =1$；相邻的珠子一定不同。另外，如果一串项链能够通过旋转得到另一串项链，那么认为他们是相同的。问有多少种不同的项链。

题解

容易发现其实题目分为几乎不相关的两部分：求珠子数量；求项链数量。

第一部分

无序不太好求，转化为求有序三元组数量然后去重，即：三个数互不相同$/6+$两个数相同$/3+$三个数都相同，容斥一下变成：三元组方案数$/6+$二元组方案数$\times \frac{1}{2}+$一元组方案数$\times \frac{1}{3}$，即 $($三元组方案数$+$二元组方案数$\times 3+$一元组方案数$\times 2)/6$。

三元组方案数为 ${\sum }_{i=1}^M{\sum }_{j=1}^M{\sum }_{k=1}^M[\gcd (i,j,k)=1]$，根据莫反套路转成 ${\sum }_{d=1}^M\mu (d)\lfloor \frac{M}{d}{\rfloor }^3$。

二元组方案数类似，为 ${\sum }_{d=1}^M\mu (d)\lfloor \frac{M}{d}{\rfloor }^2$。

一元组方案数显然是 $1$，只有 $1$ 满足 $\gcd =1$。

此时我们求出了不同的珠子个数为 $m$。

第二部分

显然求不同的环数需要套一手 $\text{burnside}$ 引理，问题变成求每个置换的不动点数量。对于一个长度为 $i$ 的置换，循环数显然是 $\gcd (i,n)$，每个循环内要求珠子一样，但是又要求相邻的珠子不能一样，即相邻的循环珠子不能一样，于是令 $f(n)$ 表示长度为 $n$ 的项链，相邻珠子不同，有多少种方案，则答案为：
$$\frac{1}{n}\sum _{i=1}^{n}f(\gcd (i,n))$$

注意，$f(n)$ 求的项链数不考虑旋转，即 $A$ 就算能旋转得到 $B$，也不认为是一样的方案，这样定义的原因比较显然，因为此时我们在求不动点数。

但是这样没法求解，设 $d=\gcd (i,n)$，考虑枚举 $d$，然后求有多少个 $i$ 满足 $\gcd (i,n)=d$，这样的 $i$ 可以表示为 $i=dk$，其中 $\gcd (\frac{n}{d},k)=1$，所以满足要求的 $k$ 有 $\phi (\frac{n}{d})$ 个，于是可以得到：
$$\sum _{i=1}^{n}f(\gcd (i,n))=\sum _{d|n}\phi (\frac{n}{d})f(d)$$

其中，$\phi$ 不能用筛法求，但是可以用 $\phi (n)={\prod }_{i=1}^k(p_i^{c_i}-p_i^{c_i-1})$ 这个式子来求。

剩下的问题是如何求 $f$，考虑往项链末尾放珠子，不难发现 $f$ 满足这个递推式：$f(n)=(m-2)f(n-1)+(m-1)f(n-2)$，第一项表示放一个 与相邻两颗珠子不同的珠子；第二项表示先放一个与第一颗珠子相同的珠子，再放一个与第一颗不同的珠子。

考虑用生成函数求通项公式，令 $F(x)={\sum }_{i=0}f(i)x^i$，将递推式写成封闭形式：
$$F-f(1)x-f(2)x^2=(m-2)(F-f(1)x)x+(m-1)F{x}^2$$

由于边界为 $f(1)=0,f(2)=m^2-m$，代入然后一顿转化可以得到：
$$F=\frac{(m^2-m)x}{1-(m-2)x+(1-m)x^2}$$

这种形式没办法展开，能展开的只有 $\frac{A}{1-ax}$ 这样的形式，于是我们想要他变成这样子的形式：
$$F=\frac{A}{1-ax}+\frac{B}{1-bx}=\frac{A(1-bx)+B(1-ax)}{1-(a+b)x+ab{x}^2}$$

观察分母，容易发现 $a,b$ 满足：
$$\{\begin{array}{l}a+b=m-2\\ ab=1-m\end{array}$$

解得 $a=m-1,b=-1$，代入到分子里，考虑求出 $A,B$：
$$\begin{array}{rl}A(1+x)+B(1-(m-1)x)& =(m^2-m)x\\ A+Ax+B+Bx-Bmx& =m^{2}x-mx\end{array}$$

因为 $A,B$ 中不含 $x$，而等式右边每一项都有 $x$，于是可以知道 $A+B=0$，即 $A,B$ 互为相反数，带上去消一下就可以得到：
$$\begin{array}{rl}-Bmx& =m^{2}x-mx\\ B& =1-m\end{array}$$

所以 $A=m-1$，于是 $F=\frac{m-1}{1-(m-1)x}-\frac{m-1}{1+x}$，展开就是 $F={\sum }_{i=0}((m-1{)}^i-(-1{)}^i(m-1))x^i$，所以 $f(i)=(m-1{)}^i-(-1{)}^i(m-1)$。

然后就做完了，但是要注意 $n$ 可能为模数的倍数，即 $p|n$，此时没法求 $n$ 的逆元。

考虑以 $p^2$ 作为模数，那么 $n$ 的逆元相当于 $\frac{n}{p}$ 在模 $p$ 意义下的逆元，这样就可以求出逆元了。但此时 $Ans$ 乘逆元求出来的是模 $p^2$ 意义下的解，由于 $Ans$ 是 $n$ 的倍数，即 $p|n|Ans$，所以 $Ans$ 可以直接除以 $p$ 得到模 $p$ 意义下的答案。

代码如下：

#include <cstdio>
#include <algorithm>
using namespace std;
#define maxn 10000010
#define ll long long

int T;
ll n,M,m,mod,Type;
const ll Mod1=1000000007,Mod2=1000000014000000049ll;
const ll inv61=166666668,inv62=833333345000000041ll;
int prime[maxn],pt=0,mu[maxn];
bool v[maxn];
void SieveInit(){
	mu[1]=1;
	for(int i=2;i<=maxn-10;i++){
		if(!v[i])prime[++pt]=i,mu[i]=-1;
		for(int j=1;j<=pt&&i*prime[j]<=maxn-10;j++){
			v[i*prime[j]]=true;
			if(i%prime[j]==0)break;
			mu[i*prime[j]]=-mu[i];
		}
		mu[i]+=mu[i-1];
	}
}
ll Mul(ll x,ll y){
	if(!Type)return x*y%mod;
	else return (x*y-(ll)(((long double)x*y+0.5)/((long double)mod))*mod+mod)%mod;
}
void solve_m(){
	m=2;
	for(int l=1,r;l<=M;l=r+1){
		r=M/(M/l);
		m=(m+Mul(Mul(Mul(M/l,M/l),M/l+3),(mu[r]-mu[l-1]+mod)%mod))%mod;
	}
	m=Mul(m,!Type?inv61:inv62);
}
ll ksm(ll x,ll y){ll re=1;for(;(y&1?re=Mul(re,x):0),y;y>>=1,x=Mul(x,x));return re;}
ll inv(ll x){x%=Mod1;ll y=Mod1-2,re=1;for(;(y&1?re=1ll*re*x%Mod1:0),y;y>>=1,x=1ll*x*x%Mod1);return re;}
ll yz[50],cnt[50],t=0,ans;
ll f(ll x){return (ksm(m-1,x)+(x&1?mod-m+1:m-1))%mod;}
void dfs(int x,ll now,ll phi){
	if(x>t){
		ans=(ans+Mul(phi%mod,f(n/now)))%mod;
		return;
	}
	for(int i=0;i<=cnt[x];i++){
		dfs(x+1,now,phi);
		now*=yz[x];
		phi*=i==0?yz[x]-1:yz[x];
	}
}
void solve(){
	ans=t=0;ll N=n;
	for(ll i=2;i*i<=N;i++)if(N%i==0){
		yz[++t]=i;cnt[t]=0;
		while(N%i==0)N/=i,cnt[t]++;
	}
	if(N>1)yz[++t]=N,cnt[t]=1;
	dfs(1,1,1);
}

int main()
{
	SieveInit();
	scanf("%d",&T);while(T--)
	{
		scanf("%lld %lld",&n,&M);
		if(n%Mod1)mod=Mod1,Type=0;
		else mod=Mod2,Type=1;
		solve_m();solve();
		if(!Type)ans=ans*inv(n%Mod1)%Mod1;
		else ans=ans/Mod1*inv(n/Mod1)%Mod1;
		printf("%lld\n",ans);
	}
}