题解[Zeta的数论题[加强版]]
题意：
记 \(f(t)=\sum\limits_{k=1}^t k[\gcd(k,t)=1]\)
求 \(\sum\limits_{i=1}^n \sum\limits_{j=1}^n (i^2+j^2+ij)f(\gcd(i,j)),n\leq 10^{10}\)
按照原题的做法，将原式莫反，得到
\(\sum\limits_{T=1}^n \left(T^2\sum\limits_{d|T}^{}f(d)\mu(T/d)\right) \sum\limits_{i=1}^{\left\lfloor\frac{n}{T}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{n}{T}\right\rfloor}i^2+j^2+ij\)
以及还有 \(f(d)=\dfrac{1}{2}\left[d=1\right]+\dfrac{1}{2}d\varphi(d)\)
这仍然以数论分块做，但由于数据范围的扩大，需要快速求出\(g(T)=T^2\sum\limits_{d|T}f(d)\mu(T/d)\) 的前缀和。
将 \(f\) 的通项带入 \(g\) ，得 \(g(T)=\dfrac{1}{2}T^2\left(\mu(T)+\sum\limits_{d|T}d\varphi(d)\mu(T/d)\right)\)
可以拆为两部分，分别求前缀和。
设 \(r(T)=T^2\mu(T),h(T)=T^2\sum\limits_{d|T}d\varphi(d)\mu(T/d)\) ，则有 \(g(T)=\dfrac{1}{2}(r(T)+h(T))\)
首先做相对简单的 \(r(T)\) 的前缀和，这可以直接杜教筛：
\(\sum\limits_{d|n}r(d)\left(\frac{n}{d}\right)^2=n^2\sum\limits_{d|n}\mu(d)=n^2[n=1]\)
然后是 \(h(T)\) ，仍然可以杜教筛，但为了构造先看看其 \(\text{dgf}\)
\(H(z)=\sum\limits_{n\geq 1}\dfrac{1}{n^{z-2}}\sum\limits_{d|n}d\varphi(d)\mu(n/d)=\left(\sum\limits_{n\geq 1}\dfrac{\varphi(n)}{n^{z-3}}\right)\left(\sum\limits_{n\geq 1}\dfrac{\mu(n)}{n^{z-2}}\right)=\dfrac{\zeta(z-4)}{\zeta(z-3)\zeta(z-2)}\)
所以 \(H(z)\) 卷上 \(\zeta(z-3)\zeta(z-2)\) 能得到 \(\zeta(z-4)\)
而 \([n^{-z}]\zeta(z-4)=n^4\) 的前缀和是 \(\frac{1}{30}n(n+1)(2n+1)(3n^2+3n-1)\) 能 \(O(1)\) 求出。
剩下就是杜教筛时 \(\zeta(z-3)\zeta(z-2)\) 的前缀和。
设 \(s(n)=[n^{-z}]\zeta(z-3)\zeta(z-2)=\sum\limits_{d|n}d^3(\frac{n}{d})^2=n^2\sum\limits_{d|n}d\)
则 \(\sum\limits_{n=1}^{N}s(n)=\sum\limits_{n=1}^{N}n^2\sum\limits_{d|n}d=\sum\limits_{d=1}^{N}d\sum\limits_{1\leq n\leq N,d|n}n^2=\sum\limits_{d=1}^{n}d^3\sum\limits_{n=1}^{\left\lfloor\frac{n}{d}\right\rfloor} n^2\)
这可以杜教筛时再套一个数论分块算出来。
时间复杂度仍然是 \(O(n^{\frac{2}{3}})\) 的，因为套的这个数论分块，仍与接下来杜教筛的递归复杂度同阶。
代码：
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=4641588,mod=998244353,hmod=1e5+7;
const ll inv2=499122177,inv4=748683265,inv6=166374059;
const ll inv12=582309206,inv30=432572553;
int p[N+10],prm,mu[N+10];bool b[N+10];char ch;
ll d[N+10],pw[N+10],pp[N+10],sum[N+10],sum_[N+10],p2[N+10],ans,n;
void write(ll x){if(x>9)write(x/10);putchar(48+x%10);}
void pre_work(){
register int i,j,k,h;
mu[1]=d[1]=sum[1]=1;
for(i=2;i<=N;++i){
if(!b[i])p[++prm]=i,d[i]=pp[i]=i+1,p2[prm]=1ll*i*i%mod,sum[i]=(p2[prm]-i-1)%mod,mu[i]=-1;
for(j=1;j<=prm&&(k=i*p[j])<=N;++j){
b[k]=1;
if(i%p[j]==0){
pp[k]=pp[i]*p[j]+1;
d[k]=d[i]/pp[i]*pp[k];
sum[k]=sum[i]*p2[j]%mod;
if((h=i/p[j])%p[j])sum[k]=(sum[k]+sum[h]*p[j]%mod)%mod;
mu[k]=0;break;
}
mu[k]=-mu[i];
d[k]=d[i]*d[p[j]]%mod;
pp[k]=pp[p[j]];
sum[k]=sum[i]*sum[p[j]]%mod;
}
}
for(i=1;i<=N;++i)d[i]=(d[i-1]+d[i]*i%mod*i%mod)%mod;
for(i=1;i<=N;++i)sum_[i]=(sum_[i-1]+1ll*i*i%mod*mu[i]%mod)%mod;
for(i=1;i<=N;++i)sum[i]=(sum[i-1]+sum[i]*i%mod*i%mod)%mod;
}
inline ll s_1(ll x){x%=mod;return x*(x+1)%mod*inv2%mod;}
inline ll s_2(ll x){x%=mod;return x*(x+1)%mod*(x<<1|1)%mod*inv6%mod;}
inline ll s_3(ll x){x%=mod;return x*(x+1)%mod*(x*x%mod+x)%mod*inv4%mod;}
inline ll s_4(ll x){x%=mod;return x*(x+1)%mod*(x<<1|1)%mod*(3*x*x%mod+3*x-1)%mod*inv30%mod;}
inline ll g(ll x){x%=mod;return x*(x+1)%mod*(11*x%mod*x%mod+7*x%mod)%mod*inv12%mod;}
struct hash_table{
int h[N+10],nextn[N+10],edg;ll from[N+10],to[N+10];
void add(int s,ll x,ll y){to[++edg]=y;from[edg]=x;nextn[edg]=h[s];h[s]=edg;}
ll inquiry(ll x){
for(int i=h[x%hmod];i;i=nextn[i])if(from[i]==x)return to[i];
return 0;
}
}H,HH,H_;
ll getsum1(ll n){
if(n<=N)return d[n];
ll tmp=HH.inquiry(n);
if(tmp)return tmp;
tmp=0;
for(register ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
tmp=(tmp+(s_3(r)-s_3(l-1))*s_2(n/l)%mod)%mod;
}
HH.add(n%hmod,n,tmp);
return tmp;
}
ll getsum(ll n){
if(n<=N)return sum[n];
ll tmp=H.inquiry(n);
if(tmp)return tmp;
tmp=s_4(n);
for(register ll l=2,r;l<=n;l=r+1){
r=n/(n/l);
tmp=(tmp-(getsum1(r)-getsum1(l-1))*getsum(n/l)%mod)%mod;
}
H.add(n%hmod,n,tmp);
return tmp;
}
ll getsum_(ll n){
if(n<=N)return sum_[n];
ll tmp=H_.inquiry(n);
if(tmp)return tmp;
tmp=1;
for(register ll l=2,r;l<=n;l=r+1){
r=n/(n/l);
tmp=(tmp-(s_2(r)-s_2(l-1))*getsum_(n/l)%mod)%mod;
}
H_.add(n%hmod,n,tmp);
return tmp;
}
main(){
pre_work();ch=getchar();
while(ch>47)n=(n<<1)+(n<<3)+(ch^48),ch=getchar();
for(register ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans=(ans+(getsum(r)-getsum(l-1)+getsum_(r)-getsum_(l-1))*g(n/l)%mod)%mod;
}
ans=(ans*inv2%mod+mod)%mod;
write(ans);
}
