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

题目描述

　　给你(n,p)，求

[sum_{i=1}^nsum_{j=1}^isum_{k=1}^igcd(i,j,k)mod p ]

　　(nleq {10}^9)

题解

[egin{align} ans&=sum_{i=1}^nsum_{j=1}^isum_{k=1}^igcd(i,j,k)\ &=sum_{i=1}^nsum_{j=1}^isum_{k=1}^isum_{d|gcd(i,j,k)}varphi(d)\ &=sum_{d=1}^nvarphi(d)sum_{d|i}^nsum_{d|j}^isum_{d|k}^i1\ &=sum_{i=1}^nvarphi(i)S_2(lfloorfrac{n}{i} floor)\ 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
",ans);
	return 0;
}