「解题报告」P9157「GLR-R4」夏至

开学了，所以我要颓废了。

这题感觉好强，做法也好神奇。

首先我们考虑 \(n=1\) 时咋做，那么就是一个积性函数前缀和。

发现这个函数满足 \(f(p) = p\)，那么考虑上 PN 筛，可以 \(O(\sqrt{m})\) 内求得答案。

然后考虑 \(n > 1\)。\(n\) 很小，所以我们可以考虑对于每一个 \(i\)，求出 \(\sum_{j=1}^m f(ij)\)。

那么设 \(F(x, n) = \sum_{i=1}^n f(ix)\)。考虑类似于 min_25 筛的思想，每次消去一个最大质因子，然后枚举 \(i\) 中 \(p\) 的次数 \(c\)。考虑容斥求恰好 \(p\) 的次数为 \(c\) 的贡献，拿至少 \(c\) 次减去至少 \(c-1\) 次，那么就可以得到：

\[F(x, n) = \sum_{i\ge 0} \left(F\left(\frac{x}{p^c}, \left\lfloor\frac{n}{p^i}\right\rfloor\right) - F\left(\frac{x}{p^{c - 1}}, \left\lfloor\frac{n}{p^{i + 1}}\right\rfloor\right)\right) f(p^{c + i}) \]

其中 \(c\) 为 \(x\) 中 \(p\) 的次数。\(\frac{x}{p^{c - 1}}\) 的原因是你钦定了还有至少一个 \(p\) 的因子，所以给 \(x\) 乘了 \(p\)。

记忆化一下，暴力跑这个东西即可。复杂度不会证，但是有效状态数大概 \(2 \times 10^6\)，转移数大概 \(3 \times 10^7\)。

然后直接跑可能会 T 一个点，可以预处理出 \(x \cdot n \le 10^6\) 的 \(F(x, n)\)，然后就能跑过去了。

这个做法很神奇啊，基本上没咋依赖 \(f(x)\) 的性质就能做了。

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1000005, P = 1000000007, N = 1000000;
int n, k;
long long m;
struct Hash {
    size_t operator()(const pair<int, long long> &p) const {
        return ((p.first * 13331 + p.second) ^ (p.second >> 16 | ((long long) p.first) << 12)) % P;
    }
};
unordered_map<pair<int, long long>, int, Hash> f;
int pri[MAXN], pcnt, mxp[MAXN], lp[MAXN], pc[MAXN];
bool vis[MAXN];
int qpow(int a, int b) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = 1ll * ans * a % P;
        a = 1ll * a * a % P;
        b >>= 1;
    }
    return ans;
}
int h[MAXN][88];
vector<pair<long long, int>> pn;
void dfs(long long x, int hx, int i) {
    if (!hx) return;
    pn.push_back({x, hx});
    for (int j = i; j <= pcnt; j++) {
        if (x > m / pri[j] / pri[j]) break;
        long long t = x * pri[j] * pri[j];
        for (int c = 2; t <= m; c++, t *= pri[j]) {
            dfs(t, 1ll * hx * h[j][c] % P, j + 1);
        }
    }
}
int G(long long n) {
    n %= P;
    return 1ll * (n + 1) * n % P * ((P + 1) / 2) % P;
}
vector<int> ff[MAXN];
int fx[MAXN];
int F(int x, long long n) {
    if (x * n <= N) {
        return ff[x][n];
    }
    if (!n) return 0;
    if (f.count({x, n})) return f[{x, n}];
    if (x == 1) {
        // PN 筛
        int ans = 0;
        for (auto p : pn) {
            if (p.first > n) break;
            ans = (ans + 1ll * p.second * G(n / p.first)) % P;
        }
        return f[{x, n}] = ans;
    }
    int p = mxp[x], x2 = lp[x];
    long long m = n;
    int ans = 0;
    for (int i = 0; m; i++, m /= p) {
        ans = (ans + 1ll * (F(x2, m) - F(x2 * p, m / p) + P) * qpow(p, __gcd(i + pc[x], k))) % P;
    }
    return f[{x, n}] = ans;
}
int main() {
    scanf("%d%lld%d", &n, &m, &k);
    fx[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!vis[i]) {
            pri[++pcnt] = i;
            mxp[i] = i;
            lp[i] = 1;
            pc[i] = 1;
        }
        for (int j = 1; j <= pcnt && i * pri[j] <= N; j++) {
            vis[i * pri[j]] = 1;
            mxp[i * pri[j]] = mxp[i];
            lp[i * pri[j]] = pri[j] == mxp[i] ? lp[i] : lp[i] * pri[j];
            pc[i * pri[j]] = pc[i] + (mxp[i] == pri[j]);
            if (i % pri[j] == 0) {
                break;
            }
        }
        fx[i] = fx[lp[i]] * qpow(mxp[i], __gcd(pc[i], k));
    }
    ff[0] = vector<int>(N + 1, 0);
    for (int i = 1; i <= N; i++) {
        ff[i].push_back(0);
        for (int j = 1; i * j <= N; j++) {
            ff[i].push_back((ff[i][j - 1] + fx[i * j]) % P);
        }
    }
    for (int i = 1; i <= pcnt; i++) {
        h[i][0] = 1;
        long long x = pri[i];
        for (int j = 1; x <= m; j++, x *= pri[i]) {
            h[i][j] = qpow(pri[i], __gcd(j, k));
            int y = pri[i];
            for (int k = 1; k <= j; k++, y = 1ll * y * pri[i] % P) {
                h[i][j] = (h[i][j] - 1ll * h[i][j - k] * y % P + P) % P;
            }
        }
    }
    dfs(1, 1, 1);
    sort(pn.begin(), pn.end());
    int ans = 0;
    for (int i = 1; i <= n; i++) {
        ans = (ans + F(i, m)) % P;
    }
    printf("%d\n", ans);
    return 0;
}