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;
}