「解题报告」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}(F(\frac{x}{p^c},\lfloor \frac{n}{p^i}\rfloor )-F(\frac{x}{p^{c-1}},\lfloor \frac{n}{p^{i+1}}\rfloor ))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;
}