正题

题目链接:https://www.luogu.com.cn/problem/P3768

---

题目大意

给出$n,p$求$$\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)\ast i\ast j$$模$p$的值。

---

解题思路

下文中定义$H_y(x)={\sum }_{i=1}^x{i}^y$
首先显然是上莫反推式子$$f(x)=\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==x]\ast i\ast j=x^2\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{x}\rfloor }ij[gcd(i,j)==1]$$
$$F(x)=\sum _{x|d}f(d)=x^2\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{x}\rfloor }ij=x^2{H}_1(\lfloor \frac{n}{x}\rfloor {)}^2$$
然后答案为$$\sum _{x=1}^{n}f(x)=\sum _{x=1}^n{x}^2\sum _{x|d}F(d)\mu (\frac{x}{d})=\sum _{d=1}^{n}F(d)\sum _{x|d}\mu (\frac{x}{d})x^2$$
展开$F(x)$有
$$\sum _{d=1}^n{H}_1(\lfloor \frac{n}{d}\rfloor {)}^2{d}^2\sum _{x|d}\mu (\frac{x}{d})x^2$$
前面可以整除分块，问题变为了如何求$d^2{\sum }_{x|d}\mu (\frac{x}{d})x^2$的前缀和

因为$n$很大，所以我们考虑使用杜教筛。

先考虑一个式子$\mu \ast I=ϵ\Rightarrow \mu \ast I\ast \phi =ϵ\ast \phi \Rightarrow \mu \ast id=\phi$

所以上面那个东西就可以变成$$\sum _{d=1}^n{H}_1(\lfloor \frac{n}{d}\rfloor {)}^2{d}^2\phi (d)$$
上杜教筛$f(n)=\phi (n)n^2,g(n)=n^2$，$$f(n)g(n)=\sum _{d|n}f(d)g(\frac{n}{d})=\sum _{d|n}\phi (d)d^2\ast \frac{n^2}{d^2}=n^2\sum _{d|n}\phi (d)=n^3$$

然后有结论$H_3(n)=H_1(n{)}^2$，具体证明就用数学归纳法把，大致是多了一个$n+1$之后前者会多一个$(n+1{)}^3$，后者会多一个$2\ast H_1(n)\ast (n+1)+(n+1{)}^2=(n+1{)}^3$。这样就可以$O(n^{\frac{2}{3}})$的筛我们需要的东西了。

然后就是套莫反了，值得一提的是$H_2(n)=\frac{n\ast (n+1)\ast (2n+1)}{6}$，大体上也是数学归纳法证明吧。

理论上时间复杂度应该是$O(n^{\frac{2}{3}}\sqrt{n})$的，但是因为杜教筛的记忆化还有各种原因所以实际上其实跑的快很多了。

---

$code$

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#define ll long long
using namespace std;
const ll N=8e6;
ll n,p,cnt,inv2,inv6;
ll phi[N],pri[N];bool v[N];
map<ll,ll> mark;
void Prime(){phi[1]=1;for(ll i=2;i<N;i++){if(!v[i])pri[++cnt]=i,phi[i]=i-1;for(ll j=1;j<=cnt&&i*pri[j]<N;j++){v[i*pri[j]]=1;if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}phi[i*pri[j]]=phi[i]*(pri[j]-1);}}for(ll i=1;i<N;i++)phi[i]=(phi[i-1]+phi[i]%p*i%p*i%p)%p;return;
}
ll power(ll x,ll b){ll ans=1;while(b){if(b&1)ans=ans*x%p;x=x*x%p;b>>=1;}return ans;
}
ll s1(ll x)
{x%=p;return x*(x+1)%p*inv2%p;}
ll s2(ll x)
{x%=p;return x*(x+1)%p*(x+x+1)%p*inv6%p;}
ll Sum(ll x){if(x<N)return phi[x];if(mark[x])return mark[x];ll ans=s1(x);ans=ans*ans%p;for(ll l=2,r;l<=x;l=r+1){r=x/(x/l);(ans-=(s2(r)-s2(l-1))*Sum(x/l)%p)%=p;}return mark[x]=(ans+p)%p;
}
int main()
{scanf("%lld%lld",&p,&n);inv2=power(2,p-2);inv6=power(6,p-2);Prime();ll ans=0;for(ll l=1,r;l<=n;l=r+1){r=n/(n/l);ll tmp=s1(n/l);tmp=tmp*tmp%p;(ans+=(Sum(r)-Sum(l-1)+p)%p*tmp%p)%=p;}printf("%lld\n",ans);
}