题目描述
给你\(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=1}^nd[\gcd(i,j,k)=d]\&=\sum_{d=1}^nd\sum_{d|l}\mu(\frac{l}{d})\sum_{i=1}^n\sum_{j=1}^i\sum_{k=1}^i[l|\gcd(i,j,k)]\&=\sum_{d=1}^nd\sum_{d|l}\mu(\frac{l}{d})S_2(\lfloor\frac{n}{l}\rfloor)\&=\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;
}