P3768 简单的数学题

https://www.luogu.com.cn/problem/P3768
推式子+杜教筛+整除分块，感觉往 \(\varphi\) 的方向推比往 \(\mu\) 上推要好

\[\begin{aligned}\sum_{i=1}^n\sum_{j=1}^n ij\gcd(i,j) &=\sum_{i=1}^n\sum_{j=1}^n ij\sum_{k|\gcd(i,j)} \varphi(k)\\ &=\sum_{k}\varphi(k)\sum_{k|i}^n\sum_{k|j}^nij\\ &=\sum_{k}k^2\varphi(k)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i\sum_{j=1}^{\lfloor\frac{n}{j}\rfloor}j\\ &=\sum_{k}k^2\varphi(k)\left(\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i\right)^2\\ \end{aligned} \]

然后这个 \(\left(\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i\right)^2\) 可以整除分块并 \(O(1)\) 求，外面再套一个 \(k^2\varphi(k)\) 用杜教筛算就行了
至于怎么杜教筛，设 \(f=\varphi \cdot \operatorname{Id}^2,g=\operatorname{Id}^2\)，那么就有 \(f*g=\sum_{d|n} \varphi(d)d^2\left(\frac{n}{d}\right)^2=n^2\sum_{d|n} \varphi(d)=n^3\)
再设 \(s\) 是 \(f\) 的前缀和，根据杜教筛的套路式子，就有：

\[\begin{aligned}f(1)s(n) &=\sum_{i=1}^ng(i)s(\lfloor\frac{n}{i}\rfloor)-\sum_{i=2}^ng(i)s(\lfloor\frac{n}{i}\rfloor)\\ &=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)s(\lfloor\frac{n}{i}\rfloor)\\ \end{aligned} \]

由于 \(g(i)\) 以及 \(\sum_{i=1}^n(f*g)(i)\) 都是很好求的，所以就能直接杜教筛了

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<map>
#define reg register
#define LL_INF (long long)(0x3f3f3f3f3f3f3f3f)
#define INT_INF (int)(0x3f3f3f3f)
inline long long read(){
	register long long x=0;register int y=1;
	register char c=std::getchar();
	while(c<'0'||c>'9'){if(c=='-') y=0;c=getchar();}
	while(c>='0'&&c<='9'){x=x*10+(c^48);c=getchar();}
	return y?x:-x;
}
#define maxn 5000000
long long n;
long long mod,inv4,inv6;
std::map<long long,long long>map;
int notprime[maxn+6],prime[maxn+6];
long long phi[maxn+6];
inline long long power(long long a,long long b){
	long long ans=1;
	while(b){
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;b>>=1;
	}
	return ans;
}
inline void pre(int n){
	phi[1]=1;
	for(reg int i=2;i<=n;i++){
		if(!notprime[i]) prime[++prime[0]]=i,phi[i]=i-1;
		for(reg int x,j=1;i*prime[j]<=n&&j<=prime[0];j++){
			x=i*prime[j];notprime[x]=1;
			if(i%prime[j]) phi[x]=phi[i]*(prime[j]-1);
			else{
				phi[x]=phi[i]*prime[j];break;
			}
		}
		phi[i]=phi[i]*i%mod*i%mod;
	}
	for(reg int i=2;i<=n;i++) phi[i]=(phi[i]+phi[i-1])%mod;
}
inline long long S2(long long n){
	n%=mod;
	return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
inline long long S3(long long n){
	n%=mod;
	return n*n%mod*(n+1)%mod*(n+1)%mod*inv4%mod;
}
long long get(long long n){
	if(n<=maxn) return phi[n];
	if(map.find(n)!=map.end()) return map[n];
	long long ans=S3(n);
	for(reg long long l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		ans=(ans-get(n/l)*(S2(r)-S2(l-1))%mod+mod)%mod;
	}
	map[n]=ans;
	return ans;
}
inline long long work(long long n){
	long long ans=0;
	for(reg long long l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		ans+=S3(n/l)*(get(r)-get(l-1)+mod)%mod;
		ans%=mod;
	}
	return ans;
}
int main(){
	mod=read();n=read();
	pre(maxn);
	inv4=power(4,mod-2);inv6=power(6,mod-2);
	printf("%lld\n",work(n));
	return 0;
}