题目链接：XJOI 1739B

给出 $n$，求：

$$ \left(\sum_{i = 1}^{n} \sum_{j = 1}^{n} \sum_{k = 1}^{n} \varphi(ij) d(jk) \mu(ki)\right) \bmod 998244353 $$

其中 $d(x)$ 表示 $x$ 的约数个数。

数据范围：$1 \le n \le 5 \cdot 10^4$。

Solution

首先，将 $\varphi, d, \mu$ 都拆开：

$$ \begin{align*} \varphi(ij) & = \frac{\varphi(i)\varphi(j)(i, j)}{\varphi((i, j))} \\ d(ij) & = \sum_{u \mid i} \sum_{v \mid j} [(u, v) = 1] = \sum_{p \mid (i, j)} \mu(p) d\left(\frac{i}{p}\right) d\left(\frac{j}{p}\right) \\ \mu(ij) & = \mu(i) \mu(j) [(i, j) = 1] = \mu(i) \mu(j) \sum_{p \mid (i, j)} \mu(p) \end{align*} $$

将所有的东西都代回去：

$$ \sum_{i = 1}^{n} \sum_{j = 1}^{n} \sum_{k = 1}^{n} \varphi(ij) d(jk) \mu(ki) \\ \sum_{i = 1}^{n} \sum_{j = 1}^{n} \sum_{k = 1}^{n} \frac{\varphi(i)\varphi(j)(i, j)}{\varphi((i, j))} \sum_{u \mid (j, k)} \mu(u) d\left(\frac{j}{u}\right) d\left(\frac{k}{u}\right) \mu(k) \mu(i) \sum_{v \mid (k, i)} \mu(v) $$

将只与 $i$ 有关的值整理成 $f(i) = \varphi(i) \mu(i)$，关于 $(i, j)$ 的值整理成 $g(x) = \frac{x}{\varphi(x)}$，得：

$$ \sum_{i = 1}^{n} \sum_{j = 1}^{n} \sum_{k = 1}^{n} f(i) \mu(k) \varphi(j) g((i, j)) \sum_{u \mid (j, k)} \mu(u) d\left(\frac{j}{u}\right) d\left(\frac{k}{u}\right) \sum_{v \mid (k, i)} \mu(v) $$

当前已经有了 $u \mid (j, k)$ 和 $v \mid (k, i)$，可以配一个 $w \mid (i, j)$。具体地，设 $h = g \ast \mu$ 则有 $g = h \ast 1$，得到：

$$ \sum_{i = 1}^{n} \sum_{j = 1}^{n} \sum_{k = 1}^{n} f(i) \mu(k) \varphi(j) \sum_{w \mid (i, j)} h(w) \sum_{u \mid (j, k)} \mu(u) d\left(\frac{j}{u}\right) d\left(\frac{k}{u}\right) \sum_{v \mid (k, i)} \mu(v) $$

变换求和顺序：

$$ \sum_{u = 1}^{n} \sum_{v = 1}^{n} \sum_{w = 1}^{n} \mu(u) \mu(v) h(w) \left(\sum_{\operatorname{lcm}(u, v) \mid k} d\left(\frac{k}{u}\right) \mu(k)\right) \left(\sum_{\operatorname{lcm}(v, w) \mid i} f(i)\right) \left(\sum_{\operatorname{lcm}(w, u) \mid j} d\left(\frac{j}{u}\right) \varphi(j)\right) $$

分别设后面的函数为 $A(u, \operatorname{lcm}(u, v)), B(\operatorname{lcm}(v, w)), C(u, \operatorname{lcm}(w, u))$，这三个函数都可以预处理求出。

套用 「SDOI 2018」旧试题 的方法进行计算。额外地，可以注意到此处的 $u, v$ 是「显性地」满足 $\mu(u) \ne 0, mu(v) \ne 0$ 的，否则无意义。此处的 $w$ 由于满足 $w \mid (i, j)$ 且原式中存在 $\mu(ki)$ 项，因此不满足 $\mu(w) \ne 0$ 的 $w$ 也是无意义的。

对于 $A, C$ 的计算，需要枚举 $u$、$u$ 的倍数 $\operatorname{lcm}(u, v)$ 和 $\operatorname{lcm}$ 的倍数 $k$，枚举次数为：

$$ \begin{align*} & \sum_{i = 1}^{n} \sum_{j \mid i} \sum_{k \mid j} 1 \\ = & \sum_{k = 1}^{n} \sum_{k \mid j} \left\lfloor\frac{n}{j}\right\rfloor \\ = & \sum_{k = 1}^{n} \sum_{j = 1}^{\left\lfloor\frac{n}{k}\right\rfloor} \left\lfloor\frac{n}{jk}\right\rfloor \\ = & \sum_{k = 1}^{n} \left\lfloor\frac{n}{k}\right\rfloor \ln \frac{n}{k} \\ \end{align*} $$

故该部分的时间复杂度为：

$$ \begin{align*} & \mathcal O\left(\sum_{i = 1}^{n} \sum_{j \mid i} \sum_{k \mid j} 1\right) \\ = & \mathcal O\left(\int_{1}^{n} \frac{n}{x}\ln\frac{n}{x} \, \mathrm{d}x\right) \\ = & \mathcal O(n \log^{2} n) \end{align*} $$

时间复杂度：$\mathcal O(n \log^{2} n + m \sqrt m)$。

Code

#include <bits/stdc++.h>
typedef unsigned int uint;
typedef unsigned long long uint64;

const uint MOD = 998244353;

