luogu P1829 [国家集训队]Crash的数字表格 / JZPTAB

题面传送门
题目就是让我们求这个东西:\(\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{m}{lcm(i,j)}}\)
然后根据小学奥数就可以化成\(\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{m}{\frac{ij}{gcd(i,j)}}}\)
套路提前枚举gcd就可以去掉分式变成\(\sum\limits_{d=1}^{n}{d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}{\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}{ij[gcd(i,j)==1]}}}\)
设\(sum(n,m)=\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{m}{ij[gcd(i,j)==1]}}\)那么原式就变为\(\sum\limits_{d=1}^{n}{d\times sum(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\)这个显然可以数论分块求解。
然后考虑\(sum\)怎么求。
看到这个东西想到经典莫反，套路上去就有\(\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{m}{ij\sum\limits_{d|gcd(i,j)}{\mu(d)}}}\)
\(\sum\limits_{d=1}^{n}{\mu(d)\times d^2\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}{\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}{ij}}}\)
我们设\(g(n,m)=\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{m}{ij}}\)
容易发现这个可以\(O(1)\)计算，具体的这个式子就是\(g(n,m)=\frac{n(n+1)}{2}\times \frac{m(m+1)}{2}\)
然后就有\(sum(n,m)=\sum\limits_{d=1}^{n}{\mu(d)\times d^2\times g(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\)
容易发现这个东西也是可以整除分块算的。然后就可以愉快地\(O(n+m)\)预处理，\(O(n^{\frac{3}{4}})\)整除分块计算了。
然而这道题数据其实很水，第二个整除分块如果暴力计算的话复杂度调和级数的\(O(nlogn)\)也可以过。
code:

#include<cstdio>
#define min(a,b) ((a)<(b)?(a):(b))
#define mod 20101009
#define ll long long
using namespace std;
int n,m,k,x,y,z,flag[10000039],pr[10000039],mu[10000039],ph;
ll q[10000039],su[10000039],ans;
inline void swap(int &x,int &y){x^=y^=x^=y;}
inline ll g(ll x,ll y){return x*(x+1)/2%mod*(y*(y+1)/2%mod)%mod;}
inline ll sum(int x,int y){
	register int i,j;ll ans=0;
	/*for(i=1;i<=x;i=j+1){
		j=min(x/(x/i),y/(y/i));
		ans+=(su[j]-su[i-1])*g(x/i,y/i)%mod;
	}*/
	for(i=1;i<=x;i++)ans+=mu[i]*(ll)i*i*g(x/i,y/i)%mod;
	return (ans%mod+mod)%mod;
}
int main(){
	freopen("1.in","r",stdin);
	register int i,j;
	scanf("%d%d",&n,&m);mu[1]=1;(n>m)&&(swap(n,m),0);
	for(i=2;i<=n;i++){
		if(!flag[i]) pr[++ph]=i,mu[i]=-1;
		for(j=1;j<=ph&&i*pr[j]<=n;j++){
			flag[i*pr[j]]=1;
			if(i%pr[j]==0) break;
			else mu[i*pr[j]]=-mu[i];
		}
	}
	for(i=1;i<=n;i++) su[i]=su[i-1]+mu[i]*(ll)i*i,su[i]=(su[i]%mod+mod)%mod;
	for(i=1;i<=n;i++) q[i]=q[i-1]+i,q[i]%=mod;
	for(i=1;i<=n;i=j+1){
		j=min(n/(n/i),m/(m/i));
		ans+=(q[j]-q[i-1])*sum(n/i,m/i)%mod;
	}
	ans=(ans%mod+mod)%mod;printf("%d\n",ans);
}