【XSY2754】求和 莫比乌斯反演 杜教筛

题目描述

　　给你\(n,p\)，求

\[\sum_{i=1}^n\sum_{j=1}^i\sum_{k=1}^i\gcd(i,j,k)\mod p \]

　　\(n\leq {10}^9\)

题解

\[\begin{align} ans&=\sum_{i=1}^n\sum_{j=1}^i\sum_{k=1}^i\gcd(i,j,k)\\ &=\sum_{i=1}^n\sum_{j=1}^i\sum_{k=1}^i\sum_{d|\gcd(i,j,k)}\varphi(d)\\ &=\sum_{d=1}^n\varphi(d)\sum_{d|i}^n\sum_{d|j}^i\sum_{d|k}^i1\\ &=\sum_{i=1}^n\varphi(i)S_2(\lfloor\frac{n}{i}\rfloor)\\ \end{align} \]

　　其中\(S_2(n)=\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6}\)

　　直接杜教筛。

　　时间复杂度：\(O(n^\frac{2}{3})\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
ll p;
int n;
ll fp(ll a,ll b)
{
	ll s=1;
	for(;b;b>>=1,a=a*a%p)
		if(b&1)
			s=s*a%p;
	return s;
}
ll inv6;
ll S(ll n)
{
	return n*(n+1)%p*(2*n+1)%p*inv6%p;
}
ll g(ll x)
{
	return S(n/x);
}
int b[10000010];
int s[10000010];
int pri[1000010];
int cnt;
int b2[10000010];
ll s2[10000010];
const int maxn=10000000;
ll f(ll x)
{
	if(x<=maxn)
		return s[x];
	if(b2[n/x])
		return s2[n/x];
	b2[n/x]=1;
	ll res=x*(x+1)/2%p;
	for(int i=2,j;i<=x;i=j+1)
	{
		j=x/(x/i);
		res=(res-(j-i+1)*f(x/i))%p;
	}
	return s2[n/x]=res;
}
int main()
{
#ifndef ONLINE_JUDGE
	freopen("b.in","r",stdin);
	freopen("b.out","w",stdout);
#endif
	scanf("%d%lld",&n,&p);
	inv6=fp(6,p-2);
	s[1]=1;
	for(int i=2;i<=maxn;i++)
	{
		if(!b[i])
		{
			pri[++cnt]=i;
			s[i]=i-1;
		}
		for(int j=1;j<=cnt&&i*pri[j]<=maxn;j++)
		{
			b[i*pri[j]]=1;
			if(i%pri[j]==0)
			{
				s[i*pri[j]]=s[i]*pri[j];
				break;
			}
			s[i*pri[j]]=s[i]*s[pri[j]];
		}
	}
	for(int i=2;i<=maxn;i++)
		s[i]=(s[i-1]+s[i])%p;
	ll ans=0;
	memset(b2,0,sizeof b2);
	for(int i=1,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans=(ans+(f(j)-f(i-1))*g(i))%p;
	}
	ans=(ans+p)%p;
	printf("%lld\n",ans);
	return 0;
}