! P3768简单的数学题

(ans=(sum_{i=1}^nsum_{j=1}^nijgcd(i,j))\%p)

预处理(n^{frac 2 3})后，时间复杂度为(O(n^{frac23}))

#include<bits/stdc++.h>
using namespace std;
#define int long long 
inline int read(){
	int x=0,f=1;char c=getchar();
	while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}
	while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	return f==1?x:-x;
}
const int N=5e6+4;
int mod,tot,inv6,ans,s[N],pri[N];
bitset<N>b;
unordered_map<int,int>ms;
inline int ksm(int x,int r){
	int ret=1;
	for(int i=0;(1ll<<i)<=r;i++){
		if((r>>i)&1)ret=ret*x%mod;
		x=x*x%mod;
	}
	return ret;
}
inline void pre(int n){
	s[1]=1;
	for(int i=2;i<=n;i++){
		if(!b[i]){pri[++tot]=i;s[i]=i-1;}
		for(int j=1;j<=tot&&pri[j]*i<=n;j++){
			b[i*pri[j]]=1;
			if(i%pri[j])s[i*pri[j]]=s[i]*(pri[j]-1)%mod;
			else{
				s[i*pri[j]]=s[i]*pri[j]%mod;
				break;
			}
		}
	}
	for(int i=2;i<=n;i++)s[i]=(s[i-1]+s[i]*i%mod*i)%mod;
}
inline int sum(int x){
	x%=mod;
	return x*(x+1)/2%mod;
}
inline int sum2(int x){
	x%=mod;
	return x*(x+1)%mod*(x*2+1)%mod*inv6%mod;
}
inline int S(int n){
	if(n<N)return s[n];
	if(ms[n])return ms[n];
	int ret=sum(n)*sum(n)%mod;
	for(int l=2,r,x;l<=n;l=r+1){
		r=n/(n/l);
		x=sum2(r)-sum2(l-1);
		ret=(ret-x*S(n/l))%mod;
	}
	ms[n]=ret;
	return ret;
}
signed main(){
	mod=read();int n=read();
	pre(N-1);
	inv6=ksm(6,mod-2);
	for(int l=1,r,pr=0,nw;l<=n;l=r+1){
		r=n/(n/l);
		nw=S(r);
		ans=(ans+sum(n/l)*sum(n/l)%mod*(nw-pr))%mod;
		pr=nw;
	}
	cout<<(ans+mod)%mod;
	return (0-0);
}