\(\text{Description}\)
\(\text{Solution}\)
现在真的干啥啥不行，吃饭第一名。
方法壹
其实就是 7.3.3.用法叁，可以得到：
\[\sum_{T=1}^n\text{sum}\left(\frac{n}{T}\right)\times T^2\sum_{d|T}\mu\left(\frac{T}{d}\right)\times d
\]
如果是普通莫比乌斯反演题，到这里就可以预处理 \(+\) 分块做了。
但这题不行。
观察 \(\sum_{d|T}\mu\left(\frac{T}{d}\right)\times d\) 的结构，其实就是 \(\mu*\text{id}=\varphi\)。
那我们只需要预处理 \(i^2\times \varphi(i)\) 的前缀和，显然这是个积性函数，试试杜教筛。
令 \(f(i)=i^2\times \varphi(i)\)，为了卷积时抵消 \(i^2\) 项，令 \(g(i)=i^2\)。
则：
\[h(n)=(f*g)(n)=\sum_{d|n}f(d)g\left(\frac{n}{d}\right)
\]
\[=\sum_{d|n}d^2\times \varphi(d)\times \left(\frac{n}{d}\right)^2
\]
\[=n^2\sum_{d|n}\varphi(d)1\left(\frac{n}{d}\right)
\]
由于 \(\varphi*1=\text{id}\)，则：
\[=n^2\times n=n^3
\]
带入杜教筛的套路式，做 \(\text{sum}\) 的分块套上杜教筛套路式的分块即可。
方法贰
当然可以
\[\sum_{i=1}^n\sum_{j=1}^ni\times j\times \text{id}(\gcd(i,j))
\]
\[=\sum_{i=1}^n\sum_{j=1}^ni\times j\times \sum_{d|i,d|j}\varphi(d)
\]
由于 \(\gcd\) 是满足同时为 \(i,j\) 的因数的最大的数，所以 \(d\) 同时为 \(i,j\) 的因数就一定是 \(\gcd\) 的因数，就可以巧妙地去掉 \(\gcd\)。
有：
\[=\sum_{d=1}^n \varphi(d)\sum_{d|i}\sum_{d|j}i\times j
\]
\[=\sum_{d=1}^n\varphi(d)\times d^2\left(\sum_{i=1}^{n/d}i\right)^2
\]
然后和 方法壹 一样分块再杜教筛即可。
\(\text{Code}\)
#include <cstdio>
#define rep(i,_l,_r) for(register signed i=(_l),_end=(_r);i<=_end;++i)
#define fep(i,_l,_r) for(register signed i=(_l),_end=(_r);i>=_end;--i)
#define erep(i,u) for(signed i=head[u],v=to[i];i;i=nxt[i],v=to[i])
#define efep(i,u) for(signed i=Head[u],v=to[i];i;i=nxt[i],v=to[i])
#define print(x,y) write(x),putchar(y)
template <class T> inline T read(const T sample) {
T x=0; int f=1; char s;
while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
return x*f;
}
template <class T> inline void write(const T x) {
if(x<0) return (void) (putchar('-'),write(-x));
if(x>9) write(x/10);
putchar(x%10^48);
}
template <class T> inline T Max(const T x,const T y) {if(x>y) return x; return y;}
template <class T> inline T Min(const T x,const T y) {if(x<y) return x; return y;}
template <class T> inline T fab(const T x) {return x>0?x:-x;}
template <class T> inline T gcd(const T x,const T y) {return y?gcd(y,x%y):x;}
template <class T> inline T lcm(const T x,const T y) {return x/gcd(x,y)*y;}
template <class T> inline T Swap(T &x,T &y) {x^=y^=x^=y;}
#include <map>
using namespace std;
typedef long long ll;
const int maxn=4e6;
int mod,inv6,p[maxn>>2],pc,phi[maxn+5];
bool is[maxn+5];
ll n;
map <ll,int> mp;
int qkpow(int x,int y) {
int r=1;
while(y) {
if(y&1) r=1ll*r*x%mod;
x=1ll*x*x%mod; y>>=1;
}
return r;
}
int calc2(ll x) {
x=x%mod;
return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
int calc3(ll x) {
x=x%mod;
return x*(x+1)/2%mod*(x*(x+1)/2%mod)%mod;
}
void init() {
phi[1]=1;
rep(i,2,maxn) {
if(!is[i]) p[++pc]=i,phi[i]=i-1;
rep(j,1,pc) {
if(i*p[j]>maxn) break;
is[i*p[j]]=1;
if(i%p[j]==0) {phi[i*p[j]]=phi[i]*p[j]; break;}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
rep(i,1,maxn) phi[i]=(0ll+phi[i-1]+1ll*i*i%mod*phi[i]%mod)%mod;
}
int S(ll x) {
if(x<=maxn) return phi[x]%mod;
if(mp.count(x)) return mp[x];
int ret=0; ll r;
ret=calc3(x);
for(ll l=2;l<=x;l=r+1) {
r=(x/(x/l));
ret=(0ll+ret-1ll*(0ll+calc2(r)-calc2(l-1)+mod)%mod*S(x/l)%mod+mod)%mod;
}
mp.insert(make_pair(x,ret));
return ret;
}
int Work() {
ll r; int ret=0;
for(ll l=1;l<=n;l=r+1) {
r=(n/(n/l));
ret=(ret+1ll*calc3(n/l)*(0ll+S(r)-S(l-1)+mod)%mod)%mod;
}
return ret;
}
int main() {
mod=read(9),n=read(9ll); inv6=qkpow(6,mod-2);
init();
print(Work(),'\n');
return 0;
}
