【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;

}

巴特西

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

题目描述

题解

代码

最新文章

热门文章

巴特西

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

题目描述

题解

代码

最新文章

热门文章

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