P3768 简单的数学题
给定 (n,p)，求 ：
[sum_{i=1}^nsum_{j=1}^nijgcd(i,j)quad (mod p)
]
(nle 10^{10},5 imes 10^{8} le ple 1.1 imes 10^{10})
话不多述直接开始推式子。
令 (gcd(i,j)=d)：
[egin{aligned}
&quadsum_{i=1}^nsum_{j=1}^nijgcd(i,j)\
&=sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^nij[gcd(i,j)=d]\
&=sum_{d=1}^nd^3sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[gcd(i,j)=1]
end{aligned}
]
设:
[f(x)=sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[gcd(i,j)=x]
]
[F(x)=sum_{x|i}f(i)
]
可以知道要求的是 (f(1))
莫比乌斯反演得到：
[egin{aligned}
&quad f(x)=sum_{x|i}mu(frac ix)F(i)\
&Rightarrow f(1)=sum_{i=1}^{n/d}mu(i)F(i)\
&ecause F(x)=sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[x|gcd(i,j)]\
&quad =x^2sum_{i=1}^{n/dx}sum_{j=1}^{n/dx}ij\
&quad =x^2(sum_{i=1}^{n/dx}i)^2\
& herefore Ans=sum_{d=1}^{n}d^3sum_{i=1}^{n/d}mu(i)i^2(sum_{j=1}^{n/id}j)^2
end{aligned}
]
设 (T=id,sum(x)=sumlimits_{i=1}^xi) ，则 (i=frac Td)，。
[egin{aligned}
Ans&=sum_{T=1}^nsum(frac nT)^2sum_{d|T}d^3(frac Td)^2mu(frac Td)\
&=sum_{T=1}^n sum(frac nT)^2cdot T^2sum_{d|T}dmu(frac Td)
end{aligned}
]
联想到 (mu *id=varphi)，我们可以得到：
[Ans=sum_{T=1}^nsum(frac nT)^2cdot T^2varphi(T)
]
前面的 (sum(frac nT)^2) 可以数论分块求，现在我们看看后面这块能不能求前缀和。
重新定义 (f(x)=x^2varphi(x))，我们发现它是一个积性函数。
先设 (S(n)=sumlimits_{i=1}^nf(i))，然后把杜教筛的式子摆在这里：
[g(1)S(n)=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(i)S(frac ni)
]
我们想要快速求 (sum f*g)。
先把 (sumlimits_{i=1}^n(f*g)(x)) 写开：
[sumlimits_{i=1}^n(f*g)(x)=sum_{i=1}^nsum_{d|i}g(frac id)d^2varphi(d)
]
我们考虑到 (varphi*I=id ext{，即} sumlimits_{d|n}varphi (d)=n)，就要将 (d^2) 因子消去，所以我们设 (g(x)=x^2)
[egin{aligned}
&quad sum_{i=1}^nsum_{d|i}g(frac id)d^2varphi(d)\
&=sum_{i=1}^ni^2sum_{d|i}varphi(d)\
&=sum_{i=1}^ni^3
end{aligned}
]
由小学奥数知识它的值为 (sum(n)^2)。
总之：
[S(n)=sum_{i=1}i^3-sum_{i=2}i^2S(frac ni)
]
答案是
[Ans=sum_{T=1}sum(frac nT)^2cdot T^2varphi(T)
]
数论分块后，前面半截可以直接求，后面半截可以用杜教筛在非线性时间里面处理出来。
时间在 (O(n^{frac 34})) 左右。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=8000010;
ll nn=N-10,mod;
ll pr[N],phi[N],tot=0;
bool mp[N];
unordered_map<ll,ll> f;
void init()
{
int n=nn;
phi[1]=mp[1]=1ll;
for(int i=2;i<=n;++i)
{
if(!mp[i]) pr[++tot]=i, phi[i]=i-1;
for(int j=1;j<=tot&&pr[j]*i<=n;++j)
{
ll tmp=1ll*i*pr[j];
mp[tmp]=1;
if(i%pr[j]) phi[tmp]=1ll*phi[i]*phi[pr[j]]%mod;
else {phi[tmp]=1ll*phi[i]*pr[j]%mod; break;}
}
}
for(int i=2;i<=n;++i)
phi[i]=(phi[i-1]+1ll*phi[i]*i%mod*i%mod)%mod;
}
ll inv2,inv6;
ll fpow(ll a,ll k,ll mod)
{
ll res=1; a%=mod;
while(k)
{
if(k&1) res=res*a%mod;
a=a*a%mod; k>>=1;
}
return (res+mod)%mod;
}
ll sum(ll x)
{
x%=mod;
return x*(x+1)%mod*inv2%mod;
}
ll sumP(ll x)
{
x%=mod;
return x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
}
ll GetSum(ll n)
{
if(n<=nn) return phi[n];
if(f[n]) return f[n];
ll res=sum(n);
res=res*res%mod;
for(ll l=2,r;l<=n;l=r+1)
{
r=n/(n/l);
ll T=(sumP(r)-sumP(l-1))%mod;
res-=GetSum(n/l)*T%mod;
res%=mod;
}
return f[n]=(res+mod)%mod;
}
int main()
{
#ifdef test
freopen("wa.txt","w",stdout);
#endif
ll n;
scanf("%lld%lld",&mod,&n);
nn=min(n,nn);
inv2=fpow(2,mod-2,mod),inv6=fpow(6,mod-2,mod);
init();
ll ans=0;
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ll Sum=sum(n/l);
Sum=Sum*Sum%mod;
ll T=(GetSum(r)-GetSum(l-1))%mod;
ans+=Sum*T%mod;
ans%=mod;
// cout<<GetSum(r)<<" "<<r<<"
";
}
// cout<<"-1 -1
";
printf("%lld",(ans+mod)%mod);
return 0;
}