LOJ6682 梦中的数论

梦中的数论

有一天 LoliconAutomaton 梦到了下面这个问题：

\[\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}[(j\mid i) \land ((j+k)\mid i)] \]

其中，\((j\mid i) \land ((j+k)\mid i)\) 指 \(j\) 整除 \(i\) 并且 \(j+k\) 也整除 \(i\)。

但是 LoliconAutomaton 的数学实在是太差啦！你能帮一帮他吗？

对于全部数据，\(1\le n\le 10^{10}\)。

题解

http://jklover.hs-blog.cf/2020/04/19/Loj-6682-梦中的数论/#more

先枚举\(i\)，答案就变成了

\[ans=\sum_{i=1}^n \binom{\sigma_0(i)}{2}\\ =\frac {\sum_{i=1}^n \sigma_0^2(i)-\sum_{i=1}^n \sigma_0(i)} 2 \\ =\frac {\sum_{i=1}^n \sigma_0^2(i)-\sum_{i=1}^n \lfloor\frac{n}{i}\rfloor} 2 \]

后面那个\(\sum\)可以直接整除分块\(O(\sqrt n)\)求出。

对于前面那个\(\sum\)，注意到\(\sigma_0^2(x)\)也是积性函数，且当\(p\)为质数时，\(\sigma_0^2(p^k)=(k+1)^2\)，容易计算。

于是可以考虑用Min_25筛求出其前缀和，时间复杂度\(O(\frac{n^{\frac 3 4}}{\log n})\)。

CO int N=2e5+10,mod=998244353;
int n,m,prime[N],num;
int pos[N],idx,ref1[N],ref2[N],H[N];

int G(int x,int j){
	if(x<=1 or prime[j]>=x) return 0;
	int ans=(H[x<=m?ref1[x]:ref2[n/x]]-4*j)%mod;
	for(int k=j+1;k<=num and prime[k]*prime[k]<=x;++k){
		int pw=prime[k];
		for(int e=1;pw<=x;++e){
			ans=(ans+(e+1)*(e+1)*((e>1)+G(x/pw,k)))%mod;
			pw*=prime[k];
		}
	}
	return ans;
}
signed main(){
	read(n),m=ceil(sqrt(n));
	for(int i=2;i<=m;++i){
		if(!prime[i]) prime[++num]=i;
		for(int j=1;j<=num and i*prime[j]<=m;++j){
			prime[i*prime[j]]=1;
			if(i%prime[j]==0) break;
		}
	}
	int ans=0;
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=(ans-(n/l)%mod*(r-l+1))%mod;
		pos[++idx]=n/l;
		n/l<=m?ref1[n/l]=idx:ref2[n/(n/l)]=idx;
	}
	for(int i=1;i<=idx;++i) H[i]=(pos[i]-1)%mod;
	for(int j=1;j<=num;++j)
		for(int i=1;i<=idx and prime[j]*prime[j]<=pos[i];++i){
			int k=pos[i]/prime[j];
			k=k<=m?ref1[k]:ref2[n/k];
			H[i]=(H[i]-(H[k]-(j-1)))%mod;
		}
	for(int i=1;i<=idx;++i) H[i]=4*H[i]%mod;
	ans=(ans+G(n,0)+1)%mod;
	ans=(mod+1)/2*ans%mod;
	printf("%lld\n",(ans%mod+mod)%mod);
	return 0;
}