题意:
给你
n
n
n和
p
p
p,求
∑
i
=
1
n
∑
j
=
1
n
i
∗
j
∗
g
c
d
(
i
,
j
)
(
m
o
d
p
)
\sum_{i=1}^n\sum_{j=1}^ni*j*gcd(i,j)(mod\ p)
∑i=1n∑j=1ni∗j∗gcd(i,j)(mod p)。
n
<
=
1
e
10
,
p
n<=1e10,p
n<=1e10,p是一个大质数。
题解:
看到有各种奇怪的推法,但是看到有一种特别快速的推法。
∑
i
=
1
n
∑
j
=
1
n
i
∗
j
∗
g
c
d
(
i
,
j
)
\sum_{i=1}^n\sum_{j=1}^ni*j*gcd(i,j)
i=1∑nj=1∑ni∗j∗gcd(i,j) 根据
∑
d
∣
n
φ
(
d
)
=
n
\sum_{d|n}\varphi(d)=n
∑d∣nφ(d)=n得:
=
∑
i
=
1
n
∑
j
=
1
n
∑
k
∣
i
,
k
∣
j
φ
(
k
)
=\sum_{i=1}^n\sum_{j=1}^n\sum_{k|i,k|j}\varphi(k)
=i=1∑nj=1∑nk∣i,k∣j∑φ(k)
=
∑
k
=
1
n
∑
k
∣
i
∑
k
∣
j
i
∗
j
=\sum_{k=1}^n\sum_{k|i}\sum_{k|j}i*j
=k=1∑nk∣i∑k∣j∑i∗j 各提出一个
k
k
k来得:
=
∑
k
=
1
n
φ
(
k
)
∗
k
2
∑
i
=
1
n
/
k
∑
j
=
1
n
/
k
i
∗
j
=\sum_{k=1}^n\varphi(k)*k^2\sum_{i=1}^{n/k}\sum_{j=1}^{n/k}i*j
=k=1∑nφ(k)∗k2i=1∑n/kj=1∑n/ki∗j
=
∑
k
=
1
n
φ
(
k
)
∗
k
2
∑
i
=
1
n
/
k
i
∗
∑
j
=
1
n
/
k
j
=\sum_{k=1}^n\varphi(k)*k^2\sum_{i=1}^{n/k}i*\sum_{j=1}^{n/k}j
=k=1∑nφ(k)∗k2i=1∑n/ki∗j=1∑n/kj
=
∑
k
=
1
n
φ
(
k
)
∗
k
2
(
∑
i
=
1
n
/
k
i
)
2
=\sum_{k=1}^n\varphi(k)*k^2(\sum_{i=1}^{n/k}i)^2
=k=1∑nφ(k)∗k2(i=1∑n/ki)2
=
∑
k
=
1
n
φ
(
k
)
∗
k
2
∑
i
=
1
n
/
k
i
3
=\sum_{k=1}^n\varphi(k)*k^2\sum_{i=1}^{n/k}i^3
=k=1∑nφ(k)∗k2i=1∑n/ki3 这一步的证明请见这里。
然后这个式子的前半部分 ∑ k = 1 n φ ( k ) ∗ k 2 \sum_{k=1}^n\varphi(k)*k^2 ∑k=1nφ(k)∗k2是可以杜教筛的,后半部分是有公式可以 O ( 1 ) O(1) O(1)计算前缀和的,于是就可以做了。
至于前面那个式子怎么用杜教筛,其实我的学习笔记里有介绍,只不过是 ϕ ⋅ i d \phi\cdot id ϕ⋅id再乘了一个 i d id id就行了。请看这里。
于是就做完了,复杂度就是杜教筛的复杂度。
代码:
#include <bits/stdc++.h>
using namespace std;
long long n,mod,p[5000010],phi[5000010],s[5000010];
long long ans,ni;
const long long maxn=5000000;
int cnt,vis[5000010];
map<long long,long long> mp;
inline long long ksm(long long x,long long y)
{
long long res=1;
while(y)
{
if(y&1)
res=res*x%mod;
x=x*x%mod;
y>>=1;
}
return res;
}
inline long long calc2(long long x)
{
x%=mod;
return x*(x+1)%mod*(2*x+1)%mod*ni%mod;
}
inline long long calc3(long long x)
{
x%=mod;
return (x*(x+1)/2)%mod*((x*(x+1)/2)%mod)%mod;
}
inline long long solve(long long x)
{
if(x<=maxn)
return s[x];
if(mp[x])
return mp[x];
long long r,res=calc3(x);
for(long long l=2;l<=x;l=r+1)
{
r=x/(x/l);
res=(res-(calc2(r)-calc2(l-1)+mod)%mod*solve(x/l)%mod+mod)%mod;
}
mp[x]=res;
return res;
}
int main()
{
scanf("%lld%lld",&mod,&n);
phi[1]=1;
for(int i=2;i<=5000000;++i)
{
if(!vis[i])
{
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*p[j]<=5000000;++j)
{
vis[i*p[j]]=1;
if(i%p[j]==0)
{
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
for(int i=1;i<=5000000;++i)
s[i]=(s[i-1]+phi[i]*i%mod*i%mod)%mod;
ni=ksm(6,mod-2);
long long r;
for(long long l=1;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+(solve(r)-solve(l-1)+mod)%mod*calc3(n/l)%mod)%mod;
}
printf("%lld\n",ans);
return 0;
}