1. 程式人生 > 實用技巧 >loj #6686. Stupid GCD 莫比烏斯反演+杜教篩

loj #6686. Stupid GCD 莫比烏斯反演+杜教篩

題目描述

傳送門

\(\sum\limits_{i=1}^ngcd(\sqrt[3]{i},i)\)

\(n \leq 10^{30}\)

分析

毒瘤題,讀入和計算要用 \(int128\)

首先列舉 \(\sqrt[3]{i}\) 的值

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

方框裡的那一個條件就等價於 \(d^3 \leq i <(d+1)^3\)

式子變為 \(\sum\limits_{d=1}^{\sqrt[3]{n}}\sum\limits_{i=d^3}^{min((d+1)^3-1,n)}gcd(i,d)\)

\(r=\sqrt[3]{n}\),把 \(min\) 拆掉

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

先考慮如何求出求 \(\sum\limits_{i=1}^ngcd(a,i)\)

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

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

後面的這一部分就可以暴力去計算了

前面的這一部分可以繼續化簡

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

這三項可以分別用杜教篩求解

\(\sum\limits_{T=1}^{r-1}T\varphi(T)\sum\limits_{d=1}^{(r-1)/T}3d^2+\sum\limits_{T=1}^{r-1}\varphi(T)\sum\limits_{d=1}^{(r-1)/T}(3d+1)\)

篩一個 \(\varphi(n)\) 再篩一個 \(n\varphi(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\n",ans);
	return 0;
}