[luogu 3768] 简单的数学题

一、题目

点此看题

二、解法

直接推式子，这几步你要有点莫比乌斯反演基础才行：

[sum_{i=1}^nsum_{j=1}^nijgcd(i,j) ]

[sum_{d=1}^nd^3sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[gcd(i,j)=1] ]

设 (sum(n)=sum_{i=1}^ni) ，我们使用莫比乌斯反演：

[sum_{d=1}^nd^3sum_{x=1}^{n/d}mu(x)x^2sum(frac{n}{dx})^2 ]

如果我们推到这一步就开始做的话，后面的东西是不好处理的，令 (T=id) ，可以把 (sum) 换到前面去：

[sum_{T=1}^nsum(frac{n}{T})^2T^2sum_{d|T}mu(frac{T}{d})d ]

因为基本性质：(mu*id=varphi) （这个乘是迪利克雷卷积的意思），所以可以把后面那部分换掉

[sum_{T=1}^nsum(frac{n}{T})^2T^2varphi(T) ]

然后就可以整除分块了，但是我们要知道 (f(x)=x^2varphi(x)) 的前缀和怎么算，不难发现这是一个杜教筛的问题。令辅助函数 (g(x)=x^2)，那么 ((f*g)(n)=sum_{x|n}x^2varphi(x)(frac{n}{x})^2=n^2sum_{x|n}varphi(x)=n^3)

然后 (g) 的前缀和也是很好算的，所以直接上杜教筛就可以了，给出一些基本的柿子（其实可以直接插值）：

[sum_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6} ]

[sum_{i=1}^ni^3=(frac{n(n+1)}{2})^2 ]

#include <cstdio>
#include <iostream>
#include <map>
using namespace std;
const int N = 8000000;
const int M = N+5;
#define int long long
int read()
{
	int x=0,f=1;char c;
	while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
	while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
int MOD,n,cnt,inv,ans,p[M],phi[M];map<int,int> sum;
void init(int n)
{
	phi[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!phi[i])
		{
			phi[i]=i-1;
			p[++cnt]=i;
		}
		for(int j=1;j<=cnt && i*p[j]<=n;j++)
		{
			if(i%p[j]==0)
			{
				phi[i*p[j]]=phi[i]*p[j];
				break;
			}
			phi[i*p[j]]=phi[i]*(p[j]-1);
		}
	}
	for(int i=2;i<=n;i++)
		phi[i]=(phi[i]*i%MOD*i+phi[i-1])%MOD;
}
int zy(int x)
{
	x%=MOD;
	return x*(x+1)/2%MOD;
}
int jzm(int x)
{
	x%=MOD;
	return x*(x+1)%MOD*(2*x+1)%MOD*inv%MOD;
}
int get(int n)
{
	if(n<=N) return phi[n];
	if(sum[n]) return sum[n];
	int ans=zy(n)*zy(n)%MOD;
	for(int l=2,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans=(ans-(jzm(r)-jzm(l-1))*get(n/l))%MOD;
	}
	return sum[n]=ans;
}
int qkpow(int a,int b)
{
	int r=1;
	while(b>0)
	{
		if(b&1) r=r*a%MOD;
		a=a*a%MOD;
		b>>=1;
	}
	return r;
}
signed main()
{
	MOD=read();n=read();
	init(N);inv=qkpow(6,MOD-2);
	for(int l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		ans=(ans+zy(n/l)*zy(n/l)%MOD*(get(r)-get(l-1)))%MOD;
	}
	printf("%lld
",(ans+MOD)%MOD);
}