牛客练习赛84F-牛客推荐系统开发之下班【莫比乌斯反演,杜教筛】

1|0正题

题目链接:https://ac.nowcoder.com/acm/contest/11174/F

1|1题目大意

给出 $n,k$ 求

\sum_{i_{1} = 1}^{n} \sum_{i_{2} = 1}^{n} . . . \sum_{i_{k} = 1}^{n} g c d (f_{i_{1}}, f_{i_{2}}, . . ., f_{i_{k}})

$\sum_{i_1=1}^n\sum_{i_2=1}^n...\sum_{i_k=1}^ngcd(f_{i_1},f_{i_2},...,f_{i_{k}})$

对 $10^9+9$ 取模

其中 $f_i$ 表示斐波那契数列的第 $i$ 项

$1\leq n,k\leq 10^8$

1|2解题思路

看上去就很莫比乌斯反演，首先把 $f$ 提出来然后直接上莫反就是

\sum_{i = 1}^{n} f_{i} \sum_{j = 1}^{⌊ \frac{n}{i} ⌋} μ (j) ⌊ \frac{n}{i j} ⌋^{k}

$\sum_{i=1}^nf_{i}\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\mu(j)\lfloor\frac{n}{ij}\rfloor^k$

这个式子其实就可以直接做了，外面 $f_i$ 是一个整除分块，然后里面的式子也是一个整除分块。需要的算法是一个矩阵乘法求斐波那契的前缀和（当然因为 $\sqrt 5$ 在 $10^9+9$ 下有二次剩余可以直接用特制方程求）还有一个杜教筛求 $\mu$ 的前缀和。

因为整除分块套整除分块的复杂度是 $O(n^{\frac{3}{4}})$ ，然后最里面套个杜教筛复杂度不会大很多所以能过。

然后有一个大佬的做法是把上面那个式子化一下

\sum_{i = 1}^{n} ⌊ \frac{n}{i} ⌋^{k} \sum_{d | i} f_{d} μ (\frac{i}{d})

$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor^k\sum_{d|i}f_{d}\mu(\frac id)$

（相当于外面枚举 $i\times j$ 然后枚举它的约数）

然后因为 $\mu*I=\epsilon$ ，所以 $f*\mu*I=f*\epsilon=f$ 。所以如果在能快速求 $f$ 的前缀和的情况下求 $f*\mu$ 的前缀和可以用杜教筛来做。

然后上面那个式子整除分块一下就好了

代码的写法是第一种

1|3code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#define ll long long
using namespace std;
const ll N=5e6+10,S=3,P=1e9+9;
struct Matrix{
	ll a[S][S];
}c,f,am;
ll n,k,ans,cnt,mu[N],pri[N/10];
bool v[N];map<ll ,ll >mp;
Matrix operator*(Matrix &a,Matrix &b){
	memset(c.a,0,sizeof(c.a));
	for(ll i=0;i<S;i++)
		for(ll j=0;j<S;j++)
			for(ll k=0;k<S;k++)
				(c.a[i][j]+=a.a[i][k]*b.a[k][j]%P)%=P;
	return c;
}
ll power(ll x,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=ans*x%P;
		x=x*x%P;b>>=1;
	}
	return ans;
}
ll Fbi(ll b){
	memset(f.a,0,sizeof(f.a));
	am.a[0][0]=am.a[0][2]=0;
	am.a[0][1]=1;f.a[2][2]=1;
	f.a[0][1]=1;f.a[1][1]=1;
	f.a[1][0]=1;f.a[1][2]=1;
	while(b){
		if(b&1)am=am*f;
		f=f*f;b>>=1; 
	}
	return am.a[0][2];
}
void Prime(){
	mu[1]=1;
	for(ll i=2;i<N;i++){
		if(!v[i])pri[++cnt]=i,mu[i]=-1;
		for(ll j=1;j<=cnt&&i*pri[j]<N;j++){
			v[i*pri[j]]=1;
			if(i%pri[j]==0)break;
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(ll i=1;i<N;i++)
		mu[i]+=mu[i-1];
	return;
}
ll GetMu(ll n){
	if(n<N)return mu[n];
	if(mp[n])return mp[n];
	ll ans=1;
	for(ll l=2,r;l<=n;l=r+1)
		r=n/(n/l),ans=ans-(r-l+1)*GetMu(n/l);
	return mp[n]=ans;
}
ll g(ll n){
	ll ans=0;
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ll tmp=GetMu(r)-GetMu(l-1);
		ans=(ans+tmp*power(n/l,k)%P)%P;
	}
	return ans;
}
signed main()
{
	Prime();
	scanf("%lld%lld",&n,&k);
	for(ll l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ll tmp=Fbi(r)-Fbi(l-1);
		ans=(ans+tmp*g(n/l)%P)%P;
	}
	printf("%lld\n",(ans+P)%P);
	return 0;
}

__EOF__

本文作者：QuantAsk
本文链接：https://www.cnblogs.com/QuantAsk/p/14879690.html
关于博主：退役OIer，GD划水选手
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！