P6384 『MdOI R2』Quo Vadis

更好的阅读体验

题目大意

给定一个 无限大 的矩阵 $A$，其中 $A_{i,j}=ij\gcd(i,j)$。

接下来有 $m$ 个操作，每行 $1$ 至 $3$ 个整数，意义如下：

• $1$：对矩阵 $A$ 进行高斯消元，使之成为一个上三角矩阵。（注意：这里，高斯消元中只允许将一行的某一个倍数加到另一行上，不允许交换任意两行，不允许将某行乘上一个倍数，保证这样之后仍然可以得到上三角矩阵，并且保证消元之后的矩阵的每个元素均为非负整数。）
• $2\ x\ y$：求出当前矩阵的 $A_{x,y}$。
• $3\ x$：求出 $\sum\limits_{i=1}^{x}\sum\limits_{j=1}^{x}A_{i,j}$。
• $4\ x$：设 $B$ 是一个 $x$ 阶矩阵，其中 $B_{i,j}=A_{i,j}$，你需要求出 $\det B$。

上述所有答案对 $998244353$ 取模。

解题思路

首先，对于矩阵 $A$ 满足 $A_{i,j}=i j \gcd(i,j)$，在进行高斯消元后，有：

$$A'_{i,j}=\begin{cases}i j \varphi(i) && i \mid j \\0 && i \nmid j\end{cases}$$

对于矩阵 $B$，满足 $B_{i,j}=\gcd(i,j)$，在进行高斯消元后，有：

$$B'_{i,j}=\begin{cases}\varphi(i) && i \mid j\\0 && i \nmid j\end{cases}$$

那么，对于 $1$ 操作前的 $2$ 操作答案为 $i j \gcd(i,j)$，对于 $1$ 操作后的 $2$ 操作答案为 $i j \varphi(i)$。

再看 $1$ 操作前的 $3$ 操作，答案为：

$$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}i j \gcd(i,j)\\ =\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}i j [\gcd(i,j=1)]\\ =\sum\limits_{d=1}^{n}d^3\sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{n}{d} \rfloor}i j [\gcd(i,j)=1]\\ =\sum\limits_{d=1}^{n}d^3 \sum\limits_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{n}{d}\rfloor}i j \sum\limits_{p|\gcd(i,j)}\mu(p)\\ =\sum\limits_{d=1}^{n}d^3\sum\limits_{p=1}^{\lfloor\frac{n}{d}\rfloor} \mu(p)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}i j\\ =\sum\limits_{d=1}^{n}d^3\sum\limits_{p=1}^{\lfloor\frac{n}{d}\rfloor}p^2 \mu(p)\sum\limits_{i=1}^{\lfloor\frac{n}{pd}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{pd}\rfloor}i j\\ =\sum\limits_{d=1}^{n}d^3\sum\limits_{p=1}^{\lfloor\frac{n}{d}\rfloor}p^2 \mu(p)\left(\sum\limits_{i=1}^{\lfloor\frac{n}{pd}\rfloor}i\right)^2\\ \sum\limits_{d=1}^{n}d\sum\limits_{p=1}^{\lfloor\frac{n}{d}\rfloor}(pd)^2 \mu(p)\left(\sum\limits_{i=1}^{\lfloor\frac{n}{pd}\rfloor}i\right)^2$$

设 $f(n)=\left(\sum\limits_{i=1}^{n}i\right)^2$，则有：

$$\sum\limits_{d=1}^{n}d\sum\limits_{p=1}^{\lfloor\frac{n}{d}\rfloor}(pd)^2 \mu(p)f(\lfloor\frac{n}{pd}\rfloor)$$

设 $T=pd$，有：

$$\sum\limits_{T}^{n}T^2f(\lfloor\frac{n}{T}\rfloor)\sum\limits_{d|T}d \mu(\frac{T}{d})$$

根据狄利克雷卷积，有：

$$\sum\limits_{T}^{n}T^2f(\lfloor\frac{n}{T}\rfloor)\varphi(T)$$

数论分块套杜教筛即可。

对于 $1$ 操作后的 $3$ 操作，设 $g(n)$ 为第 $n$ 列的和，那么有：

$$g(n)=n\sum\limits_{d|n}d\varphi(d)$$

设 $h(n)=\sum\limits_{d|n}d \varphi(d)$，显然， $h(n)$ 是个积性函数，线性筛即可。

答案即为：

$$\sum\limits_{i=1}^{n}g(i)\\ =\sum\limits_{i=1}^{n}\left(i \times h(i)\right)$$

