[BZOJ2154]Crash的数字表格

推公式（假设$n\leq m$）：$\begin{align*}\sum\limits_{i=1}^n\sum\limits_{j=1}^m[i,j]&=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\dfrac{ij}{(i,j)}\end{align*}$

枚举$d=(i,j)$，将上式化为$\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac md\right\rfloor}[(i,j)=1]ij$

设$f(n,m,k)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[(i,j)=k]ij$，则$F(n,m,k)=\sum\limits_{\substack{k|d\\d\leq n}}f(n,m,d)=k^2\dfrac{\left\lfloor\dfrac nk\right\rfloor\left(\left\lfloor\dfrac nk\right\rfloor+1\right)}2\dfrac{\left\lfloor\dfrac mk\right\rfloor\left(\left\lfloor\dfrac mk\right\rfloor+1\right)}2$

所以$f(n,m,k)=\sum\limits_{\substack{k|d\\d\leq n}}\mu\left(\dfrac dk\right)F(n,m,d)$

那么最初的式子可以化为$\sum\limits_{d=1}^ndf\left(\left\lfloor\frac nd\right\rfloor,\left\lfloor\frac md\right\rfloor,1\right)$，只要能快速求$f$就可以算出本题的答案

而$f(n,m,1)=\sum\limits_{d=1}^n\mu(d)F(n,m,d)=\sum\limits_{d=1}^nd^2\mu(d)\dfrac{\left\lfloor\dfrac nd\right\rfloor\left(\left\lfloor\dfrac nd\right\rfloor+1\right)}2\dfrac{\left\lfloor\dfrac md\right\rfloor\left(\left\lfloor\dfrac md\right\rfloor+1\right)}2$，需要预处理$d^2\mu(d)$的前缀和

于是就做完了

感谢yww指出我的sb错误==

#include<stdio.h>
#define mod 20101009
#define ll long long
int mu[10000010],pr[10000010];
bool np[10000010];
int ad(int a,int b){return(a+b)%mod;}
int mul(int a,int b){return a*(ll)b%mod;}
void sieve(int n){
int i,j,m=0;
np[1]=1;
mu[1]=1;
for(i=2;i<=n;i++){
if(!np[i]){
m++;
pr[m]=i;
mu[i]=-1;
}
for(j=1;j<=m;j++){
if(pr[j]*(ll)i>n)break;
np[i*pr[j]]=1;
if(i%pr[j]==0)break;
mu[i*pr[j]]=-mu[i];
}
}
for(i=1;i<=n;i++)mu[i]=mul(mu[i],mul(i,i));
for(i=2;i<=n;i++)mu[i]=ad(mu[i],mu[i-1]);
}
int min(int a,int b){return a<b?a:b;}
void swap(int&a,int&b){a^=b^=a^=b;}
int S(int n){return mul(mul(n,n+1),10050505);}
int f(int n,int m){
int i,nex,s=0;
for(i=1;i<=n;i=nex+1){
nex=min(n/(n/i),m/(m/i));
s=ad(s,mul((mu[nex]-mu[i-1])%mod,mul(S(n/i),S(m/i))));
}
return s;
}
int mob(int n,int m){
int i,nex,s=0;
for(i=1;i<=n;i=nex+1){
nex=min(n/(n/i),m/(m/i));
s=ad(s,mul((S(nex)-S(i-1))%mod,f(n/i,m/i)));
}
return s;
}
int main(){
int n,m;
scanf("%d%d",&n,&m);
if(n>m)swap(n,m);
sieve(n);
printf("%d",(mob(n,m)+mod)%mod);
}