BZOJ 2154.Crash 的数字表格

显然题目就是要我们求 i=1nj=1mlcm(i,j)\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \text{lcm}(i,j)
根据小学奥数,就是 i=1nj=1mijgcd(i,j)\sum\limits_{i = 1}^n \sum\limits_{j = 1}^m \dfrac{ij}{\gcd(i,j)}

TeX parse error: \begin{align*} ended with \end{aligned}

敏锐地察觉到 [n=1][n = 1] 这个东西可以替换成 ϵ(n)=dnμ(d)\epsilon(n) = \sum\limits_{d | n} \mu(d)
那么有

i=1min(n,m)ixiniyimi[gcd(x,y)=1]xy=i=1min(n,m)ixiniyimixydgcd(x,y)μ(d) \begin{aligned} & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} [\gcd(x,y) = 1]xy \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} xy \sum\limits_{d | \gcd(x,y)} \mu(d) \end{aligned}

根据套路,接着把 dd 提出来。

i=1min(n,m)ixiniyimixydgcd(x,y)μ(d)=i=1min(n,m)id=1min(ni,mi)μ(d)d2xdniydmixy=i=1min(n,m)id=1min(ni,mi)μ(d)d2sum(ni,mi) \begin{aligned} & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{xi \le \frac n i} \sum\limits_{yi \le \frac m i} xy \sum\limits_{d | \gcd(x,y)} \mu(d) \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{d = 1}^{\min(\lfloor \frac n i \rfloor,\lfloor \frac m i \rfloor)} \mu(d) d^2 \sum\limits_{xd \le \frac n i} \sum\limits_{yd \le \frac m i} xy \\ = & \sum\limits_{i = 1}^{\min(n,m)} i \sum\limits_{d = 1}^{\min(\lfloor \frac n i \rfloor,\lfloor \frac m i \rfloor)} \mu(d) d^2 \text{sum}(\lfloor \dfrac n i \rfloor,\lfloor \dfrac m i \rfloor) \end{aligned}

其中 sum(x,y)=i=1xj=1yij=x(x+1)2y(y+1)2\text{sum}(x,y) = \sum\limits_{i = 1}^x \sum\limits_{j = 1}^y ij = \dfrac{x(x + 1)} 2 \dfrac{y(y + 1)} 2
那么预处理 μ(i)i2\mu(i)i^2 的前缀和,做两次数论分块即可。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include <cstdio>
#include <algorithm>
using namespace std;
const int N = 1e7;
const long long mod = 20101009;
int n,m;
int cnt,prime[N + 5],vis[N + 5],sum[N + 5];
long long ans;
long long get_sum(int n,int m)
{
long long ret = 0;
if(n > m)
swap(n,m);
for(register int l = 1,r;l <= n;l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ret = (ret + (long long)(sum[r] - sum[l - 1] + mod) * ((long long)(n / l) * (n / l + 1) / 2 % mod) % mod * ((long long)(m / l) * (m / l + 1) / 2 % mod) % mod);
}
return ret;
}
int main()
{
sum[1] = 1;
for(register int i = 2;i <= N;++i)
{
if(!vis[i])
prime[++cnt] = i,sum[i] = -1;
for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
sum[i * prime[j]] = -sum[i];
}
}
for(register int i = 1;i <= N;++i)
sum[i] = (sum[i - 1] + (long long)i * i % mod * (sum[i] + mod)) % mod;
scanf("%d%d",&n,&m);
if(n > m)
swap(n,m);
for(register int l = 1,r;l <= n;l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ans = (ans + (long long)(l + r) * (r - l + 1) / 2 % mod * get_sum(n / l,m / l) % mod) % mod;
}
printf("%lld\n",ans);
}