[BZOJ2693]jzptab

[BZOJ2693]jzptab

试题描述

给出 \(n, m\)，求 \(\sum_{i=1}^n \sum_{j=1}^m \mathrm{lcm}(i, j)\) 对 \(100000009\) 取模后的结果。

多组询问。

输入

一个正整数 \(T\) 表示数据组数

接下来 \(T\) 行，每行两个正整数，表示 \(n\)、\(m\)

输出

\(T\) 行，每行一个整数，表示第 \(i\) 组数据的结果

输入示例

1
4 5

输出示例

122

数据规模及约定

\(T \le 10000\)

\(n, m \le 10000000\)

题解

这道题是那道题的加强版，强化的地方是多组数据而不是 \(n, m\) 的范围，所以那篇文章后面的优化做法是没有用的。

那篇博文的公式有点丑，所以在这里重推一下……

\[ \sum_{i=1}^n \sum_{j=1}^m \mathrm{lcm}(i, j) \\ = \sum_{i=1}^n \sum_{j=1}^m \frac{ij}{\gcd(i, j)} \]

现在考虑枚举最大公约数

\[ = \sum_{g=1}^{\min \{n, m\} } \frac{1}{g} \sum_{i=1}^n i \sum_{j=1}^m j \cdot [\gcd(i, j) = g] \\ = \sum_{g=1}^{\min \{n, m\} } \frac{1}{g} \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} ig \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} jg \sum_{d|i, d|j} \mu(d) \\ = \sum_{g=1}^{\min \{n, m\} } g \sum_{d=1}^{\lfloor \frac{\min \{n, m\}}{g} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{n}{gd} \rfloor} id \sum_{j=1}^{\lfloor \frac{m}{gd} \rfloor} jd \\ = \sum_{g=1}^{\min \{n, m\} } g \sum_{d=1}^{\lfloor \frac{\min \{n, m\}}{g} \rfloor} \mu(d) \cdot d^2 \frac{\lfloor \frac{n}{gd} \rfloor \left( \lfloor \frac{n}{gd} \rfloor + 1 \right)}{2} \cdot \frac{\lfloor \frac{m}{gd} \rfloor \left( \lfloor \frac{m}{gd} \rfloor + 1 \right)}{2} \\ = \frac{1}{4} \sum_{g=1}^{\min \{n, m\} } g \sum_{d=1}^{\lfloor \frac{\min \{n, m\}}{g} \rfloor} \mu(d) \cdot d^2 \lfloor \frac{n}{gd} \rfloor \left( \lfloor \frac{n}{gd} \rfloor + 1 \right) \lfloor \frac{m}{gd} \rfloor \left( \lfloor \frac{m}{gd} \rfloor + 1 \right) \]

从这里开始将和那道题不一样了，现在用 \(T\) 替换 \(gd\)

\[ = \frac{1}{4} \sum_{T=1}^{\min \{n, m\} } \lfloor \frac{n}{T} \rfloor \left( \lfloor \frac{n}{T} \rfloor + 1 \right) \lfloor \frac{m}{T} \rfloor \left( \lfloor \frac{m}{T} \rfloor + 1 \right) \sum_{d|T} \mu(d) \cdot d^2 \cdot \frac{T}{d} \\ = \frac{1}{4} \sum_{T=1}^{\min \{n, m\} } \lfloor \frac{n}{T} \rfloor \left( \lfloor \frac{n}{T} \rfloor + 1 \right) \lfloor \frac{m}{T} \rfloor \left( \lfloor \frac{m}{T} \rfloor + 1 \right) T \sum_{d|T} \mu(d) \cdot d \]

令 \(f(T) = T \sum_{d|T} \mu(d) \cdot d\)，则 \(f(T)\) 是积性函数（两个积性函数的乘积是积性函数，两个积性函数的狄利克雷卷积也是积性函数），我们可以线性筛出来。

如何线性筛？假设 \(p\) 是质数，则 \(f(p) = p(1 - p)\)，我们要解决的就是已知 \(f(T)\) 求 \(f(Tp)\)，分两种情况：

• \(p|T\)，令 \(T = T' p^k\)，且 \(T'\) 中不含因子 \(p\)，那么 \(f(Tp) = f(T') \cdot f(p^{k+1})\)，在线性筛的过程 \(f(T')\) 已知，\(f(p^k)\) 已知，那么 \(f(p^{k+1}) = f(p^k) \cdot p\)，这是因为 \(f(p^k) = p^k \sum_{i=0}^k \mu(p^i) p^i\)，后面 \(\mu(p^i)(i > 1)\) 都是 \(0\)，所以只有前面那个 \(p^k\) 会乘上一个 \(p\)，所以最后得到 \(f(Tp) = f(T) \cdot p\)；
• \(p \nmid T\)，显然 \(f(Tp) = f(T) \cdot f(p)\)。

这样我们可以 \(O(n)\) 的时间求出 \(f(T)\) 的前缀和，后面的询问就变成 \(O(\sqrt n)\) 一次啦。

注意 \(100000009\) 不是素数，但这并不影响我们求 \(4\) 的逆元。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;

#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
    int x = 0, f = 1; char c = getchar();
    while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
    while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
    return x * f;
}

#define maxn 10000010
#define MOD 100000009
#define inv4 75000007
#define LL long long

bool vis[maxn];
int prime[maxn], cp, f[maxn], sumf[maxn];

void init() {
    int n = maxn - 1;
    f[1] = sumf[1] = 1;
    rep(i, 2, n) {
        if(!vis[i]) prime[++cp] = i, f[i] = (LL)i * (1 - i + MOD) % MOD;
        for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
            vis[i*prime[j]] = 1;
            if(i % prime[j] == 0) {
                f[i*prime[j]] = (LL)f[i] * prime[j] % MOD;
                break;
            }
            f[i*prime[j]] = (LL)f[i] * f[prime[j]] % MOD;
        }
        sumf[i] = sumf[i-1] + f[i];
        if(sumf[i] >= MOD) sumf[i] -= MOD;
    }
    return ;
}

void work() {
    int n = read(), m = read(), ans = 0;
    for(int T = 1; T <= min(n, m); ) {
        int r = min(min(n / (n / T), m / (m / T)), min(n, m));
        ans += (LL)(n / T) * (n / T + 1) % MOD * (m / T) % MOD * (m / T + 1) % MOD * (sumf[r] - sumf[T-1] + MOD) % MOD;
        if(ans >= MOD) ans -= MOD;
        T = r + 1;
    }
    ans = (LL)ans * inv4 % MOD;
    printf("%d\n", ans);
    return ;
}

int main() {
    init();

    int T = read();
    while(T--) work();

    return 0;
}

转载于:https://www.cnblogs.com/xiao-ju-ruo-xjr/p/8462102.html