数论分块（二）

上部分链接：数论分块（一）

建议没看过数论分块（一）的优先去看第一部分，第二部分题目难度较高。

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

给定 \(n, m\)，求：

\[\sum_{i = 1}^n\sum_{j = 1}^m\operatorname{lcm}(i, j)\pmod {20101009} \]

• 对于 \(100\%\) 的数据，保证 \(1\le n,m \le 10^7\)。

---

将式子变形：

\[\begin{aligned} \sum_{i = 1}^n\sum_{j = 1}^m\operatorname{lcm}(i, j) &= \sum_{i = 1}^n\sum_{j = 1}^m\frac{ij}{\gcd(i, j)} \\ &= \sum_{d = 1}^{n}\sum_{i = 1}^n\sum_{j = 1}^m\frac{ij}{d}[\gcd(i, j) = d] \\ &= \sum_{d = 1}^{n}\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{m}{d}\right\rfloor}ijd[\gcd(i, j) = 1] \\ &= \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}ij\sum_{k \mid \gcd(i, j)}\mu(k) \\ &= \sum_{d = 1}^{n}d\sum_{k = 1}\mu(k)\sum_{k \mid i}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{k \mid j}^{\left\lfloor\frac{m}{d}\right\rfloor}ij \\ &= \sum_{d = 1}^{n}d\sum_{k = 1}\mu(k)\sum_{i = 1}^{\left\lfloor\frac{n}{kd}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{m}{kd}\right\rfloor}ijk^2 \\ &= \sum_{d = 1}^{n}d\sum_{k = 1}\mu(k)k^2\sum_{i = 1}^{\left\lfloor\frac{n}{kd}\right\rfloor}i\sum_{j = 1}^{\left\lfloor\frac{m}{kd}\right\rfloor}j \\ &= \sum_{d = 1}^{n}d\sum_{k = 1}\mu(k)k^2\frac{\left(\left\lfloor\frac{n}{kd}\right\rfloor + 1\right)\left\lfloor\frac{n}{kd}\right\rfloor}{2}·\frac{\left(\left\lfloor\frac{m}{kd}\right\rfloor + 1\right)\left\lfloor\frac{m}{kd}\right\rfloor}{2} \\ \end{aligned} \]

令 \(T = kd\)，按 \(k\) 枚举改为按 \(T\) 枚举。

\[\begin{aligned} &= \sum_{d = 1}^{n}d\sum_{k = 1}\mu(k)k^2\frac{\left(\left\lfloor\frac{n}{kd}\right\rfloor + 1\right)\left\lfloor\frac{n}{kd}\right\rfloor\left(\left\lfloor\frac{m}{kd}\right\rfloor + 1\right)\left\lfloor\frac{m}{kd}\right\rfloor}{4} \\ &= \sum_{d = 1}^{n}d\sum_{T = 1, d \mid T}^{n}\mu\left(\frac{T}{d}\right)\left(\frac{T}{d}\right)^2\frac{\left(\left\lfloor\frac{n}{T}\right\rfloor + 1\right)\left\lfloor\frac{n}{T}\right\rfloor\left(\left\lfloor\frac{m}{T}\right\rfloor + 1\right)\left\lfloor\frac{m}{T}\right\rfloor}{4} \\ &= \sum_{T = 1}^{n}\sum_{d \mid T}T\mu\left(\frac{T}{d}\right)\left(\frac{T}{d}\right)\frac{\left(\left\lfloor\frac{n}{T}\right\rfloor + 1\right)\left\lfloor\frac{n}{T}\right\rfloor\left(\left\lfloor\frac{m}{T}\right\rfloor + 1\right)\left\lfloor\frac{m}{T}\right\rfloor}{4} \\ &= \sum_{T = 1}^{n}T\sum_{d \mid T}\mu(d)d\frac{\left(\left\lfloor\frac{n}{T}\right\rfloor + 1\right)\left\lfloor\frac{n}{T}\right\rfloor\left(\left\lfloor\frac{m}{T}\right\rfloor + 1\right)\left\lfloor\frac{m}{T}\right\rfloor}{4} \\ \end{aligned} \]

记 \(f(T) = \sum\limits_{d \mid T}d\mu(d)\)，根据直觉 \(f(T)\) 为积性函数，我们取几个值来看看。

\[f(1) = 1 \]

