BSOJ5467 [CSPX2017#3]整数 莫比乌斯反演+杜教筛

题意简述

给你两个整数(n)(k),让你求出这个式子

[sum_{a_1=1}^n sum_{a_2=a_1}^n sum_{a_3=a_2}^n cdots sum_{a_k=a_{k-1}}^n left[ gcd {(a_1,a_2,a_3cdots,a_k)} = 1 ight] ]

做法

对于(gcd)进行莫比乌斯反演

[Ans = sum_p mu(p) sum_{a_1=1}^n sum_{a_2=a_1}^{frac{n}{p}} sum_{a_3=a_2}^{frac{n}{p}} cdots sum_{a_k=a_{k-1}}^{frac{n}{p}} 1 ]

(S(n)= sum_{a_1=1}^n sum_{a_2=a_1}^n sum_{a_3=a_2}^n cdots sum_{a_k=a_{k-1}}^n 1),可以发现

[S(n)=sum [1leq a_1 leq a_2 leq cdots leq a_kleq n]=inom {n+k-1}{k} ]

所以可得

[Ans=sum_p mu (p) Sleft(frac{n}{p} ight) = sum_p mu (p) inom {frac{n}{p}+k-1}{k} ]

杜教筛预处理(mu),整除分块计算(Sleft(frac{n}{p} ight))即可。

代码实现

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define re register int
#define db double
#define in inline 
#define rd regitser double
#define ak *
const int N=1e7+5,p=1e9+7,K=1e3+5;
char qwq;
int pri[N>>3];
ll mu[N],inv[N],fac[N],fav[N];
bool vis[N];
unordered_map<ll,ll>smu;
inline int read()
{
	re cz=0,ioi=1;qwq=getchar();
	while(!isdigit(qwq)) ioi=qwq=='-'?~ioi+1:1,qwq=getchar();
	while(isdigit(qwq)) cz=(cz<<3)+(cz<<1)+(qwq^48),qwq=getchar();
	return cz ak ioi;
}
in void get()
{
	mu[1]=1;
	for(re i=2;i<=1e7;i++)
	{
		if(!vis[i]) pri[++pri[0]]=i,mu[i]=p-1;
		for(re j=1;j<=pri[0]&&i*pri[j]<=1e7;j++)
		{
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) {mu[i*pri[j]]=0;break;}
			else mu[i*pri[j]]=p-mu[i];
		}
	}
	for(re i=2;i<=1e7;i++) mu[i]=(mu[i]+mu[i-1])%p;
}
in ll qpow(ll x,ll y,ll z=1)
{
	for(;y;y>>=1,x=x*x%p) if(y&1) z=z*x%p;
	return z;
}
in ll c(ll n,ll m)
{
	if(n<=1e7) return fac[n]*inv[n-m]%p*inv[m]%p;
	ll ans=1ll;
	for(ll i=0;i<m;i++) ans=ans*(n-i)%p;
	return ans*inv[m]%p;
}
in ll get_mu(ll n)
{
	if(n<=1e7) return mu[n];
	if(smu[n]) return smu[n];
	ll res=1;
	for(ll l=2,r;l<=n;l=r+1) 
	r=n/(n/l),res=(res-(r-l+1)*get_mu(n/l)%p+p)%p;
	return smu[n]=res;
}
int main()
{
	get();inv[0]=inv[1]=fav[0]=fac[0]=fac[1]=1;
	for(re i=2;i<=1e7;i++) inv[i]=(ll)(p-p/i)*inv[p%i]%p,fac[i]=fac[i-1]*i%p;
	for(re i=1;i<=1e7;i++) inv[i]=inv[i-1]*inv[i]%p;
	re opt=read();
	while(opt--)
	{
		ll n=read(),k=read();
		ll ans=0;
		for(ll l=1,r;l<=n;l=r+1)
		{
			r=n/(n/l); 
			ans=(ans+(ll)(get_mu(r)-get_mu(l-1)+p)%p*c(n/l+k-1,k)%p)%p;
		}
		cout<<ans<<endl;
	}
	return 0;
}
原文地址:https://www.cnblogs.com/disangan233/p/11181301.html