洛谷P3307 [SDOI2013] 項鍊

設 \(m\) 爲有多少種不一樣的珠子，先用 \(Burnside\) 引理求 \(m\)，發現對於珠子的置換有旋轉和翻轉，得其置換羣有 \(|G|=6\)，對於有 \(3\) 個不一樣數字的珠子，得其穩定化子 \(|G_x|=1\)，對於有 \(2\) 個不一樣數字的珠子，得其穩定化子 \(|G_x|=2\)，對於數字都相同的珠子，得其穩定化子 \(|G_x|=6\)。設這三類珠子的個數爲 \(S_3,S_2,S_1\)，得：

\[\large m=\frac{1}{|G|}\sum_{f\in G}c(f)=\frac{1}{6}\left( S_3+2S_2+6S_1 \right) \]

注意到：

\[\large\begin{aligned} \sum_{i=1}^a\sum_{j=1}^a\sum_{k=1}^a[\gcd(i,j,k)=1]&=S_3+S_2+S_1\\ \sum_{i=1}^a\sum_{j=1}^a[\gcd(i,j)=1]&=\frac{1}{3}S_2+S_1\\ S_1&=1 \end{aligned} \]

所以得：

\[\large m=\frac{1}{6}\left( \sum_{i=1}^a\left\lfloor\frac{a}{i}\right\rfloor^3\mu(i)+3\sum_{i=1}^a\left\lfloor\frac{a}{i}\right\rfloor^2\mu(i)+2 \right) \]

考慮求環的方案數，設 \(f_n\) 表示 \(n\) 個點的環染色的方案數，即長爲 \(n\) 的序列染色且首尾不一樣的方案數，環的置換隻有旋轉，由 \(Burnside\) 引理得方案數爲：

\[\large \frac{1}{n}\sum_{i=1}^nf_{\gcd(i,n)}=\frac{1}{n}\sum_{d\mid n}f_d\sum_{i=1}^n[\gcd(i,n)=d]=\frac{1}{n}\sum_{d\mid n}f_d\varphi\left(\frac{n}{d}\right) \]

不可貴到 \(f_n\) 的遞推式：

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

考慮第 \(1\) 個位置和第 \(n-1\) 個位置顏色是否相同便可。用特徵根求通項公式得：

\[\large\begin{aligned} x^2&-(m-2)x-(m-1)=0\\ x_1&=m-1,x_2=-1\\ f_n&=c_1x_1^n+c_2x_2^n\\ f_1&=c_1(m-1)-c_2=0\\ f_2&=c_1(m-1)^2+c_2=m(m-1)\\ f_n&=(m-1)^n+(m-1)(-1)^n\\ \end{aligned} \]

而後就能快速冪計算 \(f_n\) 了。

由於 \(n\) 達到 \(10^{14}\)，因此求 \(\varphi(n)\) 的時候應該先質因數分解 \(n\)，而後 \(dfs\) 每個因數，搜索的過程當中求出 \(\varphi(n)\)。

一樣是由於 \(n\) 很大，可能會出現 \(p \mid n\) 的狀況，這樣的話，最後除 \(n\) 時可能不存在逆元。因而統計答案的時候對 \(p^2\) 取模，由於除 \(n\) 前答案必定是 \(n\) 的倍數，那麼其也必定是 \(p\) 的倍數，因而就直接除 \(p\)，而後再除 \(\frac{n}{p}\)，由於 \(n< p^2\)，因此 \(\frac{n}{p}\) 在模 \(p\) 意義下是存在逆元的。

#include<bits/stdc++.h>
#define maxn 10000010
#define all 10000000
#define p 1000000007
#define inv6 (ll)833333345000000041
#define P ((ll)1000000007*1000000007)
using namespace std;
typedef long long ll;
template<typename T> inline void read(T &x)
{
	x=0;char c=getchar();bool flag=false;
	while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
	while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	if(flag)x=-x;
}
ll T,n,m,tot,cnt,ans;
ll pri[maxn],mu[maxn],v[maxn],k[maxn];
bool tag[maxn];
void init()
{
	mu[1]=1;
	for(int i=2;i<=all;++i)
	{
		if(!tag[i]) pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot;++j)
		{
			int k=i*pri[j];
			if(k>all) break;
			tag[k]=true;
			if(i%pri[j]) mu[k]=-mu[i];
			else break;
		}
	}
	for(int i=1;i<=all;++i) mu[i]+=mu[i-1];
}
ll mul(ll x,ll y)
{
	ll c=(long double)x*y/P+0.5;
	c=x*y-c*P;
	return c<0?c+P:c;
}
ll qp(ll x,ll y)
{
	ll v=1;
	while(y)
	{
		if(y&1) v=mul(v,x);
		x=mul(x,x),y>>=1;
	}
	return v;
}
ll inv(ll x)
{
	ll v=1,y=p-2;
	while(y)
	{
		if(y&1) v=v*x%p;
		x=x*x%p,y>>=1;
	}
	return v;
}
void get(ll n)
{
	cnt=0;
	for(int i=1;i<=tot&&(ll)pri[i]*pri[i]<=n;++i)
	{
		if(n%pri[i]) continue;
		v[++cnt]=pri[i],k[cnt]=0;
		while(n%pri[i]==0) n/=pri[i],k[cnt]++;
	}
	if(n!=1) v[++cnt]=n,k[cnt]=1;
}
void work(ll n)
{
	ll v1=0,v2=0;
	for(ll l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		v1=(v1+mul(mul(mul(n/l,n/l),n/l),(mu[r]-mu[l-1]+P)%P))%P;
		v2=(v2+mul(mul(n/l,n/l),(mu[r]-mu[l-1]+P)%P))%P;
	}
	m=mul((v1+mul(v2,3)+2)%P,inv6);
}
ll calc(ll n)
{
	return (qp(m-1+P,n)+((n&1)?1-m+P:m-1+P)%P)%P;
}
void dfs(int x,ll d,ll phi)
{
    if(x==cnt+1)
    {
        ans=(ans+mul(phi,calc(n/d)))%P;
        return;
    }
    dfs(x+1,d,phi),d*=v[x],phi*=v[x]-1,dfs(x+1,d,phi);
    for(int i=2;i<=k[x];++i) d*=v[x],phi*=v[x],dfs(x+1,d,phi);
}
int main()
{
	init(),read(T);
	while(T--)
	{
		read(n),read(m),get(n),work(m),ans=0,dfs(1,1,1);
		if(n%p) ans=ans%p*inv(n%p)%p;
		else ans=ans/p*inv(n/p)%p;
		printf("%lld\n",ans);
	}
	return 0;
}