\[f(p) = 1 - p \]

\[f(p^k) = 1 - p \]

\[f(pq) = 1 - p - q + pq = (1 - p)(1 - q) \]

那么我们可以通过线性筛筛出 \(T \times f(T)\) 的前缀和。

后面一坨做数论分块，总体时间复杂度 \(O(n + \sqrt{n})\)。

PS: 本题还有杜教筛的做法，可以做到时间复杂度 \(O(n^{\frac{2}{3}})\)

参考代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7 + 10;
const int mod = 20101009;
int n, m, cnt, primes[N];
ll f[N];
bool st[N];

void init()
{
    f[1] = 1;
    for (int i = 2; i <= n; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, f[i] = 1 - i + mod;
        for (int j = 0; primes[j] * i <= n; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0)
            {
                f[primes[j] * i] = f[i];
                break;
            }
            f[primes[j] * i] = f[i] * (1 - primes[j] + mod) % mod;
        }
    }
    for (int i = 1; i <= n; i ++ ) f[i] = (f[i] * i % mod + f[i - 1]) % mod;
}

ll calc(ll n, ll m)
{
    return ((n + 1) * n / 2 % mod) * ((m + 1) * m / 2 % mod) % mod ;
}

int main()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    init();
    ll res = 0;
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (res += (f[r] - f[l - 1]) * calc(n / l, m / l) % mod) %= mod;
    }
    cout << (res + mod) % mod;
    return 0;
}

P3327 [SDOI2015] 约数个数和

设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(n,m\)，求

\[\sum_{i=1}^n\sum_{j=1}^md(ij) \]

对于 \(100\%\) 的数据，\(1\le T,n,m \le 50000\)。

---

\[\sum_{i=1}^n\sum_{j=1}^md(ij) = \sum_{i=1}^n\sum_{j=1}^m\sum_{d \mid ij} 1 \]

设 \(xy = d\)，考虑把 \(d\) 分解有什么情况，由于所有的 \(d\) 值是唯一的，所以每一个 \(x\) 和 \(y\) 的枚举需要唯一，不妨设 \(x \mid i\)，\(y \mid j\)，考虑 \(xy\) 什么时候唯一对应 \(ij\) 的因数。

可以发现，如果 \(xy\) 不能唯一对应 \(ij\) 的因数时，说明 \(xy\) 有多种拆解方式，而 \(x \mid i\)，\(y \mid j\)，因此，如果 \(x, y\) 含有共同的因子，这 \(xy\) 是可以再分的，因此所以它的充要条件是 \([\gcd(x, y) = 1]\)，原始可继续改写。

\[\sum_{i=1}^n\sum_{j=1}^m\sum_{d \mid ij} 1 = \sum_{i=1}^n\sum_{j=1}^m\sum_{x \mid i}\sum_{y \mid j}[\gcd(x, y) = 1] \]

又有 \(x, y\) 的求和是独立的，可以改写。

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m\sum_{x \mid i}\sum_{y \mid j}[\gcd(x, y) = 1] &= \sum_{x = 1}^n\sum_{y = 1}^m\sum_{x \mid i}^n\sum_{y \mid j}^m[\gcd(x, y) = 1] ] \\ &= \sum_{x = 1}^n\sum_{y = 1}^m\bigg\lfloor\frac{n}{x}\bigg\rfloor\bigg\lfloor\frac{m}{y}\bigg\rfloor[\gcd(x, y) = 1] \\ \end{aligned} \]

考虑莫比乌斯反演 \(\sum\limits_{d \mid n} \mu(d) = [n = 1]\)。

\[\begin{aligned} \sum_{x = 1}^n\sum_{y = 1}^m\bigg\lfloor\frac{n}{x}\bigg\rfloor\bigg\lfloor\frac{m}{y}\bigg\rfloor[\gcd(x, y) = 1] &= \sum_{x = 1}^n\sum_{y = 1}^m\bigg\lfloor\frac{n}{x}\bigg\rfloor\bigg\lfloor\frac{m}{y}\bigg\rfloor\sum_{d \mid \gcd(x, y)}\mu(d) \\ &= \sum_{d = 1}^n\mu(d)\sum_{d \mid x}^n\sum_{d \mid y}^m\bigg\lfloor\frac{n}{x}\bigg\rfloor\bigg\lfloor\frac{m}{y}\bigg\rfloor \\ &= \sum_{d = 1}^n\mu(d)\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}\bigg\lfloor\frac{n}{di}\bigg\rfloor\bigg\lfloor\frac{m}{dj}\bigg\rfloor \\ &= \sum_{d = 1}^n\mu(d)\left(\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\bigg\lfloor\frac{n}{di}\bigg\rfloor\right)\left(\sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}\bigg\lfloor\frac{m}{dj}\bigg\rfloor\right) \end{aligned} \]

