题面
https://www.luogu.com.cn/problem/P3768
题解
前置知识
- 杜教筛 https://www.cnblogs.com/xh092113/p/12289949.html
- 求和中消gcd技巧 https://www.cnblogs.com/xh092113/p/12268356.html
看到gcd，开始反演
[{sumlimits_{i=1}^{n}}{sumlimits_{j=1}^{n}}ij(i,j)
]
[={sumlimits_{i=1}^{n}}{sumlimits_{j=1}^{n}}ij{sumlimits_{p|(i,j)}}{mu(p)}
]
[={sumlimits_{i=1}^{n}}{sumlimits_{j=1}^{n}}ij{sumlimits_{p|i,p|j}}{mu(p)}
]
[={sumlimits_{p}}{mu(p)}{sumlimits_{i=1}^{n}}{sumlimits_{j=1}^{n}}ij[p|i][p|j]
]
[={sumlimits_{p}}{mu(p)}*p^2*({frac{{lfloor}{frac{n}{p}}{
floor}({lfloor}{frac{n}{p}}{
floor}+1)}{2}})^2①
]
然后使用数论分块，剩下的问题就只是求解(f(p)={mu(p)}*p^2)的前缀和sf了。
令(g(p)=p^2，h(p)=p^3)（这里g,h的作用详见前置知识->杜教筛）就可以使用杜教筛求解了。值得注意的是，杜教筛可以一次性求出所有的(sf({lfloor}{frac{n}{i}}{ floor})(1{leq}i{leq}n))，而数论分块求①式的时候，也恰好只用到了这些sf值（分块时每一块的右端点R就是通过n/(n/L)得到的，所以每一个右端点都是({lfloor}{frac{n}{i}}{ floor}(1{leq}i{leq}n))中的一个），所以总体相当于只需要做一次杜教筛，就(O(n^{frac{2}{3}}))解决本题。
代码
- P.S.一定要注意爆long long的问题。虽然(p^2)并不会爆，但是(n^2)会爆。
- 如果此题的p再大一些，使得(p^2)也会爆的话，就必须要使用O(1)快速乘了。
#include<bits/stdc++.h>
using namespace std;
#define D 4641600
#define M 2160
#define ll long long
#define rg register
ll mod;
namespace ModCalc{
inline void Inc(ll &x,ll y){
x += y;if(x >= mod)x -= mod;
}
inline void Dec(ll &x,ll y){
x -= y;if(x < 0)x += mod;
}
inline ll Add(ll x,ll y){
Inc(x,y);return x;
}
inline ll Sub(ll x,ll y){
Dec(x,y);return x;
}
}
using namespace ModCalc;
ll sphi[D+M+5],pri[D/5+5];
ll n,pn,v2,v3;
bool isp[D+5];
inline ll F(ll x){
x = x % mod * ((x + 1) % mod) % mod * v2 % mod;
return x * x % mod;
}
inline void Eular(){
sphi[1] = 1;
for(rg ll i = 2;i <= D;i++)isp[i] = 1;
for(rg ll i = 2;i <= D;i++){
if(isp[i])pri[++pn] = i,sphi[i] = i - 1;
for(rg ll j = 1;i * pri[j] <= D;j++){
isp[i*pri[j]] = 0;
if(i % pri[j])sphi[i*pri[j]] = sphi[i] * sphi[pri[j]];
else{
sphi[i*pri[j]] = sphi[i] * pri[j];
break;
}
}
}
for(rg ll i = 1;i <= D;i++)sphi[i] = sphi[i] * i % mod * i % mod;
for(rg ll i = 2;i <= D;i++)Inc(sphi[i],sphi[i-1]);
}
inline ll id(ll x){
if(x <= D)return x;
return D + n / x;
}
inline ll sg(ll x){
return (x % mod) * ((x + 1) % mod) % mod * (((x << 1) | 1) % mod) % mod * v2 % mod * v3 % mod;
}
inline ll sumf(ll x){
if(sphi[id(x)])return sphi[id(x)];
ll ans = F(x),L,R = 1;
while(R < x){
L = R + 1,R = x / (x / L);
Dec(ans,Sub(sg(R),sg(L-1)) * sumf(x/L) % mod);
}
return sphi[id(x)] = ans;
}
int main(){
scanf("%lld%lld",&mod,&n);
Eular();
v2 = (mod + 1) >> 1;
v3 = mod % 3 == 1 ? ((mod << 1) | 1) / 3 : (mod + 1) / 3;
ll ans = 0,L,R = 0;
while(R < n){
L = R + 1;
R = n / (n / L);
Inc(ans,Sub(sumf(R),sumf(L-1)) * F(n/L) % mod);
}
cout << ans << endl;
return 0;
}