[国家集训队]Crash的数字表格 / JZPTAB
这个题可以低于线性，然后也可以杜教筛到 \(O(n^{2 / 3})\) 这个样子。
首先暴力推：
\[\begin{aligned} & \sum_{i=1}^{n}\sum_{j=1}^{m} \operatorname{lcm}(i,j) \\ =& \sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij}{\gcd(i,j)} \\ =& \sum_{d=1}\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij[\gcd(i,j)=d]}{d} \\ =& \sum_{d=1}\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor m /d \rfloor}\dfrac{ijd^2}{d}[\gcd(i,j)=1] \\ =& \sum_{d=1}d\sum_{i=1}^{\lfloor n / d \rfloor}\sum_{j=1}^{\lfloor m / d \rfloor}ij\sum_{k\mid \gcd(i,j)}\mu(k) \\ =& \sum_{d=1}d\sum_{k=1}\mu(k)k^2\sum_{i=1}^{\lfloor n / kd\rfloor}i\sum_{j=1}^{\lfloor m / kd\rfloor}j \end{aligned} \]设 \(h(n)=\sum_{i=1}^{n}i={n(n+1)} / {2}\)，于是这个就是：
\[\begin{aligned} =& \sum_{d=1}d\sum_{k=1}\mu(k)k^2 h\left(\left\lfloor\dfrac{n}{kd}\right\rfloor\right) h\left(\left\lfloor \dfrac{m}{kd}\right\rfloor\right) \end{aligned} \]好了，现在我们设 \(h'(d)=\sum_{k=1}\mu(k)k^2h(\lfloor n / kd\rfloor)h(\lfloor m / kd\rfloor)\)。
\[\begin{aligned}=& \sum_{d=1}dh'(d)\end{aligned} \]这个显然可以数论分块 \(O(\sqrt{n})\) 做，然后 \(h'(d)\) 又是 \(O(\sqrt{n})\) 数论分块做。加上线性筛之类的东西，这个复杂度就是 \(O(n)\) 的了。
这个题已经能做出来了。然而线性还不够 nb，我们还能有低于线性的解法。
首先设 \(T=kd\)，则 \(k= T / d\)，于是：
\[\begin{aligned} =& \sum_{T=1}\sum_{d\mid T}d\mu \left(\dfrac{T}{d}\right) \left(\dfrac{T}{d}\right)^2 h\left(\left\lfloor\dfrac{n}{T}\right\rfloor\right) h\left(\left\lfloor\dfrac{m}{T}\right\rfloor\right) \end{aligned} \]我们来研究这个比较怪异的卷积，看一看能不能杜教筛一类，先设：
\[f=(\mathbf{id}^2\cdot \mu)*\mathbf{id} \]很容易能得到它的贝尔级数：
\[f _ p ( x ) = \dfrac { 1 - p ^ 2 x } { 1 - p x } \]然后我们肯定取 \(g=\mathbf{id}^2\)，然后卷上去：
\[(\mathbf{id}^2\cdot \mu)* \mathbf{id} * \mathbf{id^2}=\mathbf{id} \]那么 \(g\) 和 \(f*g\) 的前缀和都非常好求。\(g\) 的前缀和随便给一下：
\[S_g(n)=\sum_{i=1}^ng(i)=\dfrac{n(n+1)(2n+1)}{6} \]然后我们杜教筛直接套就完了。
时间复杂度 \(O(n^{2 / 3})\)。
杜教筛不熟的我可以写个式子：
\[g(1)S_f(n)=\sum_{i=1}^n(f \ast g)(i)-\sum_{y=2}^{n}g(y)S_f\left (\left\lfloor\dfrac{n}{y}\right\rfloor\right) \]然后这里的记忆化讲实话确实难搞，所以我直接躺平，用 Hash 随便取个质数模一下就完了。然后这里我取的模数是 \(786433\)。当然用 map 也是可以的。
因为 \(2\) 和 \(6\) 都是与 \(20101009\) 互质的，所以直接求个逆元。
然而因为我预处理直接躺平，所以最后跑得不如一开始提出的线性算法。
（现在把博客内容改了，但是我还是不想重写杜教筛的代码，你们看得懂就好了。）
代码（杜教筛 \(O(n^{2 / 3 })\)）：
#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
const ll N=1e7,M=1e6,mo=20101009,inv2=10050505,inv6=16750841,hmo=786433;
ll n,m,ans,cnt;
ll prime[M+5],mu[M+5],Sf1[M+5],Sf2[M+5];
bool f[M+5],vis[M+5];
inline void Init() {
mu[1]=1;f[1]=1;
for(ll i=2;i<=M;i++) {
if(!f[i]) {prime[++cnt]=i;mu[i]=-1;}
for(ll j=1;j<=cnt&&i*prime[j]<=M;j++) {
f[i*prime[j]]=1;
if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
for(ll i=1;i<=M;i++) {
for(ll j=i;j<=M;j+=i) {
Sf1[j]=(Sf1[j]+(mu[i]*i*i+mo)%mo*(j/i)%mo)%mo;
}
}
for(ll i=1;i<=M;i++) Sf1[i]=Sf1[i-1]+Sf1[i];
}
inline ll H(ll x) {return x%hmo;}
inline ll S1(ll x) {return x*(x+1)%mo*inv2%mo;}
inline ll S2(ll x) {return x*(x+1)%mo*(2*x+1)%mo*inv6%mo;}
inline ll Sf(ll x) {
if(x<=M) return Sf1[x];ll tmp=H(x);
if(vis[tmp]) return Sf2[tmp];
vis[tmp]=1;Sf2[tmp]=S1(x);
for(ll i=2,j;i<=x;i=j+1) {
j=x/(x/i);
Sf2[tmp]=(Sf2[tmp]-(S2(j)-S2(i-1)+mo)%mo*Sf(x/i)%mo+mo)%mo;
}
return Sf2[tmp];
}
inline ll read() {
ll ret=0,f=1;char ch=getchar();
while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
return ret*f;
}
inline void write(ll x) {
static char buf[22];static ll len=-1;
if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
else {putchar(45);do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
while(len>=0) putchar(buf[len--]);
}
int main() {
n=read();m=read();if(n<m) swap(n,m);
Init();
for(ll i=1,j;i<=m;i=j+1) {
j=min(n/(n/i),m/(m/i));
ans=(ans+(Sf(j)-Sf(i-1)+mo)*S1(n/i)%mo*S1(m/i)%mo)%mo;
}
write(ans);
return 0;
}
代码（\(O(n)\)）：
#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
const ll N=1e7,mo=20101009;
ll n,m,ans,cnt;
ll prime[N+5],mu[N+5],s[N+5];
bool f[N+5];
inline void Init() {
mu[1]=1;f[1]=1;
for(ll i=2;i<=n;i++) {
if(!f[i]) {prime[++cnt]=i;mu[i]=-1;}
for(ll j=1;j<=cnt&&i*prime[j]<=n;j++) {
f[i*prime[j]]=1;
if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
for(ll i=1;i<=n;i++) {
s[i]=(s[i-1]+mu[i]*i%mo*i%mo)%mo;
}
}
inline ll S(ll x) {return (x*(x+1)/2)%mo;}
inline ll F(ll x) {
ll res=0;
for(ll i=1,j;i<=m/x;i=j+1) {
j=min((n/x)/(n/x/i),(m/x)/(m/x/i));
res=(res+(s[j]-s[i-1]+mo)*S(n/x/i)%mo*S(m/x/i)%mo)%mo;
}
return res;
}
inline ll read() {
ll ret=0,f=1;char ch=getchar();
while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
return ret*f;
}
inline void write(ll x) {
static char buf[22];static ll len=-1;
if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
else {putchar(45);do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
while(len>=0) putchar(buf[len--]);
}
inline void writeln(ll x) {write(x);putchar(10);}
int main() {
n=read();m=read();
if(n<m) swap(n,m);Init();
for(ll i=1,j;i<=m;i=j+1) {
j=min(n/(n/i),m/(m/i));
ans=(ans+(S(j)-S(i-1)+mo)*F(i)%mo)%mo;
}
write(ans);
return 0;
}