inline uint norm(const uint v) {
    return v >= MOD ? v - MOD : v;
}
struct Z {
    uint v;
    inline Z(const uint v = 0) : v(v) {}
};
inline Z operator+(const Z x, const Z y) {
    return norm(x.v + y.v);
}
inline Z operator-(const Z x, const Z y) {
    return norm(x.v + MOD - y.v);
}
inline Z operator*(const Z x, const Z y) {
    return static_cast<uint64>(x.v) * y.v % MOD;
}
inline Z &operator+=(Z &x, const Z y) {
    return x = x + y;
}
inline Z &operator-=(Z &x, const Z y) {
    return x = x - y;
}
inline Z &operator*=(Z &x, const Z y) {
    return x = x * y;
}
inline Z operator-(const Z x) {
    return x.v == 0 ? 0 : MOD - x.v;
}
inline Z power(Z x, uint64 k) {
    Z ans = 1;
    for (; k > 0; k >>= 1, x *= x) {
        if (k & 1) {
            ans *= x;
        }
    }
    return ans;
}
inline Z inver(Z x) {
    return power(x, MOD - 2);
}

const int N = 5e4 + 5, M = 5e5 + 5;
int n, tot, cnt, p[N], prime[N], deg[N], vis[N], d[N], phi[N];
Z mu[N], h[N], A[N];
std::vector<std::pair<int, int>> E[N];
std::vector<Z> B[N], C[N];
struct Edge {
    int u, v, w;
} e[M];

void init(int n) {
    std::bitset<N> flg;
    flg.set();
    mu[1] = phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (flg[i]) {
            prime[++tot] = i;
            mu[i] = MOD - 1;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= tot && i * prime[j] <= n; j++) {
            int x = i * prime[j];
            flg[x] = 0;
            if (i % prime[j] == 0) {
                mu[x] = 0;
                phi[x] = phi[i] * prime[j];
                break;
            }
            mu[x] = -mu[i];
            phi[x] = phi[i] * (prime[j] - 1);
        }
    }
    for (int i = 1; i <= n; i++) {
        if (mu[i].v) {
            p[++cnt] = i;
        }
    }
    p[cnt + 1] = INT_MAX;
    for (int i = 1; i <= n; i++) {
        Z g = i * inver(phi[i]);
        for (int j = i; j <= n; j += i) {
            h[j] += g * mu[j / i];
            A[i] += mu[j] * phi[j];
            d[j]++;
        }
    }
    for (int i = 1; i <= n; i++) {
        B[i].resize(n / i + 1);
        C[i].resize(n / i + 1);
        for (int j = i; j <= n; j += i) {
            Z sum1 = 0, sum2 = 0;
            for (int k = j; k <= n; k += j) {
                sum1 += d[k / i] * phi[k];
                sum2 += d[k / i] * mu[k];
            }
            B[i][j / i] = sum1;
            C[i][j / i] = sum2;
        }
    }
}
inline Z calc(int i, int j, int k, int ij, int jk, int ki) {
    return mu[i] * h[j] * mu[k] * A[jk] * B[i][ij / i] * C[i][ki / i];
}
int main() {
    scanf("%d", &n);
    init(n);
    int cntE = 0;
    for (int g = 1; g <= cnt; g++) {
        for (int i = 1; p[g] * p[i] <= n; i++) {
            if (mu[p[g] * p[i]].v) {
                int lim = n / (p[g] * p[i]);
                for (int j = i + 1; p[j] <= lim; j++) {
                    if (mu[p[g] * p[j]].v && std::__gcd(p[i], p[j]) == 1) {
                        int u = p[g] * p[i], v = p[g] * p[j], w = p[g] * p[i] * p[j];
                        deg[u]++;
                        deg[v]++;
                        e[++cntE] = Edge{u, v, w};
                    }
                }
            }
        }
    }
    Z ans = 0;
    for (int i = 1; i <= cnt; i++) {
        ans += h[p[i]] * A[p[i]] * B[p[i]][1] * C[p[i]][1];
    }
    for (int i = 1; i <= cntE; i++) {
        int u = e[i].u, v = e[i].v, w = e[i].w;
        ans += calc(v, u, u, w, u, w) + calc(u, v, u, w, w, u) + calc(u, u, v, u, w, w);
        ans += calc(u, v, v, w, v, w) + calc(v, u, v, w, w, v) + calc(v, v, u, v, w, w);
    }
    for (int i = 1; i <= cntE; i++) {
        int u = e[i].u, v = e[i].v;
        if (deg[u] < deg[v] || (deg[u] == deg[v] && u > v)) {
            std::swap(u, v);
        }
        E[u].emplace_back(v, e[i].w);
    }
    for (int t = 1; t <= cnt; t++) {
        int i = p[t];
        for (auto e : E[i]) {
            vis[e.first] = e.second;
        }
        for (auto e1 : E[i]) {
            int j = e1.first;
            for (auto e2 : E[j]) {
                int k = e2.first;
                if (vis[k]) {
                    int ij = e1.second, jk = e2.second, ki = vis[k];
                    ans += calc(i, j, k, ij, jk, ki) + calc(i, k, j, ki, jk, ij);
                    ans += calc(j, i, k, ij, ki, jk) + calc(j, k, i, jk, ki, ij);
                    ans += calc(k, i, j, ki, ij, jk) + calc(k, j, i, jk, ij, ki);
                }
            }
        }
        for (auto e : E[i]) {
            vis[e.first] = 0;
        }
    }
    printf("%u\n", ans.v);
    return 0;
}