I.简单的数学题
在做这题之前，我们先来见一位老朋友：
\(\sum\limits_{i=1}^n\sum\limits_{j=1}^n\gcd(i,j)\)
我们在1.v.[NOI2010]能量采集中就已经接触到了这道题。当时我们运用了\(\sum\limits_{d|n}\mu(d)=[n=1]\)的公式。现在，我们要运用另一个公式\(\sum\limits_{d|n}\varphi(d)=n\)。
\(\begin{aligned}&\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\\=&\sum_{i=1}^n\sum_{j=1}^n\sum_{x|\gcd(i,j)}\varphi(x)\\=&\sum_{i=1}^n\sum_{j=1}^n\sum_{x|i,x|j}\varphi(x)\\=&\sum_{x=1}^n\varphi(x)\left\lfloor\dfrac{n}{x}\right\rfloor^2\end{aligned}\)
前面的东西可以线性筛前缀和，后面的东西可以整除分块，这样我们也得到一种和之前运用的\(\mu\)复杂度一致的运算。
而现在这道题呢？
我们要求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\gcd(i,j)\)。
再来！
\(\begin{aligned}&\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{x|\gcd(i,j)}\varphi(x)\\=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{x|i,x|j}\varphi(x)\\=&\sum_{x=1}^n\varphi(x)x^2\sum_{i=1}^{n/x}\sum_{j=1}^{n/x}ij\\=&\sum_{x=1}^n\varphi(x)x^2(\sum_{i=1}^{n/x}i)^2\end{aligned}\)
后面的\((\sum_{i=1}^{n/x}i)^2\)可以直接套等差数列求和公式并整除分块。前面的我们考虑杜教筛。
设\(f(x)=\varphi(x)x^2\)，我们要求\(S(n)=\sum_{i=1}^n\varphi(i)i^2\)。
我们想找一个合适的\(g\)。
考虑\((f*g)(n)=\sum\limits_{d|n}\varphi(d)d^2g(\dfrac{n}{d})\)。
我们想把这个恶心的\(d^2\)消掉。这样的话，我们不如设\(g=id^2\)。这样的话，\((f*g)(n)=\sum\limits_{d|n}\varphi(d)d^2(\dfrac{n}{d})^2=\sum\limits_{d|n}\varphi(d)n^2=n^3\)！！！
我们好像找对了。
再来看式子：
\(g(1)S(n)=\sum_{d=1}^n(f*g)(d)-\sum_{d=2}^ng(d)S(\dfrac{n}{d})\)
代入，得\(S(n)=\sum_{i=1}^ni^3-\sum_{d=2}^ni^2S(\dfrac{n}{d})\)
我们有公式：
\(\sum\limits_{i=1}^ni^2=\dfrac{n(n+1)(2n+1)}{6}\)
\(\sum\limits_{i=1}^ni^3=\dfrac{n^2(n+1)^2}{4}\)
这样我们就可以杜教筛了。
我们回到一开始的式子：\(\sum_{x=1}^n\varphi(x)x^2(\sum_{i=1}^{n/x}i)^2\)。后半部分是\(O(\sqrt{n})\)的整除分块，前半部分要用杜教筛，这样，看起来复杂度为\(O(n^{\frac{2}{3}}*\sqrt{n})=O(n^{\frac{7}{6}})\)。
Emm?这东西一看就不像能过的样子呀？
首先，因为杜教筛是基于记忆化搜索实现的，多次使用复杂度明显是低于\(O(n^{\frac{2}{3}})\)的。至于低多少，就要看出题人是否是用脚造数据了。
同时，因为整除分块的过程中，杜教筛要筛的\(n\)不会总是\(n\)级别的。总复杂度应该低于\(O(\sum_{i=1}^n(\dfrac{n}{i})^{\frac{2}{3}})\)
并且，因为杜教筛是基于分块思想的，我们一开始预处理的部分也可以大于\(n^{\frac{2}{3}}\)，这样单次杜教筛复杂度就会低于\(O(n^{\frac{2}{3}})\)。这样一番优化下来，这个算法复杂度应该是亚线性，甚至直接是\(O(n^{\frac{2}{3}})\)的（其实是我不会证）
或者换一种说法，它的复杂度是\(O(\text{能过})\)。
代码：
#pragma GCC optimize(3)
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=4000000;
int p,phi[N+10],inv6,pri[N+10],n;
inline int ksm(int x,int y){
register int res=1;
for(;y;x=(1ll*x*x)%p,y>>=1)if(y&1)res=(1ll*res*x)%p;
return res;
}
inline void init(){
inv6=ksm(6,p-2);
phi[1]=1;
for(register int i=2;i<=N;i++){
if(!pri[i])pri[++pri[0]]=i,phi[i]=i-1;
for(register int j=1;j<=pri[0]&&i*pri[j]<=N;j++){
pri[i*pri[j]]=true;
if(!(i%pri[j])){phi[i*pri[j]]=(1ll*phi[i]*pri[j])%p;break;}
phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
for(register int i=1;i<=N;i++)phi[i]=(phi[i-1]+1ll*phi[i]*i%p*i%p)%p;
}
inline int sqrsum(int x){
x%=p;
return 1ll*x*(x+1)%p*(2*x+1)%p*inv6%p;
}
inline int cubsum(int x){
x%=p;
return (1ll*((1ll*x*(x+1)/2)%p)*((1ll*x*(x+1)/2)%p))%p;
}
unordered_map<int,int>mp;
inline int djs(int x){
if(x<=N)return phi[x];
if(mp[x])return mp[x];
register int res=cubsum(x);
for(register int l=2,r;l<=x;l=r+1)r=x/(x/l),res=(res-(1ll*(sqrsum(r)-sqrsum(l-1)+p)%p*djs(x/l)%p)+p)%p;
return mp[x]=res;
}
inline int solve(){
register int res=0;
for(register int l=1,r;l<=n;l=r+1)r=n/(n/l),res=(1ll*(djs(r)-djs(l-1)+p)%p*cubsum(n/l)%p+res)%p;
return res;
}
signed main(){
scanf("%lld%lld",&p,&n),init();
printf("%lld\n",solve());
return 0;
}