题意：(sum_{i=1}^n sum_{j=1}^n sigma(i imes j), n leq 1e9)
题解：惭愧。
说到约数和，我们可能能想起约数个数和：(d(i imes j)=sum_{p|i}sum_{q|j}[(x,y)==1]) 。
我们考虑把每个因子一一映射。
如果 (i imes j) 的因子 (k) 中有一个因子 (p^c)，(i) 中有因子 (p^a)，(j) 中有因子 (p^b) 。我们规定：
- 如果 (cleq a)，那么在 (i) 中选择。
- 如果 (c>a)，那么我们把 (c) 减去 (a)，在 (j) 中选择 (p^{c-a})（在 (j) 中选择 (p^e) 表示的是 (p^{a+e})）
对于 (i imes j) 的因子 (k) 的其他因子同理。于是对于任何一个 (k) 有一个唯一的映射，且每一个选择对应着唯一的 (k) 。
通过如上过程，我们发现：对于 (i imes j) 的因子 (k=prod {p_i}^{c_i}) ，我们不可能同时在 (i) 和 (j) 中选择 (p_i) (优先在 (i) 中选择，如果不够就只在 (j) 中选择不够的指数），故 (x) 和 (y) 必须互质。等式得证。
From Siyuan@luogu
有一个式子： (sigma(i imes j)=sum_{x|i}sum_{y|j} frac{i}{x}y[(x,y)==1])
证明仍然类似上面考虑，如果 (i imes j) 的因子 (k) 中有一个因子 (p^c)，(i) 中有因子 (p^a)，(j) 中有因子 (p^b) 。我们规定：
- 如果 (cleq a)，那么在 (i) 中选择，这个时候 (x) 中包含 (p^{a-c}) ，所以此时 (frac{i}{x}) 贡献了 (p^c)。
- 如果 (c>a)，那么我们把 (c) 减去 (a)，在 (j) 中选择 (p^{c-a}) ，即此时 (y) 中包含 (p^{c-a}) 。
此时对于 (i imes j) 的因子 (k=prod {p_i}^{c_i}) 仍然有唯一的映射。
所以可以化式子了：
(=sum_{i=1}^n sum_{j=1}^n sum_{x|i} sum_{y|j} frac{i}{x}y [(x,y)==1])
(=sum_{i=1}^n sum_{j=1}^n sum_{x|i} sum_{y|j} frac{i}{x}y sum_{d|gcd(x,y)} mu(d))
$=sum_{i=1}^n sum_{j=1}^n sum_{d|gcd(i,j)} mu(d) sum_{x|frac{i}{d}} sum_{y|frac{j}{d}} frac{i}{x}y $
$=sum_{d=1}^{n}d mu(d) sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{n}{d}} sum_{x|i} sum_{y|j} frac{i}{x}y $
$=sum_{d=1}^{n}d mu(d) (sum_{i=1}^{frac{n}{d}} sum_{x|i} frac{i}{x}) (sum_{j=1}^{frac{n}{d}} sum_{y|j} y) $
$=sum_{d=1}^{n}d mu(d) (sum_{i=1}^{frac{n}{d}} sigma(i))^2 $
然后 (sum_{i=1}^{frac{n}{d}} sigma(i)) 注意到ta的上界显然是可以除法分块的，此时我们再在前头套个杜教筛即可。
最后就是 (sum_{i=1}^n sigma(i) = sum_{i=1}^n ilfloor frac{n}{i}
floor) ，表示 (i) 出现了几次。
#pragma GCC optimize(2)
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#define R register int
#define ll long long
using namespace std;
namespace Luitaryi {
const int N=1000000,M=1000000007;
int n,cnt;
int mu[N+10],p[N>>1];
ll s[N+10],c[N+10];
bool vis[N+10];
inline void pre(int n) {
mu[1]=1,s[1]=1;
for(R i=2;i<=n;++i) {
if(!vis[i]) p[++cnt]=i,mu[i]=-1,s[i]=i+1,c[i]=i+1;
for(R j=1,k;j<=cnt&&i*p[j]<=n;++j) {
k=i*p[j]; vis[k]=true;
if(i%p[j]==0) {
c[k]=c[i]*p[j]+1;
s[k]=s[i]/c[i]*c[k];
break;
} mu[k]=-mu[i];
c[k]=p[j]+1;
s[k]=s[i]*(p[j]+1);
}
}
for(R i=1;i<=n;++i) mu[i]=(mu[i-1]+mu[i]*i+M)%M;
for(R i=1;i<=n;++i) s[i]=(s[i-1]+s[i])%M;
}
int memmu[64000],SZ;
#define sqr(x) (1ll*(x)*(x+1)/2%M)
inline int calmu(int x) {
if(x<=N) return mu[x];
R& ret=x>SZ?memmu[SZ+n/x]:memmu[x];
if(ret!=memmu[0]) return ret;
R ans=1,lst=1,now;
for(R l=2,r;l<=x;l=r+1) {
r=x/(x/l),now=sqr(r);
ans=(ans-1ll*(now-lst)*calmu(x/l))%M;
lst=now;
}
return ret=(ans+M)%M;
}
inline int cals(int x) {
if(x<=N) return s[x];
R ret=0,lst=0,now;
for(R l=1,r;l<=x;l=r+1) {
r=x/(x/l);
now=sqr(r);
ret=(ret+1ll*(now-lst)*(x/l))%M;
lst=now;
} return (ret+M)%M;
}
inline void main() {
pre(N); scanf("%d",&n); SZ=sqrt(n);
memset(memmu,0xcf,sizeof memmu);
R ans=0,lst=0,now,tmp;
for(R l=1,r;l<=n;l=r+1) {
r=n/(n/l),now=calmu(r),tmp=cals(n/l);
tmp=1ll*tmp*tmp%M;
ans=(ans+1ll*(now-lst)*tmp)%M;
lst=now;
} printf("%d
",(ans+M)%M);
}
} signed main() {Luitaryi::main(); return 0;}
2019.12.23