BZOJ-2154 Crash的数字表格
BZOJ-2154 Crash的数字表格
解法
请注意！下文中默认 \(n\leq m\) ！！！
一步一步推：
\[\begin{aligned}
f(n,m)=&\sum\limits_{i=1}^n\sum\limits_{j=1}^m \operatorname{lcm}(i,j)
\\f(n,m)=&\sum\limits_{i=1}^n\sum\limits_{j=1}^m \frac{i\cdot j}{\gcd(i,j)}
\\f(n,m)=&\sum\limits_{i=1}^n\sum\limits_{j=1}^m \gcd(i,j)\cdot\frac{i}{\gcd(i,j)}\cdot\frac{j}{\gcd(i,j)}
\end{aligned}\]
我们调整一下枚举顺序可得（\(d\) 表示 \(\gcd\)）：
\[\begin{aligned}
f(n,m)=&\sum\limits_{d=1}^n d\cdot \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} [\gcd(i,j)=1]i\cdot j
\end{aligned}\]
我们先不管前面的部分，设 \(g(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} [\gcd(i,j)=1]i\cdot j\) 我们看到 \([\gcd(i,j)=1]\) 自然可以想到用 \(\varepsilon\) 表示，再用 \(\mu\ast 1\) 来表示：
\[\begin{aligned}
g(n,m)=&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \varepsilon(\gcd(i,j))\cdot i\cdot j
\\g(n,m)=&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} (\mu\ast 1)(\gcd(i,j))\cdot i\cdot j
\\g(n,m)=&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d\mid \gcd(i,j)} \mu(d)\cdot i\cdot j
\end{aligned}\]
故技重施，我们还是先枚举 \(d\)：
\[\begin{aligned}
g(n,m)=&\sum\limits_{d=1}^n\mu(d)\cdot\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor} (d\cdot i)\cdot(d\cdot j)
\\g(n,m)=&\sum\limits_{d=1}^n\mu(d)\cdot d^2\cdot\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor} i\cdot j
\end{aligned}\]
我们再次剥离式子，设 \(h(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} i\cdot j\)
\[\begin{aligned}
h(n,m)=&(\sum\limits_{i=1}^{n} i)\cdot(\sum\limits_{j=1}^{m} j)
\\h(n,m)=&\frac{n*(n+1)}{2}\cdot\frac{m*(m+1)}{2}
\end{aligned}\]
好的我们的式子已经推到了可以用 \(O(1)\) 时间可以求解的式子了，我们再来一步一步带回去。
\[g(n,m)=\sum\limits_{d=1}^n\mu(d)\cdot d^2\cdot h(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)
\]
可以用数论分块解决，\(\mu(d)\cdot d^2\) 可以用前缀和来维护。
\[f(n,m)=\sum\limits_{d=1}^n d\cdot g(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)
\]
真巧，还是数论分块。总时间复杂度 \(O(n+m)\)
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!
#include<bits/stdc++.h>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=1e7+10,mod=20101009;
int n,m,mu[maxn],p[maxn/3],cnt,s[maxn];
bool is[maxn];
void sieve(int n) {
mu[1]=1;
for(int i=2;i<=n;i++) {
if(!is[i]) {
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;j++) {
if(i*p[j]>n) {
break;
}
is[i*p[j]]=1;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[p[j]*i]=-mu[i];
}
}
for(int i=1;i<=n;i++) {
s[i]=(s[i-1]+(mu[i]+mod)%mod*i%mod*i%mod)%mod;
}
}
int h(int n,int m) {
return (n*(n+1)/2)%mod*(m*(m+1)/2%mod)%mod;
}
int g(int n,int m) {
int ret=0,j;
for(int i=1;i<=n;i=j+1) {
j=min(n/(n/i),m/(m/i));
ret+=(s[j]-s[i-1]+mod)%mod*h(n/i,m/i)%mod;
ret%=mod;
}
return ret;
}
int f(int n,int m) {
int ret=0,j;
for(int i=1;i<=n;i=j+1) {
j=min(n/(n/i),m/(m/i));
ret+=(j-i+1)*(i+j)/2%mod*g(n/i,m/i)%mod;
ret%=mod;
}
return ret;
}
signed main() {
cin>>n>>m;
if(n>m) {
swap(n,m);//这一步至关重要！！！
}
sieve(n);
printf("%lld\n",f(n,m));
return 0;
}