设 \(S(n) = \sum\limits_{i = 1}^{n}\lfloor\frac{n}{i}\rfloor\)，则原式可替换为：

\[\sum_{d = 1}^n\mu(d)S\left(\bigg\lfloor\frac{n}{d}\bigg\rfloor\right)S\left(\bigg\lfloor\frac{m}{d}\bigg\rfloor\right) \]

预处理一下 \(\mu(n)\) 的前缀和与 \(S(n)\) 的值，进行数论分块即可。

参考代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e4 + 10;
int cnt, st[N], primes[N];
int n, m;
ll S[N], mu[N];

void init()
{
    mu[1] = 1;
    for (ll i = 2; i < N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
        for (int j = 0; primes[j] * i < N; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0) break;
            mu[primes[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
    {
        mu[i] += mu[i - 1];
        for (int l = 1, r; l <= i; l = r + 1)
        {
            r = i / (i / l);
            S[i] += 1ll * (r - l + 1) * (i / l);
        }
    }
}

ll calc(int n, int m)
{
    ll res = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        res += (mu[r] - mu[l - 1]) * S[n / l] * S[m / l];
    }
    return res;
}

void solve()
{
    scanf("%d%d", &n, &m);
    printf("%lld\n", calc(n, m));
}

int main()
{
    init();
    int T;
    scanf("%d", &T);
    while (T -- ) solve();
    return 0;
}

P4466 [国家集训队] 和与积

考虑 \(a + b \mid ab\) 有怎样的性质，设 \(d = \gcd(a, b), a = id, b = jd\)，则原式可化为 \(i + j \mid ijd\)，由于 \(\gcd(i, j) = \gcd(i + j, j) = \gcd(i, i + j) = 1\)，故可知 \(i + j \mid d\)。

因此我们的目标是对这样的 \(d\) 的个数求和。

\[\sum_{i=1}^n\sum_{j=1}^{i-1}\left\lfloor\frac{d}{i + j}\right\rfloor[\gcd(i, j) = 1][id \le n] \]

合法的 \(d\) 的上界为 \(\left\lfloor\frac{n}{i}\right \rfloor\)，所以

\[\sum_{i=1}^n\sum_{j=1}^{i-1}\left\lfloor\frac{\left\lfloor\frac{n}{i}\right\rfloor}{i + j}\right\rfloor[\gcd(i, j) = 1] = \sum_{i=1}^n\sum_{j=1}^{i-1}\left\lfloor\frac{n}{i(i + j)}\right\rfloor[\gcd(i, j) = 1] \]

可以发现 \(i\) 的上界不超过 \(\sqrt{n}\)，可以缩小枚举范围，并对后面来一发莫反。

\[\begin{aligned} \sum_{i=1}^{\sqrt{n}}\sum_{j=1}^{i-1}\left\lfloor\frac{n}{i(i + j)}\right\rfloor[\gcd(i, j) = 1] &= \sum_{i=1}^{\sqrt{n}}\sum_{j=1}^{i-1}\left\lfloor\frac{n}{i(i + j)}\right\rfloor\sum_{d \mid \gcd(i,j)}\mu(d) \\ &= \sum_{i=1}^{\sqrt{n}}\sum_{j=1}^{i-1}\left\lfloor\frac{n}{i(i + j)}\right\rfloor\sum_{d \mid i, d \mid j}\mu(d) \\ &= \sum_{d=1}^{\sqrt{n}}\mu(d)\sum_{i=1}^{\left\lfloor\frac{\sqrt{n}}{d}\right\rfloor}\sum_{j=1}^{i-1}\left\lfloor\frac{\left\lfloor\frac{n}{id^2}\right\rfloor}{i + j}\right\rfloor \\ &= \sum_{d=1}^{\sqrt{n}}\mu(d)\sum_{i=1}^{\left\lfloor\frac{\sqrt{n}}{d}\right\rfloor}\sum_{j=1}^{i-1}\left\lfloor\frac{\left\lfloor\frac{n}{id^2}\right\rfloor}{i + j}\right\rfloor \\ &= \sum_{d=1}^{\sqrt{n}}\mu(d)\sum_{i=1}^{\left\lfloor\frac{\sqrt{n}}{d}\right\rfloor}\sum_{j=i+1}^{2i-1}\left\lfloor\frac{\left\lfloor\frac{n}{id^2}\right\rfloor}{j}\right\rfloor \end{aligned} \]

后面可以数论分块，对复杂度积分（舍去常数）有：

\[\begin{aligned} \int_1^{\sqrt{n}}dx\int_1^{x}\sqrt{\frac{n}{x^2y}}dy \to \int_1^{\sqrt{n}}k\sqrt{\frac{n}{x}}dx \to k\sqrt{n}(\sqrt{n})^{\frac{1}{2}} \to n^\frac{3}{4} \end{aligned} \]

因此复杂度的上界为 \(O(n^\frac{3}{4})\)，可以通过。

加了一些剪枝，398ms 的最优解。

参考代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 60000;
int n;
int primes[N], mu[N], st[N], cnt;

void init(int maxn)
{
    mu[1] = 1;
    for (int i = 2; i <= maxn; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
        for (int j = 0; i * primes[j] <= maxn; j ++ )
        {
            st[i * primes[j]] = 1;
            if (i % primes[j] == 0) break;
            mu[i * primes[j]] = -mu[i];
        }
    }
}

int calc(int st, int ed, int mul)
{
    int res = 0;
    ed = min(ed, mul);
    for (int l = st, r; l <= ed; l = r + 1)
    {
        r = min(ed, mul / (mul / l));
        res += (r - l + 1) * (mul / l);
    }
    return res;
}

int main()
{
    cin >> n;
    int upper = sqrt(n);
    init(upper);
    ll res = 0;
    for (int d = 1; d <= upper; d ++ )
    {
        if (!mu[d]) continue;
        for (int i = 1; i * d <= upper; i ++ )
            res += mu[d] * calc(i + 1, (i << 1) - 1, n / (d * d * i));
    }
    cout << res;
    return 0;
}

P5572 [CmdOI2019] 简单的数论题

给出 \(n,m\) 求下列式子的值 ：

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m \varphi\left(\dfrac{{\rm lcm}(i,j)}{\gcd(i,j)}\right) \bmod 23333 \]

对于所有测试点， \(T\leq 3\times 10^4,\ m\leq n\leq 5\times 10^4\)。

---

先看看这个式子的真面目。

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m \varphi\left(\dfrac{{\rm lcm}(i,j)}{\gcd(i,j)}\right) &= \sum_{i=1}^n\sum_{j=1}^m \varphi\left(\dfrac{ij}{\gcd^{2}(i,j)}\right) \\ &= \sum_{d=1}^{n}\sum_{i=1}^n\sum_{j=1}^m \varphi\left(\dfrac{ij}{d^2}\right)[\gcd(i, j)=d] \\ &= \sum_{d=1}^{n}\sum_{i=1}^n\sum_{j=1}^m \varphi\left(\dfrac{i}{d}\right)\varphi\left(\dfrac{i}{d}\right)[\gcd(i, j)=d] \\ &= \sum_{d=1}^{n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor} \varphi(i)\varphi(j)[\gcd(i, j)=1] \\ &= \sum_{d=1}^{n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor} \varphi(i)\varphi(j)\sum_{k\mid{\gcd(i,j)}}\mu(k) \\ &= \sum_{k=1}^n\mu(k)\sum_{d=1}^{n}\sum_{k \mid i}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{k \mid j}^{\left\lfloor\frac{m}{d}\right\rfloor} \varphi(i)\varphi(j) \\ &= \sum_{k=1}^n\mu(k)\sum_{d=1}^{n}\sum_{i = 1}^{\left\lfloor\frac{n}{dk}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{m}{dk}\right\rfloor} \varphi(ki)\varphi(kj) \\ &= \sum_{k=1}^n\mu(k)\sum_{d=1}^{n}\left(\sum_{i = 1}^{\left\lfloor\frac{n}{dk}\right\rfloor} \varphi(ki)\right)\left(\sum_{j = 1}^{\left\lfloor\frac{m}{dk}\right\rfloor}\varphi(kj)\right) \\ &= \sum_{k=1}^n\sum_{d=1}^n\mu(d)\left(\sum_{i = 1}^{\left\lfloor\frac{n}{dk}\right\rfloor} \varphi(di)\right)\left(\sum_{j = 1}^{\left\lfloor\frac{m}{dk}\right\rfloor}\varphi(dj)\right) \end{aligned} \]

