自动搬运

查看原文

来自洛谷，原作者为

	ForgotMe 花与剑无痕，高挂一轮明灯。银蝶染旧梦，生死莫再问。

搬运于2025-08-24 22:28:20，当前版本为作者最后更新于2021-01-20 15:08:40，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

写个比较好懂的题解，看官方题解理解了半天。。。。

先写一下怎么推的柿子。（不想看的可以跳过）

$$\sum_{i=1}^n\sum_{j=1}^nf(\gcd(i,j))$$ $$\sum_{d=1}^nf(d)\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d] $$$$\sum_{d=1}^nf(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[\gcd(i,j)=1] $$$$\sum_{d=1}^nf(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[\gcd(i,j)=1] $$$$\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[\gcd(i,j)=1] $$

这一坨指数很经典，就是仪仗队。

所以答案为

$$\sum_{d=1}^nf(d)(2\sum_{i=1}^{\frac{n}{d}}\varphi(i)-1) $$

右边这一坨直接整除分块套杜教筛即可，考虑怎么处理 $f$ 函数的前缀和。

这里讲两个思路。

Way-1

$f(n) = \max(m-g(n)+1,0) n^k$ 中这坨 $\max(m-g(n)+1,0)$ 指数很恶心，让人无从下手，但发现 $g(n)$ 的最大值最多只为 $\log_2(n)$，于是考虑枚举 $g(n)$ 的值 $u$。

问题传化为怎么快速求

$$\sum_{i=1}^n[g(i)=u]i^k$$

这东西不好求，考虑转化。

$$\sum_{i=1}^n[g(i)=u]i^k=\sum_{i=1}^n[g(i)\le u]i^k-\sum_{i=1}^n[g(i)\lt u]i^k $$

这东西考虑用 min_25 直接求，简直就是 min_25 板子题，直接改一改 S 求和函数里的限制条件就行，然后差分直接乘上系数答案就出来了，理论时间复杂度 $\mathcal{O}(\sqrt{n}+kn+n^{\frac{3}{4}}+n^{\frac{2}{3}})$。

但是，这个东西我写了一下，跑 $2\times 10^9$ 本地都要跑 15S 左右，$10^{10}$ 直接 T 爆，反正我的代码是过不了的 （可能是我写的太丑了），有兴趣的同学可以去写一写，要是过了 @ 我一下。

不过这玩意不知道能不能用新版的 min_25 筛去搞，要是能行的话应该能过，反正我没试过，理论复杂度应该会变为 $\mathcal{O(n^{\frac{2}{3}}}\log n)$。

Way-2

还是不要搞着这些歪门邪道，要想正解，这题怎么可能这么简单。

同样的跟上面的思路一样，枚举 $g(n)$ 的值 $u$。

将 $f(n)$ 的系数搞一下事情，$f(n)$ 变为 $\max((m+1)-g(n),0)n^k$。

这启发我们可以转化问题，先算出所有的贡献，在减去不应该有的贡献。

形象化的表达即

$$\sum_{i=1}^n f(i)=(m+1)\sum_{i=1}^ni^k-\sum_{u=1}^{\min(m+1,\log_2(n))}\sum_{i=1}^n[u\le g(i)]i^k $$

前面的柿子就是自然数幂和，伯努利数，第二类斯特林数直接搞，推一推后面的柿子。

$$\sum_{u=1}^{\min(m+1,\log_2(n))}\sum_{i=1}^n[u\le g(i)]i^k $$

任意一个正整数的 $g(i)$ 都是大于等于 $1$（题目中特殊规定了 $g(1)=1$），所以可以拆一下柿子。

$$\sum_{u=2}^{\min(m+1,\log_2(n))}\sum_{i=1}^n[u\le g(i)]i^k+\sum_{i=1}^ni^k $$

然后来看前面这一坨柿子。

$$\sum_{u=2}^{\min(m+1,\log_2(n))}\sum_{i=1}^n[u\le g(i)]i^k $$

考虑做容斥，枚举一些质因子强制大于等于 $u$。

$$\sum_{u=2}^{\min(m+1,\log_2(n))}\sum_{i=1}^n[u\le g(i)]i^k=-\sum_{u=2}^{\min(m+1,\log_2(n))}\sum_{p=2}^{\sqrt[u]{n}}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{p^u}\rfloor}(ip^u)^k $$

然后移项。

$$\sum_{u=2}^{\min(m+1,\log_2(n))}\sum_{p=2}^{\sqrt[u]{n}}\mu(p)p^u\sum_{i=1}^{\lfloor\frac{n}{p^u}\rfloor}i^k $$

发现这东西可以直接算，后面那一坨就是自然数幂和，我们可以通过对 $n$ 整除分块得到对于任意一个 $i$，$\lfloor\dfrac{n}{i}\rfloor$ 的值，然后套上求自然数幂和的板子直接预处理出来，然后每次 $\mathcal{O(1)}$ 用就行了。

至于每次求前缀和的时间复杂度大约在 $\mathcal{O(\sqrt{n})}$，至于证明我不会，我可是个时间复杂度的白痴。

关于为啥这个式子是正确的。

$$\sum_{i=1}^n[u\le g(i)]i^k=-\sum_{p=2}^{\sqrt[u]{n}}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{p^u}\rfloor}(ip^u)^k $$

我来解释一下，因为莫比乌斯函数就相当于是个容斥，当 $p$ 有多个重复的质因子是值为 $0$，于是相当于枚举的不同质因子的乘积。

