其他分享
首页 > 其他分享> > SDOI2013 项链

SDOI2013 项链

作者:互联网

首先我们可以求出来满足条件的珠子的总数 \(m\)。有以下结论:

我是硬推出来这个式子的2333

具体来说,设 \(f(n)\) 表示大小为 \(n\) 的环的答案,考虑容斥:

  • 先不管第 \(1\) 个元素和第 \(n\) 个元素之间的限制,那么可以发现答案就是 \(m\times (m-1)^{n-1}\)。这是因为第 \(1\) 个元素有 \(m\) 种选择,后面的元素都只有 \(m-1\) 种选择。
  • 再去掉第一个元素和第 \(n\) 个元素相同时的答案,此时可以把这两个元素看成一个,那么方案数就是 \(f(n-1)\)。
  • 但需要特判 \(n\le 2\) 的情况:\(f(1)=m,f(2)=m\times(m-1)\)。

我们有 \(f(n)=m\times (m-1)^{n-1}-f(n-1)\),其中 \(n\ge 3\)。展开得到

\[f(n)=m\times \sum_{i=1}^{n-1}(-1)^{n+1-i}(m-1)^i=m\times\left(\sum_{i=0}^{n-1}(-1)^{n+1-i}(m-1)^i\right)-m \]

把后面个求和里面的 \((-1)^{n+1}\) 提出来得到

\[f(n)=m\times (-1)^{n+1}\times \sum_{i=0}^{n-1}(-1)^{-i}(m-1)^i=m\times (-1)^{n+1}\times \sum_{i=0}^{n-1}(1-m)^i-(-1)^{n+1}m \]

由等比数列求和公式有

\[\begin{aligned} f(n)&=m\times (-1)^{n+1}\times\dfrac{1-(1-m)^{n}}{1-(1-m)}=(-1)^{n+1}\times (1-(1-m)^n)-(-1)^{n+1}m\\ &=(-1)^{n+1}+(m-1)^n-(-1)^{n+1}m=(-1)^{n}(m-1)+(m-1)^n \end{aligned} \]

这就证明了结论。

考虑怎么算出 \(m\)。再用一次 \(\text{Burnside}\) 引理,我们设

\[S_1=1\\ S_2=\sum_{i=1}^a\sum_{j=1}^a[\gcd(i,j)=1]\\ S_3=\sum_{i=1}^a\sum_{j=1}^a\sum_{k=1}^a[\gcd(i,j,k)=1] \]

这些都可以 \(O(\sqrt{a})\) 求解。那么 \(m=\frac{1}{6}(2S_1+3S_2+S_3)\)。

本题中置换群由 \(n\) 种旋转构成。

如果转了 \(k\) 格,那么会形成 \(\gcd(k,n)\) 个置换环,答案就是 \((m-1)^{\gcd(k,n)}+(-1)^{\gcd(k,n)+1}\)。显然可以 \(O(\sqrt{n}\log n)\) 求出来。

看上去 \(n=10^{14}\) 非常不可过,但是你发现这个算法的复杂度实际上是 \(O(d(n)\log n)\),我们都知道 \(d(n)\sim O(\sqrt[3]{n})\),于是这个算法也就可以通过本题了。

#include<bits/stdc++.h>

#define int long long

using namespace std;

int read(){
	int x=0,f=1;char c=getchar();
	for(;(!(c>='0'&&c<='9'));c=getchar())if(c=='-')f=-1;
	for(;(c>='0'&&c<='9');c=getchar())x=x*10+c-'0';
	return x*f;
}

const int N=1e7;
int pri[N+5],cnt=0,mu[N+5];
bool vis[N+5];

void init(){
	for(int i=2;i<=N;i++){
		if(!vis[i])pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pri[j]<=N;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
			mu[i*pri[j]]=-mu[i];
		}
	}
	mu[1]=1;for(int i=1;i<=N;i++)mu[i]+=mu[i-1];
}

const int mod=1e9+7;

int mul(int x,int y,int p){
	__int128 xx=x,yy=y,pp=p;
	return (int)(xx*yy%pp);
}

const int MOD=mod*mod;

int ksm(int x,int y,int p=MOD){
	int ans=1;
	for(int i=y;i;i>>=1,x=mul(x,x,p))if(i&1)ans=mul(ans,x,p);
	return ans;
}

int inv(int x,int p=MOD){
	return ksm(x,p-2,p)%p;
}

const int iv=833333345000000041;

int calc(int a){
	int res=0;
	for(int l=1,r=0;l<=a;l=r+1){
		r=a/(a/l);int w=a/l;
//		cout<<l<<" "<<r<<endl;
		res=(res+mul((mul(w,mul(w,w,MOD),MOD)+mul(3,mul(w,w,MOD),MOD)),(mu[r]-mu[l-1]+MOD)%MOD,MOD))%MOD;
//		cout<<res<<endl;
	}
	res+=2;
//	cout<<"res = "<<res<<endl;
	return mul(res,iv,MOD);
}

int p[30],c[30],r[30],k,n,m,ans=0;
void dfs(int now){
	if(now==k+1){
		int phi=1,d=1;
		for(int i=1;i<=k;i++){
			if(r[i]>0)phi=phi*ksm(p[i],r[i]-1)*(p[i]-1);
			d=d*ksm(p[i],c[i]-r[i]);
		}
//		cout<<d<<" | "<<phi<<endl;
		ans=(ans+(mul((ksm(m-1,d)%MOD+mul(d%2==0?1:-1,m-1,MOD)+MOD)%MOD,phi%MOD,MOD)))%MOD;
		return ;
	}
	for(int i=0;i<=c[now];i++)r[now]=i,dfs(now+1),r[now]=0;
}

void solve(){
	k=0;ans=0;
	n=read();int a=read();
	m=calc(a);int div=n;
//	cout<<m<<endl;
	for(int i=2;i*i<=n;i++){
		if(n%i==0){
			p[++k]=i,c[k]=0;
			while(n%i==0)n/=i,c[k]++;
		}
	}
	if(n>1)p[++k]=n,c[k]=1;
//	for(int i=1;i<=k;i++)cout<<p[i]<<" ^ "<<c[i]<<endl;
	dfs(1);
	if(div%mod!=0)cout<<ans%mod*inv(div%mod,mod)%mod<<endl;
	else cout<<(ans/mod)*inv(div/mod,mod)%mod<<endl;
}

signed main(void){

#ifndef ONLINE_JUDGE
	freopen("in.in","r",stdin);
#endif

	init();
	int tt=read();
	while(tt--)solve();

	return 0;
}

标签:gcd,int,sum,元素,times,SDOI2013,项链,ksm
来源: https://www.cnblogs.com/YunQianQwQ/p/16594455.html