【2020杭电多校round1 1010】HDU6760 Math is Simple

题目大意

题目链接

(T)次询问。每次给定(n)。求

[sum_{substack{1leq a<bleq n\gcd(a,b)=1\a+bgeq n}}frac{1}{ab} ]

输出答案在(mod 998244353)意义下的值。

数据范围：(1leq Tleq 10^4)，(1leq nleq 10^8)。

本题题解

记(n)的答案为(f_n)。即：

[f_n=sum_{substack{1leq a<bleq n\gcd(a,b)=1\a+bgeq n}}frac{1}{ab} ]

考虑递推(f_n)。那么：

[f_n=f_{n-1}+sum_{substack{1leq a<bleq n\gcd(a,b)=1\a+b=n}}frac{1}{ab}-sum_{substack{1leq a<bleq n-1\gcd(a,b)=1\a+b=n-1}}frac{1}{ab} ]

发现后面两坨东西长得很像，只是(n)变成了(n-1)。于是我们设：

[g_n=sum_{substack{1leq a<bleq n\gcd(a,b)=1\a+b=n}}frac{1}{ab} ]

则上面的递推式可以写为：

[egin{align} f_n&=f_{n-1}+g_n-g_{n-1}\ &=(f_{n-2}+g_{n-1}-g_{n-2})+g_n-g_{n-1}=f_{n-2}-g_{n-2}+g_n\ &=dots\ &=f_{2}-g_2+g_n end{align} ]

又因为(f_2=frac{1}{2})，(g_{2}=0)（这可以手算出来），所以(f_n=frac{1}{2}+g_n)。

于是问题就转化为了求(g_n)。因为(a+b=n)。我们考虑只枚举(a)，则(b=n-a)。发现(gcd(a,n-a)=1)的充分必要条件是(gcd(a,n)=1)。再观察(frac{1}{ab})，可以写成(frac{1}{a(n-a)}=frac{1}{n}cdot frac{n}{a(n-a)}=frac{1}{n}cdot(frac{1}{a}+frac{1}{n-a}))。因此：

[g_n=frac{1}{n}sum_{substack{1leq aleq n\gcd(a,n)=1}}frac{1}{a} ]

这看起来比较简洁。因此很有希望推出正确的做法，我们继续推：

[egin{align} g_n&=frac{1}{n}sum_{a=1}^{n}frac{1}{a}[gcd(a,n)=1]\ &=frac{1}{n}sum_{a=1}^{n}frac{1}{a}sum_{d|gcd(a,n)}mu(d)\ &=frac{1}{n}sum_{d|n}mu(d)sum_{d|a}frac{1}{a}\ &=frac{1}{n}sum_{d|n}mu(d)sum_{i=1}^{frac{n}{d}}frac{1}{icdot d}\ &=frac{1}{n}sum_{d|n}mu(d)frac{1}{d}sum_{i=1}^{frac{n}{d}}frac{1}{i} end{align} ]

发现(sum_{i=1}^{m}frac{1}{i})这个东西是可以预处理出来的。记为(S(m))。则：

[g_n=frac{1}{n}sum_{d|n}mu(d)frac{1}{d}S(frac{n}{d}) ]

因为(n)的约数数量是(O(sqrt{n}))的。所以我们可以用(O(n+Tsqrt{n}))的时间复杂度解决本题。

但是因为(n)太大，空间也很卡。具体来说，空间限制是(512 ext{MB})，一个大小为(10^8)的( exttt{int})型数组，大小为(382 ext{MB})，所以我们只能开一个这样的数组。前面说过要预处理出(S)。那么(mu)就不能用线性筛预处理出来。我们可以先分解出(n)的所有质因子，然后用( ext{dfs})来枚举每个质因子的次幂，以此枚举出约数(d)，顺便求出(mu)。并且，(S)数组也只能开一个，那么有一种先求阶乘，再推逆元的做法就不管用了，必须利用下面这个公式求逆元：

[i^{-1}equiv -lfloorfrac{P}{i} floor(Pmod i)^{-1}pmod{P}quad(igeq2) ]

边界是(1^{-1}=1)。

参考代码：

//problem:HDU6760
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;

template<typename T>inline void ckmax(T& x,T y){x=(y>x?y:x);}
template<typename T>inline void ckmin(T& x,T y){x=(y<x?y:x);}

const int MOD=998244353;
inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
inline void add(int& x,int y){x=mod1(x+y);}
inline void sub(int& x,int y){x=mod2(x-y);}
inline int pow_mod(int x,int i){int y=1;while(i){if(i&1)y=(ll)y*x%MOD;x=(ll)x*x%MOD;i>>=1;}return y;}

const int MAXN = 1e8;
const int MAXC = 30;
int n,inv_2,inv_n;
int S[MAXN+5];
void init(){
	inv_2 = pow_mod(2,MOD-2);
	S[1]=1;
	for(int i=2;i<=MAXN;++i)
		S[i]=mod2( - (ll)(MOD/i) * S[MOD%i] % MOD);
	for(int i=2;i<=MAXN;++i)
		add(S[i],S[i-1]);
}
int p[MAXC+5],cnt,ans;
void dfs(int idx,int d,int mu){
	if(idx == cnt+1){
		mu=mod2(mu);
		add(ans,(ll)mu * mod2(S[d]-S[d-1]) %MOD * S[n/d] %MOD);
		return;
	}
	dfs(idx+1,d,mu);
	dfs(idx+1,d*p[idx],-mu);
	// 次数>=2时,mu=0,对答案没有贡献,不用递归
}
void solve_case(){
	cin>>n;
	if(n==2){
		cout<<inv_2<<endl;
		return;
	}
	inv_n = pow_mod(n,MOD-2);
	int tmp=n;
	cnt=0;
	for(int i=2;i*i<=n;++i)
		if(n%i==0){
			p[++cnt]=i;
			while(n%i==0)
				n/=i;
		}
	if(n!=1)
		p[++cnt]=n;
//	for(int i=1;i<=cnt;++i)
//		cerr<<p[i]<<" ";
//	cerr<<endl;
	n=tmp;
	ans=0;
	dfs(1,1,1);
	ans=(ll)ans*inv_n%MOD;
	add(ans,inv_2);
	cout<< ans <<endl;
}
int main() {
	init();
	int T;cin>>T;while(T--){
		solve_case();
	}
	return 0;
}

补充一下，最后推(g_n)那段，我们主要用了(sum_{d|n}mu(d)=[n=1])这个结论。这只能说是用到了莫比乌斯函数，而不是真正的莫比乌斯反演。

不过用反演也能推。设(F(x)=sum_{i=1}^{n}frac{1}{i}[gcd(i,n)=x])，则(g_n=frac{1}{n}F(1))，我们就是要求(F(1))。

设(G(x)=sum_{x|d}F(d))，可以推出(G(x)=[x|n]sum_{i=1}^{frac{n}{x}}frac{1}{ix})。

再套用莫比乌斯反演的结论，(F(x)=sum_{x|d}G(d)mu(frac{d}{x}))，可得(F(1)=sum_{d|n}mu(d)frac{1}{d}sum_{i=1}^{frac{n}{d}}frac{1}{i})。这与我们推出的式子是一样的，相当于一个验证。