直到这里，步骤和上一题几乎完全一致，我们考虑用 \(T = dk\) 换元。

\[\sum_{k=1}^n\sum_{d=1}^n\mu(d)\left(\sum_{i = 1}^{\left\lfloor\frac{n}{dk}\right\rfloor} \varphi(di)\right)\left(\sum_{j = 1}^{\left\lfloor\frac{m}{dk}\right\rfloor}\varphi(dj)\right) = \sum_{T=1}^{n}\sum_{d\mid T}\mu(d)\left(\sum_{i = 1}^{\left\lfloor\frac{n}{T}\right\rfloor} \varphi(di)\right)\left(\sum_{j = 1}^{\left\lfloor\frac{m}{T}\right\rfloor}\varphi(dj)\right) \]

观察式子后两坨满足 \(f(d, s) = \sum\limits_{i=1}^s\varphi(di)\)，预处理的复杂度满足调和级数，总复杂度 \(O(n\ln{n})\)。

记 \(g(x, y, z) = \sum\limits_{T=1}^{z}\sum\limits_{d\mid T}\mu(d)\left(\sum\limits_{i = 1}^{x} \varphi(di)\right)\left(\sum\limits_{j = 1}^{y}\varphi(dj)\right)\)，把答案用端点 \(\left\lfloor\frac{n}{i}\right\rfloor\)，\(\left\lfloor\frac{m}{i}\right\rfloor\) 拆分成若干线段 \([l, r]\)，一个线段的答案贡献为 \(g(\left\lfloor\frac{n}{i}\right\rfloor, \left\lfloor\frac{m}{i}\right\rfloor, r) - g(\left\lfloor\frac{n}{i}\right\rfloor, \left\lfloor\frac{m}{i}\right\rfloor, l - 1)\)，但是令人悲伤的是，预处理复杂度是 \(O(n^2m\log{n})\) 的，无法接受，没办法直接上数论分块。

