求 $\sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)$
考虑欧拉反演： $\sum_{d|n}\varphi(d)=n$
$\Rightarrow \sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{d|gcd(i,j)}\varphi(d)$
$\Rightarrow \sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{d|i,d|j}\varphi(d)$
$\Rightarrow \sum_{d=1}^{n}\varphi(d)\sum_{d|i}\sum_{d|j}ij$
$\Rightarrow\sum_{d=1}^{n} \varphi(d)d^2\sum_{i=1}^{\frac{n}{d}}i\sum_{j=1}^{\frac{n}{d}}j$
对于 $\sum_{i=1}^{\frac{n}{d}}i\sum_{j=1}^{\frac{n}{d}}j$，直接 $O(1)$求
令 $calc(n,m)=\frac{n(n+1)}{2}\times \frac{m(m+1)}{2}$
将 $\frac{n}{d}$ 带入即可.
原式可化为 $\sum_{d=1}^{n} \varphi(d)d^2\times calc(\frac{n}{d},\frac{n}{d})$
这个复杂度是 $O(\sqrt n)$ 的.
然而，有一个问题：我们正常只能求出 $[1,1e7]$ 的欧拉函数值.
于是，要用杜教筛优化一下.
令 $f(x)=\varphi(x)x^2$
搬出杜教是公式： $S(n)=\frac{\sum_{i=1}^{n}(f*g)(i)-\sum_{i=2}^{n}g(i)S(\frac{n}{i})}{g(1)}$
$f(x)$ 是固定的，现在只需选一个合适的 $g(x)$，使得可以快速算出 $\sum_{i=1}^{n}(f*g)(i)$
$(f*g)(i)=\sum_{d|i}f(d)g(\frac{i}{d})$
$\Rightarrow \sum_{d|i}\varphi(d)d^2g(\frac{i}{d})$
$d^2$ 有点讨厌，要是能消掉就好了.
令 $g(x)=x^2$
$\Rightarrow \sum_{d|i}\varphi(d)d^2(\frac{i}{d})^2$
$\Rightarrow \sum_{d|i}\varphi(d)i^2$
$\Rightarrow i^2\times\sum_{d|i}\varphi(d)$
$\Rightarrow i^3$
将 $(f*g)(i)$ 的结果带入杜教筛公式中.
即当 $f(x)=\varphi(x)x^2$, $g(x)=x^2$ 时，
$S(n)=\sum_{i=1}^{n}i^3-\sum_{i=2}^{n}i^2\times S(\frac{n}{i})$
搬出高中数学公式：
$\sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}{6}$
$\sum_{i=1}^{n}i^3=\frac{n^2(n+1)^2}{4}$
回到原式 $\sum_{d=1}^{n} \varphi(d)d^2\times calc(\frac{n}{d},\frac{n}{d})$
直接用杜教筛算 $\varphi(d)d^2$ 的前缀和并用整除分块算出结果即可.
#include<bits/stdc++.h>
#define maxn 10200006
#define ll long long
#define M 10000007
using namespace std;
int cnt;
ll sumv[maxn], rev4, rev6, mod, rev2;
bool vis[maxn];
ll phi[maxn], prime[maxn];
map<ll,ll>ansphi;
void setIO(string s)
{
string in=s+".in";
freopen(in.c_str(),"r",stdin);
}
ll qpow(ll base, ll k)
{
ll tmp=1;
while(k)
{
if(k&1) tmp=tmp*base%mod;
base=base*base%mod;
k>>=1;
}
return tmp;
}
void init()
{
int i,j;
rev4=qpow(4ll, mod-2), rev6=qpow(6ll, mod-2), rev2=qpow(2ll, mod-2);
phi[1]=1;
for(i=2;i<=M;++i)
{
if(!vis[i]) prime[++cnt]=i, phi[i]=i-1;
for(j=1;j<=cnt&&1ll*i*prime[j]<=M;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
for(i=1;i<=M;++i) sumv[i]=(sumv[i-1]+(1ll*phi[i]*i%mod*i%mod))%mod;
}
// 平方
ll cal1(ll i)
{
i%=mod;
ll re=i%mod;
re=re*(i+1)%mod;
re=re*(i+i+1)%mod;
re=(re*rev6)%mod;
return re;
}
// 立方
ll cal2(ll i)
{
i%=mod;
ll re=i%mod;
re=(re*i)%mod;
re=(re*(i+1))%mod;
re=re*(i+1)%mod;
re=(re*rev4)%mod;
return re;
}
ll get(ll n)
{
if(n<=M) return sumv[n];
if(ansphi[n]) return ansphi[n];
ll i,j,re=cal2(n),tmp;
for(i=2;i<=n;i=j+1)
{
j=n/(n/i);
tmp=(cal1(j)-cal1(i-1)+mod)%mod;
tmp=(tmp*get(n/i))%mod;
re=(re-tmp+mod)%mod;
}
return ansphi[n]=re;
}
ll calc(ll n)
{
n%=mod;
return (((n*(n+1))%mod)*(rev2%mod))%mod ;
}
int main()
{
// setIO("input");
ll n,i,j,re=0,tmp=0;
scanf("%lld%lld",&mod,&n);
init();
for(i=1;i<=n;i=j+1)
{
j=n/(n/i);
tmp=(calc(n/i)*calc(n/i)%mod*(get(j)-get(i-1)+mod)%mod)%mod;
re=(re+tmp+mod)%mod;
}
printf("%lld\n",re);
return 0;
}
