【牛客练习赛84 F】牛客推荐系统开发之下班

题目

题目链接：https://ac.nowcoder.com/acm/contest/11174/F
给定 \(n,k\)，求

\[\sum^{n}_{i_1=1}\sum^{n}_{i_2=1}\cdots \sum^{n}_{i_k=1}\gcd(\text{fib}_{i_1},\text{fib}_{i_2},\cdots,\text{fib}_{i_k}) \]

对 \(10^9+9\) 取模后的值。

\(n,k\leq 10^8\)。

思路

直接上莫比乌斯反演。

\[\sum^{n}_{i_1=1}\sum^{n}_{i_2=1}\cdots \sum^{n}_{i_k=1}\gcd(\text{fib}_{i_1},\text{fib}_{i_2},\cdots,\text{fib}_{i_k}) \]

\[=\sum^{n}_{i_1=1}\sum^{n}_{i_2=1}\cdots \sum^{n}_{i_k=1}\text{fib}_{\gcd(i_1,i_2,\cdots,i_k)} \]

\[=\sum^{n}_{i=1}\text{fib}_i\left ( \sum_{i|j}\mu\left (\frac{j}{i}\right )\lfloor \frac{n}{j} \rfloor^k \right) \]

\[=\sum^{n}_{i=1}\text{fib}_i\left ( \sum^{\lfloor\frac{n}{i}\rfloor}_{j=1}\mu(j)\lfloor\frac{n}{ij}\rfloor^k \right) \]

令 \(g(i)=\sum^{i}_{j=1}\mu(j)\lfloor\frac{n}{j}\rfloor^k\)，则原式

\[=\sum^{n}_{i=1}\text{fib}_i\times g\left( \frac{n}{i} \right) \]

考虑到 \(g\) 只有 \(\sqrt{n}\) 种取值，\(g(i)\) 可以通过整除分块 + 杜教筛求出，考虑直接暴力求出 \(g\)。事实上暴力求出所有 \(g\) 的复杂度是 \(O(n^{\frac{3}{4}})\)，但是我并不会证明。
然后原式进行整除分块后我们需要快速求斐波那契数列区间和。有一个经典的结论：

\[\sum^{n}_{i=1}\text{fib}_i=\text{fib}_{i+2}-1 \]

那么问题又转化为求 \(O(\sqrt{n})\) 个斐波那契数的值。可以通过矩阵乘法求出。
时间复杂度 \(O(n^{\frac{3}{4}})\)。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=10000010,MOD=1e9+9;
int n,m,k,mu[N],prm[N];
bool v[N];
ll ans;
map<int,int> smu;

void findprm(int n)
{
	mu[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++m]=i,mu[i]=-1;
		for (int j=1;j<=m;j++)
		{
			if (i>n/prm[j]) break;
			v[i*prm[j]]=1; mu[i*prm[j]]=-mu[i];
			if (!(i%prm[j])) { mu[i*prm[j]]=0; break; }
		}
	}
	for (int i=2;i<=n;i++) mu[i]+=mu[i-1];
}

ll fpow(ll x,ll k)
{
	ll res=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) res=res*x%MOD;
	return res;
}

ll getmu(int n)
{
	if (n<N) return mu[n];
	if (smu[n]) return smu[n];
	ll res=1;
	for (int i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		res=(res-getmu(n/i)*(j-i+1))%MOD;
	}
	return smu[n]=res;
}

struct Matrix
{
	ll a[3][3];
	Matrix() { memset(a,0,sizeof(a)); }
	
	friend Matrix operator *(Matrix a,Matrix b)
	{
		Matrix c;
		for (int i=1;i<=2;i++)
			for (int j=1;j<=2;j++)
				for (int k=1;k<=2;k++)
					c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%MOD;
		return c;
	}
}f,a;

ll fib(int n)
{
	if (!n) return 0;
	if (n<=2) return 1;
	f.a[1][1]=f.a[1][2]=1;
	a.a[1][1]=a.a[1][2]=a.a[2][1]=1; a.a[2][2]=0;
	n-=2;
	for (;n;n>>=1,a=a*a)
		if (n&1) f=f*a;
	return f.a[1][1];
}

ll calc(int n)
{
	ll res=0;
	for (int i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		res=(res+(getmu(j)-getmu(i-1))*fpow(n/i,k))%MOD;
	}
	return res;
}

int main()
{
	findprm(N-1);
	scanf("%d%d",&n,&k);
	for (int i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans=(ans+(fib(j+2)-fib(i+1))*calc(n/i))%MOD;
	}
	cout<<(ans%MOD+MOD)%MOD;
	return 0;
}