loj #6686. Stupid GCD 莫比乌斯反演+杜教筛

题目描述

传送门

求 (sumlimits_{i=1}^ngcd(sqrt[3]{i},i))

(n leq 10^{30})

分析

毒瘤题，读入和计算要用 (int128)

首先枚举 (sqrt[3]{i}) 的值

(sumlimits_{d=1}^{sqrt[3]{n}}sumlimits_{i=1}^n[sqrt[3]{i}=d]gcd(i,d))

方框里的那一个条件就等价于 (d^3 leq i <(d+1)^3)

式子变为 (sumlimits_{d=1}^{sqrt[3]{n}}sumlimits_{i=d^3}^{min((d+1)^3-1,n)}gcd(i,d))

令 (r=sqrt[3]{n})，把 (min) 拆掉

(sumlimits_{d=1}^{r-1}sumlimits_{i=d^3}^{(d+1)^3-1}gcd(i,d)+sumlimits_{i=r^3}^{n}gcd(r,i))

先考虑如何求出求 (sumlimits_{i=1}^ngcd(a,i))

(egin{aligned} sumlimits_{i=1}^ngcd(a,i) &=sumlimits_{d|a}dsumlimits_{i=1}^{n/d}[gcd(i,a/d)=1] \&= sumlimits_{d|a}dsumlimits_{i=1}^{n/d}sumlimits_{k|gcd(i,a/d)}mu(k)\&= sumlimits_{d|a}dsumlimits_{k|(a/d)}mu(k)frac{n}{kd}\&=sumlimits_{T|a}frac{n}{T}sumlimits_{d|T}mu(d)frac{T}{d}\ &=sumlimits_{T|a}frac{n}{T}varphi(T) end{aligned})

(egin{aligned}& sumlimits_{d=1}^{r-1}sumlimits_{i=d^3}^{(d+1)^3-1}gcd(i,d)+sumlimits_{i=r^3}^{n}gcd(r,i)\ &=sumlimits_{d=1}^{r-1}sumlimits_{T|d}varphi(T)(frac{(d+1)^3-1}{T}-frac{d^3-1}{T})+sumlimits_{T|r}varphi(T)(frac{n}{T}-frac{r^3-1}{T}) end{aligned})

后面的这一部分就可以枚举约数，对于每一个约数暴力计算 (varphi) 值

能过但是复杂度不太对

正确的做法是把质因子筛出来，然后进行 (dfs)，同时计算欧拉函数

因为暴力能过所以就没有去写

前面的这一部分可以继续化简

(egin{aligned}&sumlimits_{d=1}^{r-1}sumlimits_{T|d}varphi(T)(frac{(d+1)^3-1}{T}-frac{d^3-1}{T})\ &=sumlimits_{d=1}^{r-1}sumlimits_{T|d}varphi(T)(frac{(d+1)^3-1}{T}-frac{d^3-1}{T})\&=sumlimits_{T=1}^{r-1}varphi(T)sumlimits_{d=1}^{(r-1)/T}(frac{(dT+1)^3-1}{T}-frac{(dT)^3-1}{T})\&=sumlimits_{T=1}^{r-1}varphi(T)sumlimits_{d=1}^{(r-1)/T}3Td^2+3d+1end{aligned})

最后一步左边的部分是能整除 (T) 的，右边向下取整后会少一个 (1)，前面有一个减号又变成了加 (1)

这三项可以分别用杜教筛求解

(sumlimits_{T=1}^{r-1}Tvarphi(T)sumlimits_{d=1}^{(r-1)/T}3d^2+sumlimits_{T=1}^{r-1}varphi(T)sumlimits_{d=1}^{(r-1)/T}(3d+1))

筛一个 (varphi(n)) 再筛一个 (nvarphi(n)) 即可

因为洛谷 (IDE) 被禁了,而本机编译不了(int128)，所以推荐在这个网站编译

要注意的是如果用系统自带的 (pow) 函数开三次根号精度可能会出问题，所以要用两个循环判一下