对于操作 $4$，行列式即消元后主对角线上元素乘积，所以答案即为：

$$\prod\limits_{i=1}^{n}i^2\varphi(i)$$

CODE

#include <bits/stdc++.h>

using namespace std;

#define int long long

inline int read()
{
    int x = 0, f = 1;
    char c = getchar();
    while (c < '0' || c > '9')
    {
        if (c == '-')
            f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9')
    {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * f;
}

const int MOD = 998244353, _ = 5e6 + 51, INV6 = 166374059, INV4 = 748683265, MAXM = 4e5 + 7;

unordered_map<int, int> resP;

int n, op, flag, tot;

int x, y;

int prime[MAXM], phi[_], vis[_], p[_], low[_], g[_];

int invPr[MAXM], prefix[_], prod[_];

inline int qpow(int x, int y)
{
    int res = 1;
    while (y)
    {
        if (y & 1)
            res = res * x % MOD;
        x = x * x % MOD;
        y >>= 1;
    }
    return res;
}

inline void solve(int n)
{
    int pp;
    phi[1] = g[1] = prod[0] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!vis[i])
        {
            low[i] = prime[++tot] = i, phi[i] = i - 1, g[i] = (i * i % MOD - i + 1 + MOD) % MOD;
            invPr[tot] = qpow(prime[tot] + 1, MOD - 2);
        }
        for (int j = 1; j <= tot && i * prime[j] <= n; j++)
        {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0)
            {
                low[i * prime[j]] = low[i] * prime[j];
                phi[i * prime[j]] = phi[i] * prime[j];
                if (low[i] == i)
                {
                    pp = i * prime[j];
                    g[pp] = (pp * pp % MOD * prime[j] % MOD + 1) * qpow(prime[j] + 1, MOD - 2) % MOD;
                }
                else
                {
                    g[i * prime[j]] = g[i / low[i]] * g[low[i] * prime[j]] % MOD;
                }
                break;
            }
            low[i * prime[j]] = prime[j], g[i * prime[j]] = g[i] * g[prime[j]] % MOD;
            phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
    for (int i = 1; i <= n; i++)
    {
        p[i] = (p[i - 1] + i * i % MOD * phi[i] % MOD) % MOD;
        prefix[i] = (prefix[i - 1] + i * g[i] % MOD) % MOD;
        prod[i] = prod[i - 1] * i % MOD * i % MOD * phi[i] % MOD;
    }
}


inline int get_phi(int x)
{
    int res = x;
    for (int i = 2; i * i <= x; i++)
        if (x % i == 0)
        {
            res = res / i * (i - 1);
            while (x % i == 0)
                x /= i;
        }
    if (x != 1)
        res = res / x * (x - 1);
    return res;
}

inline int calc2(int x)
{
    return x * (x + 1) % MOD * (2 * x + 1) % MOD * INV6 % MOD;
}

inline int calc3(int x)
{
    return x * (x + 1) % MOD * x % MOD * (x + 1) % MOD * INV4 % MOD;
}

inline int prefixP(int num)
{
    if (num <= 5e6)
        return p[num];
    if (resP[num])
        return resP[num];
    int res = calc3(num);
    for (int l = 2, r; l <= num; l = r + 1)
    {
        r = num / (num / l);
        res = (res - prefixP(num / l) * (calc2(r) - calc2(l - 1) + MOD) % MOD + MOD) % MOD;
    }
    return resP[num] = res;
}

inline int calc(int num)
{
    int res = 0;
    for (int l = 1, r; l <= num; l = r + 1)
    {
        r = num / (num / l);
        res = (res + calc3(num / l) * (prefixP(r) - prefixP(l - 1) + MOD) % MOD) % MOD;
    }
    return res;
}

signed main()
{
    solve(5000010);
    n = read();
    while(n--)
    {
        op = read();
        if (op == 1)
            flag = 1;
        if (op == 2)
        {
            x = read(), y = read();
            if (flag)
                printf("%lld\n", y % x ? 0 : x * y % MOD * get_phi(x) % MOD);
            else
                printf("%lld\n", (x % MOD) * (y % MOD) % MOD * (__gcd(x, y) % MOD) % MOD);
        }
        if (op == 3)
        {
            x = read();
            if (flag)
                printf("%lld\n", prefix[x]);
            else
                printf("%lld\n", calc(x));
        }
        if (op == 4)
        {
            x = read();
            printf("%lld\n", prod[x]);
        }
    }
    return 0;
}