正题

题目链接:https://www.luogu.com.cn/problem/P3307

---

题目大意

$n$个珠子的一个环形项链，每个珠子有三个$1\sim k$的整数。

• 两个珠子不同当且仅当它们不能通过翻转或者旋转得到
• 两个项链不同当且仅当它们不能通过旋转得到
• 珠子要求上面的数字互质
• 项链要求相邻珠子不同

求方案数。

---

解题思路

珠子的计数和项链计数是没什么关系的，所以两个分开求。

先考虑求珠子有多少种。相当与求有多少种三元组，因为三个数字的排列通过以上操作可以得到任何其他排列。
首先答案包括了三个不同的，有两个相同的，全部都相同的。

也就是设$S_i$表示$i$元组的数目的话答案就是$S_3+3S_2+2S_1$（容斥一下可以得到的）

可以先不考虑顺序问题，最后除上一个六就好了。

要求互质所以可以直接上莫反就变为求${\sum }_{i=1}^k\mu (i)F(i)$。对于每个$F(i)$来说就是有$\lfloor \frac{k}{i}\rfloor$个数字不要求互质的答案。
对于三元组来说$F(i)=\lfloor \frac{k}{i}{\rfloor }^3$，二元组就是$F(i)=\lfloor \frac{k}{i}{\rfloor }^2$可以和三元组一起解决，一元组就是$1$。
这样整除分块就可以$\sqrt{k}$的处理单次询问了。

现在考虑项链的方案，根据$\text{Burnside}$引理的话就是$\frac{1}{n}{\sum }_{i=1}^{n}f(gcd(i,n))$，其中$f(i)$表示大小为$i$的不动点数量（也就是转$i$次的话置换中环的数量）。化一下就是$\frac{1}{n}{\sum }_{d|n}f(d)\phi (\frac{n}{d})$，需要考虑如何求$f$。

因为是不动点，所以每一个环中的颜色都得相同，然后因为每个环都是相互隔开的，也就是可以视为有$gcd(n,i)$个珠子头尾相连要涂色，不能有相邻颜色相同。
这个$f$是有递推式$f(i)=f(i-1)\ast (m-2)+f(i-2)\ast (m-1)$（据说法是第一个是直接在末端插入一个颜色的珠子，第二个是先在末端插入一个和首部相同颜色的珠子再插入一个不同颜色的将它们隔开，这个是第一个无法处理的情况）。

这个可以用矩阵乘法快速递推，但是好像有更简单的方法可以一试，考虑将式子化为$f(x)-af(x-1)=bf(x-1)+abf(x-2)$的形式，这样可以将答案化为一个带幂的函数。

那就有方程$\{\begin{array}{c}b-a=m-2\\ ab=m-1\end{array}$，然后可以解出$a=1$或者$a=1-m$。带个$1$进去做。也就是有$f(x)-f(x-1)=(m-1)(f(x-1)+f(x-2))$也就是乘上$m-1$就往后退一步，所以有$f(x)-f(x-1)=(m-1{)}^{x-2}(f(1)+f(2))$直接带$f(1)=0,f(2)=m(m-1)$进去就有$f(x)=m(m-1{)}^{x-1}-f(x-1)$
设一个$A$满足
$$f(x)+(m-1)A=-(f(x-1)+A)$$
$$f(x)+mA=-f(x-1)$$
又因为有$f(x)=m(m-1{)}^{x-1}-f(x-1)$
$$m(m-1{)}^{x-1}+mA=0\Rightarrow A=(m-1{)}^{x-1}$$
所以把$A$带回最初式子，把$f(x-1)$换成$f(x)$就有式子
$$f(x)+(m-1{)}^x=-(f(x)-(m-1{)}^{x-1})$$
$$\Rightarrow f(x)=(-1{)}^x(m-1)+(m-1{)}^x$$
（这个好像是叫特征方程的玩意）

现在就可以用快速幂算$f(x)$了，$\phi$比较大，不能直接$\sqrt{n}$求单个，可以先把$n$质因数分解了然后$dfs$下去就好了。然后还有因为$n$很大，有可能是模数$P$的倍数，这个时候我们求答案模$P^2$，然后除以一个$P$在和$\frac{n}{P}$搞逆元就好了。

时间复杂度$O(k+T(\sqrt{n}+\sqrt{m}))$

---

code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const ll N=1e7+10,XJQ=1e9+7,inv6=833333345000000041ll;
ll T,n,m,k,ans,P,cnt,num;
ll pri[N/10],mu[N],p[110],c[110];
bool v[N];
void init(){mu[1]=1;for(ll i=2;i<N;i++){if(!v[i])pri[++num]=i,mu[i]=-1;for(ll j=1;j<=num&&i*pri[j]<N;j++){v[i*pri[j]]=1;dans}for(ll i=1;i<N;i++)mu[i]+=mu[i-1];return;
}
ll ksc(ll a,ll b,ll p){a%=p;b%=p;ll tmp=(long double)a*b/p;long double ans=a*b-tmp*p;if(ans>=p)ans-=p;else if(ans<0)ans+=p;return ans;
}
ll power(ll x,ll b,ll p){ll ans=1;while(b){if(b&1)ans=ksc(ans,x,p);x=ksc(x,x,p);b>>=1;}return ans;
}
void Get_M(ll n){m=2;for(ll l=1,r;l<=n;l=r+1){r=n/(n/l);ll k=n/l;(m+=ksc(ksc(ksc(k,k,P),k+3,P),mu[r]-mu[l-1],P))%=P;}m=ksc(m,inv6,P);return;
}
ll f(ll n)
{return power(m-1,n,P)+((n&1)?(P-m+1):(m-1))%P;}
void dfs(ll dep,ll x,ll val){if(dep>cnt){(ans+=ksc(f(n/x),val,P))%=P;return;}dfs(dep+1,x,val);val*=p[dep]-1;x*=p[dep];dfs(dep+1,x,val);for(ll i=2;i<=c[dep];i++)val=val*p[dep],x*=p[dep],dfs(dep+1,x,val);return;
}
signed main()
{scanf("%lld",&T);init();while(T--){scanf("%lld%lld",&n,&k);if(n%XJQ!=0)P=XJQ;else P=XJQ*XJQ;Get_M(k);// printf("%lld",m);ll x=n;cnt=ans=0;for(ll i=1;i<=num&&pri[i]*pri[i]<=x;i++){if(x%pri[i]==0){p[++cnt]=pri[i];c[cnt]=0;while(x%pri[i]==0)c[cnt]++,x/=pri[i];}}if(x>1)p[++cnt]=x,c[cnt]=1;dfs(1,1,1);if(n%XJQ==0)printf("%lld\n",ans/XJQ*power(n/XJQ,XJQ-2,XJQ)%XJQ);elseprintf("%lld\n",ans*power(n,XJQ-2,XJQ)%XJQ);}return 0;
}