SP20173题解

膜拜 rqy。

题意：

求：

\[\sum_{i=1}^n \sigma_0(i^2) \]

首先我们知道 \(\sigma_0((p^k)^2)=2 \times k + 1=k+(k+1)=\sigma_0(p^k)+\sigma_0(p^{k-1})=(\mu^2 * \sigma_0)(p^k)\)

那么我们大力展开：

\[\sum_{i=1}^n\sum_{d|i}\mu^2(d) \times \sigma_0(\frac n d) \]

\[\sum_{d=1}^n \mu^2(d)\sum_{i=1}^{\lfloor \frac n d \rfloor}\sigma_0(\frac n d) \]

众所周知 \(\sum_{i=1}^n \sigma_0(i)=\sum_{i=1}^n \lfloor \frac n i \rfloor,\sum_{i=1}^n\mu^2(i)=\sum_{i=1}^{\sqrt n}\mu(i)\lfloor \frac n {i^2} \rfloor\)，那么这里再使用杜教筛的数据分治即可做到 \(O(n^{\frac 2 3}+Tn^{\frac 2 3})\)。

稍微卡下常数就能 11.53s 了。实际上容易发现瓶颈其实是线性筛。。。

#include<unordered_map>
#include<cstdio>
#include<bitset>
#include<cmath>
typedef unsigned ui;
typedef __uint128_t L;
typedef unsigned long long ull;
const ui M=1e8+5;
ui T,lim,d[M],s[M];ui top,pri[5761456];ull n[10005];std::bitset<M>vis,mu1,mu2;
std::unordered_map<ull,ull>s1,d1;
struct FastDiv{
	ull b,m;
	FastDiv(const ull&x=1):b(x),m(ull((L(1)<<64)/x)){}
	friend inline ull operator/(const ull&a,const FastDiv&mod){
		if(mod.b==1)return a;ull r=1+(L(mod.m)*a>>64);return r-(a-r*mod.b>>63);
	}
}F[1000005],G[1000005];
inline void sieve(const ui&n){
	ui i,j,x;d[1]=1;vis[1]=1;mu1[1]=1;
	for(i=1;i<=1000001;++i)F[i]=FastDiv(i),G[i]=FastDiv(1ull*i*i);
	for(i=2;i<=n;++i){
		if(!vis[i])mu2[pri[++top]=i]=true,d[i]=s[i]=2;
		for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
			if(vis[x]=true,!(i%pri[j])){
				d[x]=ui(d[i]/F[s[i]])*(s[x]=s[i]+1);break;
			}
			d[x]=d[i]<<1;s[x]=2;(mu1[i]||mu2[i])&&(mu1[x]=mu1[i]^1,mu2[x]=mu2[i]^1);
		}
	}
	for(i=1;i<=n;++i)s[i]=mu1[i]||mu2[i],s[i]+=s[i-1],d[i]+=d[i-1];
}
inline ull S(const ull&n){
	if(n<lim)return d[n];if(d1.find(n)!=d1.end())return d1[n];
	ui i;ull ans(0);for(i=1;1ull*i*i<=n;++i)ans+=n/F[i];return--i,d1[n]=(ans<<1)-1ull*i*i;
}
inline ull Sum(const ull&n){
	if(n<lim)return s[n];if(s1.find(n)!=s1.end())return s1[n];ui i;ull ans(0);
	for(i=1;1ull*i*i<=n;++i)mu1[i]&&(ans+=n/G[i]),mu2[i]&&(ans-=n/G[i]);return s1[n]=ans;
}
inline ull Solve(const ull&n){
	ull L,R,ans(0);ui x;
	for(L=1;1ull*lim*L<=n;++L)if(mu1[L]||mu2[L])ans+=S(n/F[L]);
	for(;1ull*L*L<=n;++L)if(mu1[L]||mu2[L])ans+=d[n/F[L]];
	for(x=n/F[L];L<=n;L=R+1)1ull*x*L>n&&--x,ans+=d[x]*(Sum(R=n/F[x])-Sum(L-1));
	return ans;
}
signed main(){
	ull mx(0);scanf("%u",&T);for(ui i=1;i<=T;++i)scanf("%llu",n+i),n[i]>mx&&(mx=n[i]);sieve(lim=pow(mx,2./3));
	for(ui i=1;i<=T;++i)printf("%llu\n",Solve(n[i]));
}