[BZOJ2154]Crash的数字表格
题目大意：
给定$n,m(n,m\leq10^7)$，求$\displaystyle\sum_{x=1}^n\sum_{y=1}^m{\rm lcm}(x,y)$。
思路：
令$d=\gcd(x,y),x=ad,y=bd$，则：
$$
\begin{align*}
原式&=\sum_{d=1}^{\min(n,m)}d\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}ab[\gcd(ab)=1]\\
&=\sum_{d=1}^{\min(n,m)}d\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}ab\sum_{p|\gcd(a,b)}\mu(p)
\end{align*}
$$
令$p=\gcd(a,b),a=jp,b=kp$，则：
$$
\begin{align*}
原式&=\sum_{d=1}^{\min(n,m)}d\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)p^2\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{dp}\rfloor}jk\\
&=\sum_{d=1}^{\min(n,m)}d\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)p^2\frac{(\lfloor\frac{n}{dp}\rfloor+1)\lfloor\frac{n}{dp}\rfloor}{2}\cdot\frac{(\lfloor\frac{m}{dp}\rfloor+1)\lfloor\frac{m}{dp}\rfloor}{2}\\
\end{align*}
$$
预处理$\mu(p)p^2$的前缀和，暴力枚举$d$。对于$\lfloor\frac{n}{dp}\rfloor$和$\lfloor\frac{m}{dp}\rfloor$分别相等的$p$可以分块计算。
1 #include<cstdio> 2 #include<cctype> 3 #include<algorithm> 4 typedef long long int64; 5 inline int getint() { 6 register char ch; 7 while(!isdigit(ch=getchar())); 8 register int x=ch^'0'; 9 while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0'); 10 return x; 11 } 12 const int N=10000001,M=664580,mod=20101009; 13 bool vis[N]; 14 int mu[N],p[M],sum[N]; 15 inline void sieve(const int lim) { 16 mu[1]=1; 17 for(register int i=2;i<=lim;i++) { 18 if(!vis[i]) { 19 p[++p[0]]=i; 20 mu[i]=-1; 21 } 22 for(register int j=1;j<=p[0]&&i*p[j]<=lim;j++) { 23 vis[i*p[j]]=true; 24 if(i%p[j]==0) { 25 mu[i*p[j]]=0; 26 break; 27 } 28 mu[i*p[j]]=-mu[i]; 29 } 30 } 31 for(register int i=1;i<=lim;i++) { 32 sum[i]=(sum[i-1]+(int64)mu[i]*i%mod*i%mod)%mod; 33 } 34 } 35 int main() { 36 const int n=getint(),m=getint(),lim=std::min(n,m); 37 sieve(lim); 38 int ans=0; 39 for(register int d=1;d<=lim;d++) { 40 int tmp=0; 41 const int x=n/d,y=m/d,lim=std::min(x,y); 42 for(register int i=1,j;i<=lim;i=j+1) { 43 j=std::min(x/(x/i),y/(y/i)); 44 (tmp+=(sum[j]-sum[i-1]+mod)%mod*((int64)(x/i+1)*(x/i)/2%mod)%mod*((int64)(y/i+1)*(y/i)/2%mod)%mod)%=mod; 45 } 46 (ans+=((int64)tmp*d)%mod)%=mod; 47 } 48 printf("%d\n",ans); 49 return 0; 50 }
