做题重心转移到 LOJ 了。
至于为什么，如果你知道“……”的密码，就去看吧。
LOJ 上用户自创题大多数都不可做，今天看到个可做题（而且还是个水题），就来做了一发。
明显枚举立方根。（以下令 $m=lfloorsqrt[3]{n} floor$）
$$sumlimits_{i=1}^msumlimits_{j=i^3}^{min(n,(i+1)^3-1)}gcd(i,j)$$
由于 $i=m$ 比较特殊，我们把它拎出来：（其实就是把 $min$ 拆开）
$$sumlimits_{i=1}^{m-1}sumlimits_{j=i^3}^{(i+1)^3-1}gcd(i,j)+sum_{j=m^3}^ngcd(m,j)$$
现在考虑 $sumlimits_{i=1}^ngcd(x,i)$ 怎么求。
$$sumlimits_{i=1}^nsumlimits_{d|(x,i)}varphi(d)$$
$$sumlimits_{d|x}varphi(d)lfloorfrac{n}{d} floor$$
那么套回去：
$$sumlimits_{i=1}^{m-1}sumlimits_{d|i}varphi(d)(lfloorfrac{(i+1)^3-1}{d} floor-lfloorfrac{i^3-1}{d} floor)+sumlimits_{d|m}varphi(d)(lfloorfrac{n}{d} floor-lfloorfrac{m^3-1}{d} floor)$$
后面那部分可以直接枚举因数，暴力求欧拉函数。我这里预处理了前 $10^7approx m^{frac{2}{3}}$ 个欧拉函数，时间复杂度是远远小于 $O(m^{frac{2}{3}})$ 的。
接下来继续推前面的：
$$sumlimits_{i=1}^{m-1}sumlimits_{d|i}varphi(d)(lfloorfrac{(i+1)^3-1}{d} floor-lfloorfrac{i^3-1}{d} floor)$$
$$sumlimits_{d=1}^{m-1}varphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}(lfloorfrac{(id+1)^3-1}{d} floor-lfloorfrac{(id)^3-1}{d} floor)$$
$$sumlimits_{d=1}^{m-1}varphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}(lfloorfrac{(id)^3+3(id)^2+3(id)}{d} floor-lfloorfrac{(id)^3-1}{d} floor)$$
$$sumlimits_{d=1}^{m-1}varphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}((i^3d^2+3i^2d+3i)-(lfloor i^3d^2-frac{1}{d} floor))$$
$$sumlimits_{d=1}^{m-1}varphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}((i^3d^2+3i^2d+3i)-(i^3d^2-1))$$
$$sumlimits_{d=1}^{m-1}varphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}(3i^2d+3i+1)$$
$$sumlimits_{d=1}^{m-1}dvarphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}3i^2+sumlimits_{d=1}^{m-1}varphi(d)sumlimits_{i=1}^{lfloorfrac{m-1}{d} floor}(3i+1)$$
分块杜教筛即可。
时间复杂度：由于每次分块杜教筛都跑得过 $10^{10}$，我就姑且当它是 $O(m^{frac{2}{3}})$ 的。所以时间复杂度 $O(n^{frac{2}{9}})$。
upd: 已经会证复杂度是 $O(m^{frac{2}{3}})$ 了。如果用 map 是可以做到这个复杂度（忽略那个 $log$ 啦），但如果用其他什么奇怪的记忆化方法复杂度可能退化到 $O(m^{frac{5}{6}})$。
代码很丑。
#include<bits/stdc++.h> using namespace std; typedef long long ll; typedef __int128 lll; const int maxn=10001000,mod=998244353,inv6=166374059; #define FOR(i,a,b) for(int i=(a);i<=(b);i++) #define ROF(i,a,b) for(int i=(a);i>=(b);i--) #define MEM(x,v) memset(x,v,sizeof(x)) inline lll read(){ char ch=getchar();lll x=0,f=0; while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar(); while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar(); return f?-x:x; } lll n,nn; int pr[maxn],pl,phi[maxn],pre1[maxn],pre2[maxn],ans; bool vis[maxn]; map<ll,int> p1,p2; ll cb,cb1; void init(int upr){ phi[1]=1; FOR(i,2,upr){ if(!vis[i]) phi[pr[++pl]=i]=i-1; FOR(j,1,pl){ if(i*pr[j]>upr) break; vis[i*pr[j]]=true; if(i%pr[j]) phi[i*pr[j]]=phi[i]*(pr[j]-1); else{ phi[i*pr[j]]=phi[i]*pr[j];break; } } } pre1[1]=pre2[1]=1; FOR(i,2,upr){ pre1[i]=(pre1[i-1]+phi[i])%mod; pre2[i]=(pre2[i-1]+1ll*phi[i]*i)%mod; } } int s1(ll n){ int nn=n%mod; return 1ll*nn*(nn+1)/2%mod; } int s1(ll l,ll r){return (s1(r)-s1(l-1)+mod)%mod;} int s2(ll n){ int nn=n%mod; return 1ll*nn*(nn+1)%mod*(2*nn+1)%mod*inv6%mod; } int s2(ll l,ll r){return (s2(r)-s2(l-1)+mod)%mod;} int S1(ll n){ if(n<=10000000) return pre1[n]; if(p1.count(n)) return p1[n]; int nn=n%mod,s=1ll*nn*(nn+1)/2%mod; for(ll l=2,r;l<=n;l=r+1){ r=n/(n/l); s=(s-1ll*(r-l+1)%mod*S1(n/l)%mod+mod)%mod; } return p1[n]=s; } int S2(ll n){ if(n<=10000000) return pre2[n]; if(p2.count(n)) return p2[n]; int nn=n%mod,s=1ll*nn*(nn+1)%mod*(2*nn+1)%mod*inv6%mod; for(ll l=2,r;l<=n;l=r+1){ r=n/(n/l); s=(s-1ll*s1(l,r)*S2(n/l)%mod+mod)%mod; } return p2[n]=s; } int calc_phi(ll n){ if(n<=10000000) return phi[n]; ll s=n; for(ll i=2;i*i<=n;i++) if(n%i==0){ s=s/i*(i-1); while(n%i==0) n/=i; } if(n>1) s=s/n*(n-1); return s%mod; } int calc(ll x){ int s=0; for(ll i=1;i*i<=x;i++) if(x%i==0){ s=(s+1ll*calc_phi(i)*((n/i)%mod-(nn/i)%mod+mod)%mod)%mod; if(i*i!=x) s=(s+1ll*calc_phi(x/i)*((n/(x/i))%mod-(nn/(x/i))%mod+mod)%mod)%mod; } return s; } int main(){ n=read(); cb=pow(n,1.0/3);cb1=cb-1; nn=(lll)cb*cb*cb-1; init(min(cb,10000000ll)); for(ll l=1,r;l<=cb1;l=r+1){ r=cb1/(cb1/l); ans=(ans+3ll*s2(1,cb1/l)*(S2(r)-S2(l-1)+mod)%mod+(3ll*s1(1,cb1/l)+cb1/l)%mod*(S1(r)-S1(l-1)+mod)%mod)%mod; } printf("%d ",(ans+calc(cb))%mod); }
upd：
被出题人针对着卡，结果真卡掉了……
其实是开立方根出精度问题了，改一改就好了。
#include<bits/stdc++.h> using namespace std; typedef long long ll; typedef __int128 lll; const int maxn=10001000,mod=998244353,inv6=166374059; #define FOR(i,a,b) for(int i=(a);i<=(b);i++) #define ROF(i,a,b) for(int i=(a);i>=(b);i--) #define MEM(x,v) memset(x,v,sizeof(x)) inline lll read(){ char ch=getchar();lll x=0,f=0; while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar(); while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar(); return f?-x:x; } lll n,nn; int pr[maxn],pl,phi[maxn],pre1[maxn],pre2[maxn],ans; bool vis[maxn]; map<ll,int> p1,p2; ll cb,cb1; void init(int upr){ phi[1]=1; FOR(i,2,upr){ if(!vis[i]) phi[pr[++pl]=i]=i-1; FOR(j,1,pl){ if(i*pr[j]>upr) break; vis[i*pr[j]]=true; if(i%pr[j]) phi[i*pr[j]]=phi[i]*(pr[j]-1); else{ phi[i*pr[j]]=phi[i]*pr[j];break; } } } pre1[1]=pre2[1]=1; FOR(i,2,upr){ pre1[i]=(pre1[i-1]+phi[i])%mod; pre2[i]=(pre2[i-1]+1ll*phi[i]*i)%mod; } } int s1(ll n){ int nn=n%mod; return 1ll*nn*(nn+1)/2%mod; } int s1(ll l,ll r){return (s1(r)-s1(l-1)+mod)%mod;} int s2(ll n){ int nn=n%mod; return 1ll*nn*(nn+1)%mod*(2*nn+1)%mod*inv6%mod; } int s2(ll l,ll r){return (s2(r)-s2(l-1)+mod)%mod;} int S1(ll n){ if(n<=10000000) return pre1[n]; if(p1.count(n)) return p1[n]; int nn=n%mod,s=1ll*nn*(nn+1)/2%mod; for(ll l=2,r;l<=n;l=r+1){ r=n/(n/l); s=(s-1ll*(r-l+1)%mod*S1(n/l)%mod+mod)%mod; } return p1[n]=s; } int S2(ll n){ if(n<=10000000) return pre2[n]; if(p2.count(n)) return p2[n]; int nn=n%mod,s=1ll*nn*(nn+1)%mod*(2*nn+1)%mod*inv6%mod; for(ll l=2,r;l<=n;l=r+1){ r=n/(n/l); s=(s-1ll*s1(l,r)*S2(n/l)%mod+mod)%mod; } return p2[n]=s; } int calc_phi(ll n){ if(n<=10000000) return phi[n]; ll s=n; for(ll i=2;i*i<=n;i++) if(n%i==0){ s=s/i*(i-1); while(n%i==0) n/=i; } if(n>1) s=s/n*(n-1); return s%mod; } int calc(ll x){ int s=0; for(ll i=1;i*i<=x;i++) if(x%i==0){ s=(s+1ll*calc_phi(i)*((n/i)%mod-(nn/i)%mod+mod)%mod)%mod; if(i*i!=x) s=(s+1ll*calc_phi(x/i)*((n/(x/i))%mod-(nn/(x/i))%mod+mod)%mod)%mod; } return s; } int main(){ n=read(); cb=pow(n,1.0/3); while((lll)cb*cb*cb<=n) cb++; while((lll)cb*cb*cb>n) cb--; cb1=cb-1; nn=(lll)cb*cb*cb-1; init(min(cb,10000000ll)); for(ll l=1,r;l<=cb1;l=r+1){ r=cb1/(cb1/l); ans=(ans+3ll*s2(1,cb1/l)*(S2(r)-S2(l-1)+mod)%mod+(3ll*s1(1,cb1/l)+cb1/l)%mod*(S1(r)-S1(l-1)+mod)%mod)%mod; } printf("%d ",(ans+calc(cb))%mod); }