洛谷 3768简单的数学题（莫比乌斯反演+杜教筛）

题目链接

(Description)
求$$sum_{i-1}^nsum_{j=1}^nijgcd(i,j)mod p$$

(n<=10^{10})

---

(Solution)

[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] ]

发现跟这道题只差了一个(d^3)。最后化简得到

[ans=sum_{T=1}^n(frac{frac{n}{T} imes(frac{n}{T}+1)}{2})^2phi(T)T^2 ]

令(f(i)=phi(i)i^2)。
已知(phi*1=Id),所以可以令(g(i)=i^2),这样(f*g(x)=sum_{d|x}phi(d)d^2frac{x^2}{d^2}=x^3),(f*g)的前缀和就可以(O(1)计算了)。

#include<complex>
#include<cstdio>
#include<map>
#define LL long long
using namespace std;
const int N=7e6+7;
int tot,ans,mod,nn,div6;
LL n;
int prime[N],phi[N];
bool check[N];
map<LL,int>mp;
void Init()
{
	check[1]=phi[1]=1;
	nn=min(n,(LL)N-1);
	for(int i=2;i<=nn;i++)
	{
		if(!check[i])prime[++tot]=i,phi[i]=i-1;
		for(int j=1;j<=tot && i*prime[j]<=nn;j++)
		{
			check[i*prime[j]]=1;
			if(i%prime[j])phi[i*prime[j]]=phi[i]*phi[prime[j]];
			else
			{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
		}
	}
	for(int i=1;i<=nn;i++)
		phi[i]=(phi[i-1]+1ll*i*i%mod*phi[i]%mod)%mod;
}
int calc1(LL x)
{
	x%=mod;
	return x*(x+1)/2%mod;
}
int calc2(LL x)
{
	x%=mod;
	return x*(x+1)%mod*(x+x+1)%mod*div6%mod;
}
LL Sum(LL x)
{
	if(x<=nn)return phi[x];
	if(mp[x])return mp[x];
	LL res=calc1(x);res=res*res%mod;
	for(LL l=2,r;l<=x;l=r+1)
	{
		r=x/(x/l);
		res=(res-(calc2(r)-calc2(l-1))*Sum(x/l)%mod)%mod;
	}
	return mp[x]=(res+mod)%mod;
}
int Fpow(LL b,int p)
{
	LL res=1;
	for(;p;p>>=1,b=b*b%mod)
		if(p&1)res=res*b%mod;
	return res;
}
int main()
{
	scanf("%d%lld",&mod,&n);
	div6=Fpow(6,mod-2);
	Init();
	for(LL l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		LL tmp=calc1(n/l);
		ans=(ans+tmp*tmp%mod*(Sum(r)-Sum(l-1))%mod)%mod;
	}
	printf("%d
",(ans+mod)%mod);
	return 0;
}