洛谷P1829 [国家集训队] Crash的数字表格 / JZPTAB 题解

题意：

给定 \(n,m\) 求 \[ \sum_{i=1}^n \sum_{j=1}^m \operatorname{lcm}(i, j) \] 模 \(20101009\) 的结果。

数据范围：\(1\le n,m \le 10^7\) 。

来讲一个杜教筛的解法吧。

首先先规定一下 \(n \le m\) ，然后开始推式子 \[ \begin{aligned} \mathrm{Ans} &= \sum_{i=1}^n \sum_{j=1}^m \frac{i \cdot j}{\gcd(i, j)} \\[6pt]&= \sum_{i=1}^n \sum_{j=1}^m \sum_{d|i, d| j}\left[\gcd\left(\frac{i}{d}, \frac{j}{d}\right)=1\right]\cdot \frac{i \cdot j}{d} \\[6pt]&= \sum_{d=1}^n d\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\operatorname{gcd}(i, j)=1] \cdot i \cdot j \end{aligned} \] 不妨设后面一堆东西为 \(g\left(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor\right)\) ，那么 \[ \begin{aligned} g(n,m) &= \sum_{i=1}^n \sum_{j=1}^m[\operatorname{gcd}(i, j)=1] \cdot i \cdot j \\[6pt] &= \sum_{d=1}^n \mu(d) \cdot d^2 \cdot \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor} i \cdot j \\[6pt] &= \sum_{d=1}^n \mu(d) \cdot d^2 \cdot S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)S\left(\left\lfloor\frac{m}{d}\right\rfloor\right) && S(n) = \frac{n(n+1)}{2} \end{aligned} \] 那么代回原式可得 \[ \begin{aligned} \mathrm{Ans} &=\sum_{d=1}^n d \cdot g\left(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor\right) \\[6pt]&=\sum_{i=1}^n i\sum_{d=1}^{\left\lfloor\frac{n}{i}\right\rfloor} \mu(d) \cdot d^2 \cdot S\left(\left\lfloor\frac{n}{di}\right\rfloor\right)S\left(\left\lfloor\frac{m}{di}\right\rfloor\right) \\[6pt]&=\sum_{T=1}^nS\left(\left\lfloor\frac{n}{T}\right\rfloor\right)S\left(\left\lfloor\frac{m}{T}\right\rfloor\right) \sum_{d \mid T}d^2\mu(d)\frac{T}{d} && \textbf{Let } T = di. \end{aligned} \] 设 \(f=(\mu\, \mathrm{id}^2) *\mathrm{id}\) ，也就是 \(S\) 后面那个东西。

考虑构造一个 \(g\) 以快速计算 \(f\) 的前缀和。这里直接说了，就是 \(g=\mathrm{id}^2\) ，那么 \[ \begin{aligned} (f*g)(n) &= (\mu\, \mathrm{id}^2) *\mathrm{id} * \mathrm{id}^2 \\[6pt]&=(\mu\, \mathrm{id}^2) *\mathrm{id}^2 * \mathrm{id} \\[6pt]&=\left(\sum_{d \mid n}\mu(d)d^2\cdot \frac{n^2}{d^2}\right)*\mathrm{id} \\[6pt]&=\epsilon * \mathrm{id} \\[6pt]&=\sum_{d \mid n}\epsilon(d)\frac{n}{d} \\[6pt]&= n = \mathrm{id}(n) \end{aligned} \] 这里第一行利用了狄利克雷卷积的交换律。

其实到这里我们就可以在 \(\mathcal{O}(n^{\frac{3}{4}})\) 的复杂度内求解出答案了

但是还能更快，前提是我们能够提前筛出前 \(\mathcal{O}(n^{\frac{2}{3}})\) 个答案。

由于 \(f(T) = \sum_{d \mid T}d\mu(d)\) 是一个积性函数，不妨考察为 \(T\) 添加质因子 \(p\) 得到 \(Tp\) 时

若 \(T\) 有质因子 \(p\) ，乘以 \(p\) 不会产生任何贡献，故此时 \(f(Tp) = f(T)\) 。
若 \(T\) 没有质因子 \(p\) ，乘以 \(p\) 会使 \(\mu(Tp)=-\mu(T)\) ，故此时 \(f(Tp)=f(T) -pf(T) = f(p)f(T)\)

考虑线性筛筛出前 \(\mathcal{O}(n^{\frac{2}{3}})\) 个答案即可。总时间复杂度为 \(\mathcal{O}(n^{\frac{2}{3}})\)

代码：（非常快！）

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
const int mod = 20101009, inv6 = 16750841;
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define rep(i, a, b) for(int i = (a), i##END = (b); i <= i##END; i++)
#define Rep(i, a, b) for(int i = (a), i##END = (b); i >= i##END; i--)
#define N ((int)(5e6 + 15))

char ck[N]; int M, pcnt, p[N], f[N], sum[N];
unordered_map<int, int> vis;
void add(int &x, int y) { (x += y) >= mod ? x -= mod : 0; }
int qpow(int a, int b)
{
    int r = 1;
    for(a %= mod; b > 0; b >>= 1, a = a * a % mod)
        if(b & 1) r = r * a % mod;
    return r;
}
void init(const int _ = N - 5)
{
    f[1] = 1;
    rep(i, 2, _)
    {
        if(!ck[i]) { p[++pcnt] = i; f[i] = (1 + mod - i); }
        for(int j = 1; j <= pcnt && i * p[j] <= _; j++)
        {
            const int pos = i * p[j]; ck[pos] = true;
            if(i % p[j] == 0) { f[pos] = f[i]; break; }
            else { f[pos] = f[i] * f[p[j]] % mod; }
        }
    }
    rep(i, 1, _) add(sum[i] = sum[i - 1], f[i] * i % mod);
}
int fun1(int x) { x %= mod; return (x * (x + 1) / 2) % mod; }
int fun2(int x) { x %= mod; return x * (x + 1) % mod * (2 * x % mod + 1) % mod * inv6 % mod; }
int getsum(int n)
{
    if(n <= M) return sum[n];
    if(vis[n]) return vis[n];
    int res = fun1(n);
    for(int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        add(res, mod - (fun2(r) - fun2(l - 1) + mod) % mod * getsum(n / l) % mod);
    }
    return vis[n] = res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n, m; cin >> n >> m; init(M = min((int)powl(n, 0.67) + 10, N - 5));
    if(n > m) swap(n, m); int res = 0;
    for(int l = 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        add(res, (getsum(r) - getsum(l - 1) + mod) % mod * fun1(n / l) % mod * fun1(m / l) % mod);
    }
    cout << res << '\n';
    return 0;
}