【SDOI2013】项链【莫比乌斯反演】【Polya定理】【递推式求通项】【数论】

题意：$T$ 组数据，每组给定 $n,a$，求满足下列条件的项链数量：

• 有 $n$ 个珠子。
• 每个珠子上有三个 $[1,a]\cap \mathbb{Z}$ 的数，且三个数 $\gcd$ 为 $1$。
• 相邻两个珠子不同。

珠子旋转、翻转同构，项链旋转同构。对 $10^9+7$ 取模。

$n\le 10^{14},a\le 10^7$

有毒……

显然是道缝合怪题，先考虑求有多少个不同的珠子。

设 $f(k)$ 表示 $k$ 个 $[1,a]$ 的数 $\gcd$ 为 $1$ 的方案数，用 Polya 定理数一下可知

$$ans=\frac{f(3)+3f(2)+2f(1)}{6}$$

然后就是个简单的反演

$$f(k)=\sum _{i=1}^a\mu (i){\lfloor \frac{a}{i}\rfloor }^k$$

整除分块算就可以了（其实直接暴力也可以？）。

然后第二步，数项链。设第一步得到的珠子种数为 $k$，众所周知，$F(n)$ 为 $n$ 个珠子不循环同构的方案，答案就是 $F$ 和 $\phi$ 的狄利克雷卷积除以 $n$。

考虑 $F$ 怎么算。

考虑一个 $n$ 个珠子的合法方案，我们删掉编号为 $1$ 的点。如果两边的点不同，就对应了一种 $F(n-1)$ 的方案，乘上 $1$ 号珠子的方案 $(k-2)$。如果相同，任意删掉一个后就对应一种 $F(n-2)$ 的方案，乘上 $1$ 号珠子的方案 $(k-1)$。所以

$$F(n)=(k-2)F(n-1)+(k-1)F(n-2)$$

注意边界是 $F(1)=0,F(2)=k(k-1)$,倒退出来 $F(0)=m$。

列出生成函数：

$$F=(k-2)xF+(k-1)x^{2}F+k-k(k-2)x$$

这里减 $k(k-2)x$ 是为了平衡 $(k-2)xF$ 对一次项的影响。

即

$$F=\frac{-m(m-2)x+m}{(1+x)(1-(m-1)x)}$$

直接待定系数拆掉

$$F=\frac{a}{1+x}+\frac{b}{1-(k-1)x}$$

列出方程：

$$a-a(k-1)x+b+bx=-k(k-2)x+k$$

解得 $a=k-1,b=1$

所以

$$F=\frac{k-1}{1+x}+\frac{1}{1-(k-1)x}$$

$$F(n)=(k-1{)}^n+(-1{)}^n(k-1)$$

然后就可以算了。

这里有点卡，就不用 Polya 模板那个垃圾做法，而是把因子预处理出来 dfs，可以做到 $\mathrm{O}(d(n))$ 的优秀复杂度。

然后因为 $n$ 很大，最后除 $n$ 的时候又会有喜闻乐见的除 $0$ 问题。所以要对 $(10^9+7{)}^2$ 取模，最后再判一下。

总复杂度 $\mathrm{O}(n+T(\sqrt{a}+d(n)\log ))$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#define int long long
using namespace std;
const int mod=1e9+7,MOD=mod*mod;
inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}
inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}
inline int mul(int a,int b){return (a*b-(int)((long double)a/MOD*b+0.5)*MOD+MOD)%MOD;}
inline int qpow(int a,int p)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=mul(ans,a);
		a=mul(a,a),p>>=1;
	}
	return ans;
}
int INV6=833333345000000041ll;
int n,a;
const int N=1e7,MAXN=N+5;
bool np[MAXN];
signed pl[MAXN],cnt,mu[MAXN];
void init()
{
	np[1]=1,mu[1]=1;
	for (int i=2;i<=N;i++)
	{
		if (!np[i]) mu[pl[++cnt]=i]=-1;
		for (int j=1,x;(x=i*pl[j])<=N;j++)
		{
			np[x]=1;
			if (i%pl[j]==0) break;
			mu[x]=-mu[i];
		}
	}
	for (int i=2;i<=N;i++) mu[i]+=mu[i-1];
}
int calc()
{
	int ans=0;
	for (int l=1,r;l<=a;l=r+1)
	{
		r=a/(a/l);
		int t=mu[r]-mu[l-1];
		(t<0)&&(t+=MOD);
		ans=add(ans,mul(t,add(mul(a/l,mul(a/l,a/l)),3*mul(a/l,a/l)%MOD)));
	}
	return ans;
}
int k,ans;
int lis[20005],siz[20005],tot;
void dfs(int pos,int phi,int res)
{
	if (pos>tot)
	{
		int t=qpow(k,res);
		if (res&1) t=dec(t,k);
		else t=add(t,k);
		ans=add(ans,mul(phi,t));
		return;
	}
	int t=qpow(lis[pos],siz[pos]);
	dfs(pos+1,phi,res*t);
	t/=lis[pos];
	dfs(pos+1,phi*=(lis[pos]-1),res*t);
	t/=lis[pos];
	if (!t) return;
	for (;t;dfs(pos+1,phi*=lis[pos],res*t),t/=lis[pos]);
}
signed main()
{
	init();
	int T;
	cin>>T;
	while (T--)
	{
		cin>>n>>a;
		k=dec(mul(calc()+2,INV6),1);
		tot=ans=0;
		int t=n;
		for (int i=1;i<=cnt;i++)
			if (t%pl[i]==0)
			{
				lis[++tot]=pl[i],siz[tot]=0;
				while (t%pl[i]==0) t/=pl[i],++siz[tot];
			}
		if (t>1) lis[++tot]=t,siz[tot]=1;
		dfs(1,1,1);
		if (n%mod==0) n/=mod,ans/=mod;
		ans=mul(ans,qpow(n,mod-2));
		cout<<ans%mod<<'\n';
	}
	return 0;
}