如果我们暴力计算式子，则对于每次询问我们的复杂度是 \(O(n\log{n})\) 的，也无法通过，于是我们想到用根号分治来解决这个问题。

设置阈值 \(B\)，对 \(\left\lfloor\frac{n}{B}\right\rfloor\) 以内的 \(z\) 暴力遍历，时间复杂度是 \(O\left(\frac{n}{B}\log{\frac{n}{B}}\right)\)，其余的 \(g(x, y, z)\) 使用数论分块，将预处理时间分摊为 \(O(B^2n\log{n})\)，总的时间复杂度应为 \(O\left(T\frac{n}{B}\log{\frac{n}{B}}+T\sqrt{B}\log{n} + B^2n\log{n}\right)\)，取 \(B\) 取 \(20\) 到 \(100\) 中任意一个数都可以，取 \(50\) 比较快。

参考代码

#include<bits/stdc++.h>
using namespace std;
const int N = 5e4 + 10, B = 50, mod = 23333;
int n, m;
int st[N], prime[N], phi[N], mu[N], cnt;
vector<int> factor[N], f[N], g[B + 5][B + 5];

int calc(int x, int y, int z)
{
    int res = 0;
    for (auto t : factor[z])
        (res += mu[t] * f[t][x] % mod * f[t][y] % mod + mod) %= mod;
    return res;
}

