Luogu P3768 简单的数学题

题目链接：Click here

Solution:

看到gcd，我们先尝试常规转化一下

\[\sum_{i=1}^n \sum_{j=1} ^n ij \,gcd(i,j)\\ \sum_{d=1}^n d^3 \sum_{i=1}^{\lfloor{n \over d}\rfloor} \sum_{j=1} ^{\lfloor{n \over d}\rfloor} ij [gcd(i,j)=1]\\ \sum_{d=1}^n d^3 \sum_{i=1}^{\lfloor{n \over d}\rfloor} \sum_{j=1} ^{\lfloor{n \over d}\rfloor} ij \sum_{t|i,t|j} \mu(t)\\ s(n)=\sum_{i=1}^n i\\ \sum_{d=1}^n d^3 \sum_{t=1} ^{n\over d} \mu(t) t^2s(\lfloor \frac{n}{dt} \rfloor )^2\\ T=dt\\ \sum_{T=1} ^n s(\lfloor \frac{n}{T} \rfloor) ^2\sum_{d|T} d^3\mu({T\over d}) {T \over d} ^2\\ \sum_{T=1} ^n s(\lfloor \frac{n}{T} \rfloor) ^2T^2\sum_{d|T} d\mu({T\over d}) \\ \]

做到这里，我们发现有点麻烦，前面可以数论分块，但后面部分线筛的时间我们是接受不了的

我们把后面那部分拿出来，看看怎么来处理

\[f(T)=T^2\sum_{d|T} d\mu({T\over d}) \\ \]

后面是一个卷积的形式，于是我们有

\[(id*\mu)=\varphi\\ f(n)=n^2\varphi(n)\\ \]

现在我们要求\(f(n)\)的前缀和，考虑杜教筛

\[S(n)=\sum_{i=1}^n f(i)\\ g(1)S(n)=\sum_{i=1}^n (g*f)(i)-\sum_{i=2}^n g(i)S(\lfloor{n\over i} \rfloor)\\ \]

我们让\(g(n)=n^2\)

\[(g*f)(n)=\sum_{d|n}f(d)g({n\over d})\\ (g*f)(n)=\sum_{d|n}d^2 \varphi(d){n^2 \over d^2}\\ (g*f)(n)=n^2\sum_{d|n}\varphi(d)\\ \]

我们知道\(\varphi(n)\)有很多美妙的性质

\[\sum_{d|n} \varphi(d)=n\\ (g*f)(n)=n^3\\ g(1)S(n)=\sum_{i=1}^n i^3-\sum_{i=2}^n i^2S(\lfloor{n\over i} \rfloor)\\ \]

对于幂和，我们很多结论

\[\sum_{i=1}^n i^3=s(n)^2\\ \sum_{i=1}^n i^2={n(n+1)(2n+1)\over 6}\\ \]

那么我们就可以利用杜教筛得到\(S(n)\)

\[\sum_{T=1} ^n s(\lfloor \frac{n}{T} \rfloor) ^2f(T)\\ \]

那么我们只需要对上面的式子数论分块即可

Code:

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=4e6+11;
bool vis[N];
int un,cnt,mod,phi[N],p[N];
int inv2,inv6,ans;
map<int,int> vphi;
int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-48;ch=getchar();}
	return x*f;
}
int qpow(int a,int b){
    int re=1;
    while(b){
        if(b&1) re=re*a%mod;
        b>>=1;a=a*a%mod;
    }return re;
}
void prepare(){phi[1]=1;
    for(int i=2;i<N;i++){
        if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*p[j]<N;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j]%mod;
                break;
            }phi[i*p[j]]=phi[i]*phi[p[j]]%mod;
        }
    }
    for(int i=1;i<N;i++){
        phi[i]=(i*i%mod)*phi[i]%mod;
        phi[i]=(phi[i]+phi[i-1])%mod;
    }
}
int sum(int x){x%=mod;return x*(x+1)%mod*inv2%mod;}
int sum2(int x){x%=mod;return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;}
int sum3(int x){x%=mod;return sum(x)*sum(x)%mod;}
int getS(int n){
    if(n<N) return phi[n];
    if(vphi[n]) return vphi[n];
    int re=sum3(n);
    for(int i=2,j;i<=n;i=j+1){
        j=n/(n/i);
        int tmp=(sum2(j)-sum2(i-1))%mod;
        re=re-(tmp*getS(n/i)%mod);re=re%mod;
    }return vphi[n]=(re+mod)%mod;
}
signed main(){
    mod=read(),un=read();
    prepare();
    inv2=qpow(2,mod-2);
    inv6=qpow(6,mod-2);
    for(int i=1,j;i<=un;i=j+1){
        j=un/(un/i);
        int tmp=getS(j)-getS(i-1);
        int uv=sum(un/i);uv=uv*uv%mod;
        tmp=tmp%mod;ans=ans+(uv*tmp%mod);
        ans=ans%mod;
    }printf("%lld\n",(ans+mod)%mod);
    return 0;
}