代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#define rg register
const int maxn=8e6+5;
#define ll __int128
const int mod=998244353,inv2=499122177,inv6=166374059,Mod=4999999;
template <class Tp>
void read(Tp &x) {
	static char ch;
	static bool neg;
	for (ch = neg = 0; ch < '0' || ch > '9'; neg |= (ch == '-'), ch = getchar());
	for (x = 0; ch >= '0' && ch <= '9'; (x *= 10) += (ch ^ 48), ch = getchar());
	neg && (x = -x);
}
struct has{
	struct asd{
		int nxt,val;
		ll num;
	}b[maxn];
	int h[maxn],tot;
	void insert(rg ll num,rg int val){
		rg int now=num%Mod;
		b[tot].nxt=h[now];
		b[tot].val=val;
		b[tot].num=num;
		h[now]=tot++;
	}
	int cx(rg ll num){
		rg int now=num%Mod;
		for(rg int i=h[now];i!=0;i=b[i].nxt){
			if(b[i].num==num) return b[i].val;
		}
		return -1;
	}
}ans1,ans2;
int getsum1(rg ll now){
	now%=mod;
	return (now+1)*now%mod*inv2%mod;
}
int getsum2(rg ll now){
	now%=mod;
	return (2*now+1)%mod*(now+1)%mod*now%mod*inv6%mod;
}
bool not_pri[maxn];
int phi1[maxn],phi2[maxn],sum1[maxn],sum2[maxn],mmax,pri[maxn],ans;
void pre(){
	not_pri[0]=not_pri[1]=1;
	phi1[1]=1;
	for(rg int i=2;i<=mmax;i++){
		if(!not_pri[i]){
			pri[++pri[0]]=i;
			phi1[i]=i-1;
		}
		for(rg int j=1;j<=pri[0] && 1LL*i*pri[j]<=mmax;j++){
			not_pri[i*pri[j]]=1;
			if(i%pri[j]==0){
				phi1[i*pri[j]]=1LL*phi1[i]*pri[j]%mod;
				break;
			} else {
				phi1[i*pri[j]]=1LL*phi1[i]*phi1[pri[j]]%mod;
			}
		}
	}
	for(rg int i=1;i<=mmax;i++) phi2[i]=1LL*phi1[i]*i%mod;
	for(rg int i=1;i<=mmax;i++){
		sum1[i]=(sum1[i-1]+phi1[i])%mod;
		sum2[i]=(sum2[i-1]+phi2[i])%mod;
	}
}
int get_ans1(rg ll now){
	if(now<=mmax) return sum1[now];
	if(ans1.cx(now)!=-1) return ans1.cx(now);
	rg int nans=getsum1(now);
	for(rg ll l=2,r;l<=now;l=r+1){
		r=(now/(now/l));
		nans-=1LL*(r-l+1)*get_ans1(now/l)%mod;
		nans=(nans+mod)%mod;
	}
	ans1.insert(now,nans);
	return nans;
}
int get_ans2(rg ll now){
	if(now<=mmax) return sum2[now];
	if(ans2.cx(now)!=-1) return ans2.cx(now);
	rg int nans=getsum2(now);
	for(rg ll l=2,r;l<=now;l=r+1){
		r=(now/(now/l));
		nans-=1LL*(getsum1(r)-getsum1(l-1)+mod)%mod*get_ans2(now/l)%mod;
		nans=(nans+mod)%mod;
	}
	ans2.insert(now,nans);
	return nans;
}
ll n,sqr=1,tot;
int get_phi(rg ll now){
	if(now<=mmax) return phi1[now];
	rg ll cs=now,nans=now;
	for(rg ll i=2;i*i<=now;i++){
		if(cs%i==0){
			nans=nans/i*(i-1);
			while(cs%i==0) cs/=i;
			if(cs==1) return nans%mod;
		}
	}
	if(cs>1) nans=nans/cs*(cs-1);
	return nans%mod;
}
int solve(){
	rg ll now;
	rg int nans=0;
	for(rg ll i=1;i*i<=sqr;i++){
		if(sqr%i==0){
			now=i;
			nans=(nans+1LL*(n/now-(tot-1)/now)%mod*get_phi(now)%mod)%mod;
			if(i*i!=sqr){
				now=sqr/i;
				nans=(nans+1LL*(n/now-(tot-1)/now)%mod*get_phi(now)%mod)%mod;
			}
		}
	}
	return nans;
}
int main(){
	read(n);
	mmax=8e6;
	pre();
	sqr=std::pow(n,1.0/3);
	while((sqr+1)*(sqr+1)*(sqr+1)<=n) sqr++;
	while((sqr-1)*(sqr-1)*(sqr-1)>=n) sqr--;
	sqr--;
	rg ll cs,tmp;
	for(rg ll l=1,r;l<=sqr;l=r+1){
		r=(sqr/(sqr/l));
		cs=sqr/l;
		tmp=(get_ans1(r)-get_ans1(l-1)+mod)%mod;
		ans=(ans+1LL*(get_ans2(r)-get_ans2(l-1)+mod)%mod*3LL%mod*getsum2(cs)%mod)%mod;
		ans=(ans+3LL*tmp%mod*getsum1(cs)%mod)%mod;
		ans=(ans+tmp*cs%mod)%mod;
	}
	sqr++;
	tot=sqr*sqr*sqr;
	ans=(ans+solve())%mod;
	printf("%d
",ans);
	return 0;
}