洛谷P3768 简单的数学题 杜教筛

看样子就是个数学题。
其实这篇博客也是参考了部分题解的。
我们要求的是
$$ans=\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$$
这个肯定需要进行相应的变换
$$ans=\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^{n}ij[gcd(i,j)=d]（枚举gcd）$$
$$=\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor n/d\rfloor }\sum _{j=1}^{\lfloor n/d\rfloor }ij[gcd(i,j)=1](把d提取出来)$$
$$=\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor n/d\rfloor }\sum _{j=1}^{\lfloor n/d\rfloor }ij\sum _{x|gcd(i,j)}\mu (x)(莫比乌斯反演)$$
$$=\sum _{d=1}^n{d}^3\sum _{x=1}^{\lfloor n/d\rfloor }\mu (x)x^2\sum _i^{\lfloor n/dx\rfloor }\sum _j^{\lfloor n/dx\rfloor }ij(把x提前)$$
$$=\sum _{d=1}^n{d}^3\sum _{x=1}^{\lfloor n/d\rfloor }\mu (x)x^2(1+2+...+\lfloor \frac{n}{xd}\rfloor {)}^2(后面两个相等，可以自己理解一下)$$
我们令T=xd, Sum(x)=${\sum }_{i=1}^{x}i$
$$ans=\sum _{d=1}^n{d}^3\sum _{x=1}^{\lfloor n/d\rfloor }\mu (x)x^{2}Sum(\lfloor \frac{n}{T}\rfloor )$$
$$=\sum _{T=1}^{n}Su{m}^2(\lfloor \frac{n}{T}\rfloor )\sum _{d|T}{d}^3\ast (\frac{T}{d}{)}^2\ast \mu (\frac{T}{d})$$
$$=\sum _{T=1}^{n}Su{m}^2(\lfloor \frac{n}{T}\rfloor )\ast T^2\sum _{d|T}d\ast \mu (\frac{T}{d})$$
由于${\sum }_{d|T}d\ast \mu (T/d)=\varphi (T)$(不懂的可以看这里），因此我们有
$$ans=\sum _{T=1}^{n}Su{m}^2(\lfloor \frac{n}{T}\rfloor )\ast T^2\varphi (T)$$
由杜教筛的结论，我们可以对前面对部分采用数论分块，并可以在低于线性时间内求出$\varphi (T)$的前缀和，这题就这么做完了。
由于本人比较菜，并不会算复杂度，数论分块的复杂度是$\sqrt{n}$，杜教筛复杂度为$n^{2/3}$,最后复杂度应该是低于线性的吧

附上代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6+5;
ll pri[N],vis[N],phi[N],tot,mod,inv2,inv6;
unordered_map<ll,ll>vv,mp;
ll power(ll a,ll n){
	ll res=1;
	for(;n;n>>=1,a=a*a%mod)
		if(n&1)
			res=res*a%mod;
	return res;
}
void init(){
	inv2=power(2,mod-2);
	inv6=power(6,mod-2);
	phi[1]=1;
	for(int i=2;i<N;i++){
		if(!vis[i]){
			pri[tot++]=i;
			phi[i]=i-1;
		}
		for(int j=0;j<tot&&pri[j]*i<N;j++){
			int x=i*pri[j];
			vis[x]=1;
			if(i%pri[j]==0){
				phi[x]=phi[i]*pri[j];
				break;
			}
			phi[x]=phi[i]*(pri[j]-1);
		}
	}
	for(int i=2;i<N;i++){
		phi[i]=(phi[i-1]+1ll*i*i%mod*phi[i]%mod)%mod;
	}
}
ll sum1(ll x){
	x%=mod;
	return x%mod*(x+1)%mod*inv2%mod; 
}
ll sum2(ll x){
	x%=mod;
	return x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
}
ll cal(ll x){
	if(x<N)
		return phi[x];
	if(vv[x])
		return mp[x];
	ll res=0;
	for(ll i=2,j;i<=x;i=j+1){
		j=x/(x/i);
		res=(res+cal(x/i)*(sum2(j)-sum2(i-1)+mod)%mod)%mod;
	}
	res=(sum1(x)*sum1(x)%mod-res+mod)%mod;
	vv[x]=1,mp[x]=res;
	return res;
}
int main(){
	ll n,ans=0;
	scanf("%lld%lld",&mod,&n);
	init();
	for(ll i=1,j;i<=n;i=j+1){
		j=n/(n/i);
		ans=(ans+sum1(n/i)*sum1(n/i)%mod*(cal(j)-cal(i-1)+mod)%mod)%mod;
	}
	printf("%lld\n",ans);
	return 0;
}