void init()
{
    phi[1] = mu[1] = 1;
    for (int i = 2; i < N; i ++ )
    {
        if (!st[i]) prime[cnt ++ ] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 0; prime[j] * i < N; j ++ )
        {
            st[prime[j] * i] = 1;
            if (i % prime[j] == 0)
            {
                phi[prime[j] * i] = phi[i] * prime[j];
                break;
            }
            phi[prime[j] * i] = phi[i] * (prime[j] - 1);
            mu[prime[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
        for (int j = i; j < N; j += i)
            factor[j].push_back(i);
    for (int i = 1; i < N; i ++ )
    {
        f[i].push_back(0);
        for (int j = 1; i * j < N; j ++ )
            f[i].push_back((f[i][j - 1] + phi[i * j]) % mod);
    }
    for (int i = 1; i <= B; i ++ )
    {
        for (int j = 1; j <= B; j ++ )
        {
            g[i][j].push_back(0);
            for (int k = 1; k * max(i, j) < N; k ++ )
                g[i][j].push_back((g[i][j][k - 1] + calc(i, j, k)) % mod);
        }
    }
}

void solve()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1; i * B <= m; i ++ ) (ans += calc(n / i, m / i, i)) %= mod;
    for (int l = m / B + 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (ans += g[n / l][m / l][r] - g[n / l][m / l][l - 1] + mod) %= mod;
    }
    cout << ans << endl;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL), cout.tie(NULL);
    init();
    int T;
    cin >> T;
    while (T -- ) solve();
    return 0;
}

P4240 毒瘤之神的考验

\(T\) 次询问，每次给定 \(n, m\)，求：

\[\left( \sum_{i=1}^n \sum_{j=1}^m \varphi(ij) \right) \bmod 998244353 \]

对于 \(100\%\) 的数据，\(1 \le T \le {10}^4\)，\(1 \le n, m \le {10}^5\)。

---

对于 \(\varphi(ij)\) 有（\(p\) 为质数）：

\[\begin{aligned} \varphi(ij) &= ij\prod_{p \mid ij}\frac{p - 1}{p} \\ &= ij\frac{\prod_{p \mid i}\frac{p - 1}{p}\prod_{p \mid j}\frac{p - 1}{p}}{\prod_{p \mid i, p \mid j}\frac{p - 1}{p}} \\ &= \frac{i\prod_{p \mid i}\frac{p - 1}{p}j\prod_{p \mid j}\frac{p - 1}{p}}{\prod_{p \mid \gcd(i, j)}\frac{p - 1}{p}} \\ &= \frac{\varphi(i)\varphi(j)\gcd(i, j)}{\varphi(\gcd(i, j))} \\ \end{aligned} \]

一如既往的化式子：

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m\varphi(ij) &= \sum_{i=1}^n\sum_{j=1}^m\frac{\varphi(i)\varphi(j)\gcd(i, j)}{\varphi(\gcd(i, j))} \\ &= \sum_{d=1}^{n}\sum_{d\mid i}^n\sum_{d\mid j}^m\frac{\varphi(i)\varphi(j)d}{\varphi(d)}[\gcd(i, j)=d] \\ &= \sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\varphi(id)\varphi(jd)[\gcd(i, j)=1] \\ &= \sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\varphi(id)\varphi(jd)\sum_{k\mid\gcd(i,j)}\mu(k) \\ &= \sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{k=1}^n\mu(k)\sum_{k\mid i}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{k\mid j}^{\left\lfloor\frac{m}{d}\right\rfloor}\varphi(id)\varphi(jd) \\ &= \sum_{d=1}^{n}\frac{d}{\varphi(d)}\sum_{k=1}^n\mu(k)\sum_{i=1}^{\left\lfloor\frac{n}{kd}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{kd}\right\rfloor}\varphi(idk)\varphi(jdk) \\ \end{aligned} \]

记 \(T = kd\) 有：

\[\sum_{T=1}^{n}\sum_{d\mid T}\frac{d}{\varphi(d)}\mu\left(\frac{T}{d}\right)\sum_{i=1}^{\left\lfloor\frac{n}{T}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{T}\right\rfloor}\varphi(Ti)\varphi(Tj) \]

这和上一个题也是类似的，记 \(f(d, s) = \sum\limits_{i=1}^s\varphi(di)\)，用 \(O(n\ln{n})\) 的时间复杂度预处理出。

记 \(g(x, y, z) = \sum\limits_{T=1}^z\sum\limits_{d\mid T}\frac{d}{\varphi(d)}\mu\left(\frac{T}{d}\right)\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{y}\varphi(Ti)\varphi(Tj)\)

