luogu P3768 简单的数学题

背景：

以下图片均来自我的$PDF$文件，谢绝转载。

题目传送门：

https://www.luogu.org/problemnew/show/P3768

题意：

思路：

代码：

#include<cstdio>
#include<cstring>
#include<map>
#include<algorithm>
#define LL long long
#define ma 7000000
#define R register
#define I inline
using namespace std;
map<LL,LL> w_phi;
	LL n,mod,inv6;
	bool bz[ma+10];
	int prime[ma+10];
	LL phi[ma+10],sum_phi[ma+10];
LL dg(LL x,LL k)
{
	if(!k) return 1;
	LL o=dg(x,k>>1);
	return k&1?o*o%mod*x%mod:o*o%mod;
}
I void init()
{
	inv6=dg(6,mod-2);
	bz[0]=bz[1]=true;
	phi[1]=1;
	int t=0;
	for(R int i=2;i<=ma;i++)
	{
		if(!bz[i])
		{
			prime[++t]=i;
			phi[i]=i-1;
		}
		for(R int j=1;j<=t&&prime[j]*i<=ma;j++)
		{
			bz[i*prime[j]]=true;
			if(!(i%prime[j]))
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
			phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
	}
	for(R int i=1;i<=ma;i++)
		sum_phi[i]=(sum_phi[i-1]+(LL)i*i%mod*phi[i]%mod)%mod;
}
I LL calc_sigma2(LL x)
{
	x%=mod;
	return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
I LL calc_sigma3(LL x)
{
	x%=mod;
	return ((x*(x+1)/2)%mod)*(x*(x+1)/2%mod)%mod;
}
LL calc_phi(LL x)
{
	if(x<=ma) return sum_phi[x];
	if(w_phi[x]) return w_phi[x];
	LL ans=calc_sigma3(x);
	for(R LL l=2,r;l<=x;l=r+1)
	{
		r=x/(x/l);
		ans=((ans-(LL)(calc_sigma2(r)-calc_sigma2(l-1)+mod)%mod*calc_phi(x/l)%mod)+mod)%mod;
	}
	return w_phi[x]=ans;
}
I LL solve(LL n)
{
	LL sum=0;
	for(R LL l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		sum=(sum+(calc_phi(r)-calc_phi(l-1)+mod)%mod*calc_sigma3(n/l)%mod)%mod;
	}
	return sum;
}
int main()
{
	scanf("%lld %lld",&mod,&n);
	init();
	printf("%lld",solve(n));
}