LOJ6682 梦中的数论

梦中的数论

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

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

其中,((jmid i) land ((j+k)mid i))(j) 整除 (i) 并且 (j+k) 也整除 (i)

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

对于全部数据,(1le nle 10^{10})

题解

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

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

[ans=sum_{i=1}^n inom{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 lfloorfrac{n}{i} floor} 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
",(ans%mod+mod)%mod);
	return 0;
}
原文地址:https://www.cnblogs.com/autoint/p/12777400.html