举个通俗易懂的栗子来说明一下这个式子，$\mu(2)=-1$，$\mu(3)=-1$，$\mu(6)=1$，当 $p=2$ 时，我减去了有质因子 $2$ 的次数大于等于 $u$ 的数的贡献，当 $p=3$ 时，我减去了有质因子 $3$ 的次数大于等于 $u$ 的数的贡献，但是你会发现我减多了，对于形如 $i6^u$ 的数我们会减去两次，于是需要加回来。

如果还不懂的话可以去搜一搜莫比乌斯函数与容斥，包看包懂。

代码

#include <cstdio>
#include <map>
#include <iostream>
#include <algorithm>
#include <queue>
#include <cstring>
#include <cmath>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
using namespace std;
using namespace __gnu_pbds;
#define LL long long
#define int long long
const int mod = 998244353;
const int inv2 = 499122177;
int dec(int x, int y) { return x >= y ? x - y : x + mod - y; }
int mul(int x, int y) { return 1ll * x * y % mod; }
int add(int x, int y) {
    if (x + y >= mod)
        return x + y - mod;
    return x + y;
}
int qkpow(LL a, int b) {
    int res = 1;
    for (; b; b >>= 1, a = mul(a, a))
        if (b & 1)
            res = mul(res, a);
    return res;
}
int sq, size, cnt;
int m, k, prime[10000005], id1[1000005], id2[1000005], Str[105][105], fk[10000005];
int h[1000005], hprime[2000005], phi[10000005], mn[10000005], mx[10000005], pw[10000005], f[10000005],
    inv[205];
int mu[200005], fs[200005];
LL n, num[1000005];
bool vis[10000005];

gp_hash_table<LL, int> Phi;

int Id(LL x) {
    if (x == 0)
        return 0;
    return (x <= sq ? id1[x] : id2[n / x]);
}

int S(LL n) {
    if (n == 1)
        return m;
    if (n <= 10000000)
        return f[n];
    if (fs[Id(n)])
        return fs[Id(n)];
    int maxg = log2(n) + 1;
    int tot = mul(m, h[Id(n)]);
    for (int p = 2; 1ll * p * p <= n; p++) {
        if (mu[p]) {
            int pkx = mul(pw[p], pw[p]);
            LL anopkx = 1ll * p * p;
            for (int u = 2; anopkx <= n && u <= min(m + 1, maxg);
                 u++, pkx = mul(pkx, pw[p]), anopkx = 1ll * anopkx * p) {
                int now = mul(mu[p], pkx);
                tot = add(tot, mul(now, h[Id(n / anopkx)]));
            }
        }
    }
    return fs[Id(n)] = tot;
}

int calcphi(LL n) {
    if (n <= 10000000)
        return phi[n];
    if (Phi[n])
        return Phi[n];
    int res = mul(n % mod, add(n, 1) % mod);
    res = mul(res, inv2);
    for (LL i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        res = dec(res, mul(dec(add(j % mod, 1) % mod, i % mod), calcphi(n / i)));
    }
    return Phi[n] = res;
}

signed main() {
    scanf("%lld %lld %lld", &n, &m, &k);
    for (int i = 1; i <= k + 2; i++) inv[i] = qkpow(i, mod - 2);
    Str[0][0] = 1;
    for (int i = 1; i <= k; i++)
        for (int j = 1; j <= k; j++) Str[i][j] = add(Str[i - 1][j - 1], mul(j, Str[i - 1][j]));
    sq = sqrt(n * 1.0);
    phi[1] = 1, f[1] = m, pw[1] = 1;
    for (int i = 2; i <= 10000000; i++) {
        if (!vis[i])
            prime[++cnt] = i, phi[i] = i - 1, pw[i] = qkpow(i, k), mn[i] = mx[i] = 1;
        for (int j = 1; j <= cnt && prime[j] * i <= 10000000; j++) {
            vis[i * prime[j]] = 1;
            pw[i * prime[j]] = mul(pw[i], pw[prime[j]]);
            if (i % prime[j] == 0) {
                mn[i * prime[j]] = mn[i] + 1, mx[i * prime[j]] = max(mx[i], mn[i] + 1);
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            mn[i * prime[j]] = 1, mx[i * prime[j]] = mx[i];
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
        f[i] = add(f[i - 1], mul(pw[i], max(m - mx[i] + 1, 0ll))) % mod;
        phi[i] = add(phi[i], phi[i - 1]);
    }
    for (int i = 1; i <= sq; i++) vis[i] = 0;
    cnt = 0;
    mu[1] = 1;
    for (int i = 2; i <= sq; i++) {
        if (!vis[i])
            mu[i] = -1, prime[++cnt] = i, fk[cnt] = pw[i], hprime[cnt] = add(hprime[cnt - 1], fk[cnt]);
        for (int j = 1; j <= cnt && prime[j] * i <= sq; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= sq; i++) mu[i] = (mu[i] == -1 ? mod - 1 : mu[i]);
    for (LL l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        num[++size] = n / r;
        LL Num = n / r % mod;
        int prod = 1;
        for (int i = 0; i <= k; ++i)
            prod = mul((Num + 1 - i) % mod, prod),
            h[size] = add(h[size], mul(mul(Str[k][i], prod), inv[i + 1]));
        if (n / r <= sq)
            id1[n / r] = size;
        else
            id2[r] = size;
    }
    int tot = 0, lstans = 0;
    for (LL l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        int phisum = dec(mul(2, calcphi(n / r)), 1), ano, fsum = S(r);
        ano = fsum;
        fsum = dec(fsum, lstans);
        lstans = ano;
        tot = add(tot, mul(fsum, phisum));
    }
    printf("%lld\n", tot);
    return 0;
}