测试地址:简单的数学题
做法:本题需要用到杜教筛+数论分块。
首先推式子,因为
∑d|nφ(d)=n
∑
d
|
n
φ
(
d
)
=
n
,所以有:
∑ni=1∑nj=1ijgcd(i,j)
∑
i
=
1
n
∑
j
=
1
n
i
j
gcd
(
i
,
j
)
=∑ni=1i∑nj=1j∑d|gcd(i,j)φ(d)
=
∑
i
=
1
n
i
∑
j
=
1
n
j
∑
d
|
gcd
(
i
,
j
)
φ
(
d
)
=∑nd=1φ(d)∑⌊nd⌋i=1id∑⌊nd⌋j=1jd
=
∑
d
=
1
n
φ
(
d
)
∑
i
=
1
⌊
n
d
⌋
i
d
∑
j
=
1
⌊
n
d
⌋
j
d
=∑nd=1φ(d)⋅d2⋅(⌊nd⌋(⌊nd⌋+1)2)2
=
∑
d
=
1
n
φ
(
d
)
⋅
d
2
⋅
(
⌊
n
d
⌋
(
⌊
n
d
⌋
+
1
)
2
)
2
显然如果我们能求
f(n)=φ(n)⋅n2
f
(
n
)
=
φ
(
n
)
⋅
n
2
的前缀和的话,就可以数论分块做了。由于
f
f
是积性函数,考虑杜教筛,令,那么有
(f∗g)(n)=n3
(
f
∗
g
)
(
n
)
=
n
3
,这两个函数都能
O(1)
O
(
1
)
算出前缀和,那么我们就解决了这一题。
要注意的一点是,因为数据范围很大,达到了
1010
10
10
,有的时候直接计算乘法会溢出,一定要在乘上一个
1010
10
10
大小的数前先对这个数取模。
以下是本人代码:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll hashsiz=8000009;
ll n,mod,limit,hashlist[hashsiz+10]={0},hashval[hashsiz+10];
ll phi[5000010],sum[5000010],prime[5000010],inv[7];
bool vis[5000010]={0};
void calc()
{
phi[1]=1;
prime[0]=0;
for(ll i=2;i<=limit;i++)
{
if (!vis[i])
{
prime[++prime[0]]=i;
phi[i]=i-1;
}
for(ll j=1;j<=prime[0]&&i*prime[j]<=limit;j++)
{
vis[i*prime[j]]=1;
if (i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
sum[0]=0;
for(ll i=1;i<=limit;i++)
sum[i]=(sum[i-1]+i*i%mod*phi[i]%mod)%mod;
}
ll sumg(ll n)
{
return (n%mod)*((n+1)%mod)%mod*((2*n+1)%mod)%mod*inv[6]%mod;
}
ll sumfg(ll n)
{
return (n%mod)*(n%mod)%mod*((n+1)%mod)%mod*((n+1)%mod)%mod*inv[4]%mod;
}
ll hashfind(ll x)
{
ll pos=x%hashsiz;
while(hashlist[pos]&&hashlist[pos]!=x) pos++;
if (hashlist[pos]) return pos;
else return -1;
}
void hashinsert(ll x,ll v)
{
ll pos=x%hashsiz;
while(hashlist[pos]&&hashlist[pos]!=x) pos++;
hashlist[pos]=x;
hashval[pos]=v;
}
ll solve(ll n)
{
ll pos=hashfind(n);
if (n<=limit) return sum[n];
if (pos!=-1) return hashval[pos];
ll ans=sumfg(n);
for(ll i=n;i>=2;i=n/(n/i+1))
{
ll l=max(n/(n/i+1)+1,2ll),r=i;
ans-=solve(n/i)*(sumg(r)-sumg(l-1));
ans=(ans+mod)%mod;
}
hashinsert(n,ans);
return ans;
}
ll F(ll x)
{
ll ans=(x%mod)*((x+1)%mod)%mod*inv[2]%mod;
ans=ans*ans%mod;
return ans;
}
int main()
{
scanf("%lld%lld",&mod,&n);
inv[1]=1;
for(ll i=2;i<=6;i++)
inv[i]=(mod-mod/i)*inv[mod%i]%mod;
for(ll i=1;i*i*i<=n;i++)
if ((i+1)*(i+1)*(i+1)>n)
{
limit=i*i;
break;
}
calc();
ll ans=0;
for(ll i=n;i>=1;i=n/(n/i+1))
{
ll l=n/(n/i+1)+1,r=i;
ans=((ans+F(n/i)*(solve(r)-solve(l-1)))%mod+mod)%mod;
}
printf("%lld",ans);
return 0;
}