考虑根号分治，设阈值为 \(B\)，爆算 \(\left\lfloor\frac{n}{B}\right\rfloor\) 以内的 \(z\)，预处理出剩下的部分，复杂度和上一题完全一致，甚至式子都是这么的形似，时间复杂度 \(O\left(T\frac{n}{B}\log{\frac{n}{B}}+T\sqrt{B}\log{n} + B^2n\log{n}\right)\)，取 \(B\) 为 \(45\) 左右即可。

参考代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5 + 10, B = 45, mod = 998244353;
int n, m, st[N], primes[N], phi[N], mu[N], inv[N], cnt;
vector<int> factor[N], f[N], g[B + 5][B + 5];

int calc(int x, int y, int z)
{
    int res = 0;
    for (auto t : factor[z])
        (res += 1ll * t * inv[phi[t]] * mu[z / t] % mod) %= mod;
    res = (1ll * res * f[z][x] % mod * f[z][y] % mod + mod) % mod;
    return res;
}

void init()
{
    inv[1] = phi[1] = mu[1] = 1;
    for (int i = 2; i < N; i ++ ) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
    for (int i = 2; i < N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 0; 1ll * primes[j] * i < N; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0)
            {
                phi[primes[j] * i] = phi[i] * primes[j];
                break;
            }
            phi[primes[j] * i] = phi[i] * (primes[j] - 1);
            mu[primes[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
        for (int j = i; j < N; j += i)
            factor[j].push_back(i);
    for (int i = 1; i < N; i ++ )
    {
        f[i].push_back(0);
        for (int j = 1; i * j < N; j ++ )
            f[i].push_back((f[i][j - 1] + phi[i * j]) % mod);
    }
    for (int i = 1; i <= B; i ++ )
    {
        for (int j = 1; j <= B; j ++ )
        {
            g[i][j].push_back(0);
            for (int k = 1; k * max(i, j) < N; k ++ )
                g[i][j].push_back((g[i][j][k - 1] + calc(i, j, k)) % mod);
        }
    }
}

void solve()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1; i * B <= m; i ++ ) (ans += calc(n / i, m / i, i)) %= mod;
    for (int l = m / B + 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (ans += (g[n / l][m / l][r] - g[n / l][m / l][l - 1]) % mod) %= mod;
    }
    cout << (ans + mod) % mod << '\n';
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);cout.tie(NULL);
    init();
    int T;
    cin >> T;
    while (T -- ) solve();
    return 0;
}

我们发现其实不用每次都去枚举质因数，存储进一个系数数组即可，这样时间复杂度少一个 \(\log(n)\) 级别。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5 + 10, B = 45, mod = 998244353;
int n, m, st[N], primes[N], phi[N], mu[N], inv[N], coff[N], cnt;
vector<int> f[N], g[B + 5][B + 5];

void init()
{
    inv[1] = phi[1] = mu[1] = 1;
    for (int i = 2; i < N; i ++ ) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
    for (int i = 2; i < N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 0; primes[j] * i < N; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0)
            {
                phi[primes[j] * i] = phi[i] * primes[j];
                break;
            }
            phi[primes[j] * i] = phi[i] * (primes[j] - 1);
            mu[primes[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
        for (int j = i; j < N; j += i)
            (coff[j] += (1ll * mu[j / i] * i * inv[phi[i]] % mod + mod) % mod) %= mod;
    for (int i = 1; i < N; i ++ )
    {
        f[i].push_back(0);
        for (int j = 1; i * j < N; j ++ )
            f[i].push_back((f[i][j - 1] + phi[i * j]) % mod);
    }
    for (int i = 1; i <= B; i ++ )
    {
        for (int j = 1; j <= B; j ++ )
        {
            g[i][j].push_back(0);
            for (int k = 1; k * max(i, j) < N; k ++ )
                g[i][j].push_back((g[i][j][k - 1] + 1ll * coff[k] * f[k][i] % mod * f[k][j] % mod) % mod);
        }
    }
}

void solve()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1; i * B <= m; i ++ ) (ans += 1ll * coff[i] * f[i][n / i] % mod * f[i][m / i] % mod) %= mod;
    for (int l = m / B + 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (ans += ((g[n / l][m / l][r] - g[n / l][m / l][l - 1]) % mod + mod) % mod) %= mod;
    }
    cout << ans << '\n';
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);cout.tie(NULL);
    init();
    int T;
    cin >> T;
    while (T -- ) solve();
    return 0;
}