Luogu9157「GLR-R4」夏至
抢到最优解了，UOJ 校验码上 80pts 过不去。/kk
这里是官方题解的简化。
首先考虑 \(n=1\) 怎么做，相当于对 \(m\le 10^{10}\) 筛出 \(f\) 的前缀和。由于 \(f(p)=p\)，直接构造函数 \(g(n)=n\) 然后 PN 筛 \(O(\sqrt m)\) 求即可。
然后考虑 \(n>1\)，由于 \(n\) 比较小，考虑对每一个 \(i\le n\)，求出 \( \sum\limits_{j=1}^mf(ij)\)。
令 \(F(x,y)=\sum\limits_{i=1}^yf(ix)\)，答案就是 \(\sum\limits_{i=1}^nF(i,m)\)。类似 Min25 筛的想法，消去 \(x\) 的最大质因子 \(p^c\)，然后枚举 \(y'\le y\) 中 \(y'\) 中的 \(p\) 的最大次数为 \(i\)，用容斥即 「至少出现了 \(i\) 次的贡献」减去「至少出现了 \(i+1\) 次的贡献」，由于钦定至少出现 \(i+1\) 次比至少出现 \(i\) 次少了一个 \(p\) 的贡献，将 \(p\) 乘到 \(\frac{x}{p^c}\) 上面：
\[F(x,y)=\sum\limits_{i\ge 0}\left(F\left(\frac{x}{p^c},\left\lfloor \frac{y}{p^i}\right\rfloor\right)-F\left(\frac{x}{p^{c-1}},\left\lfloor\frac{y}{p^{i+1}}\right\rfloor\right)\right)f(p^{c+i})
\]
记忆化一下，边界 \(n=1\) 直接 PN 筛，否则暴力递归即可。预处理 \(xy\le 10^6\) 的所有 \(F(x,y)\) 后，根据官方题解有效状态数大概只有 \(10^6\) 左右，转移数大概 \(1.5\times 10^7\)。
// Problem: P9157 「GLR-R4」夏至
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P9157
// Memory Limit: 512 MB
// Time Limit: 4000 ms
//
// Powered by CP Editor (https://cpeditor.org)
#include <bits/stdc++.h>
#include <bits/extc++.h>
#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2")
#define pb emplace_back
#define mt make_tuple
#define mp make_pair
#define fi first
#define se second
using namespace std;
using namespace __gnu_cxx;
using namespace __gnu_pbds;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<ll, ll> pi;
typedef tuple<int, int, int> tu;
bool Med;
const int C = 40;
const int N = 1e6 + 100;
const int M = 1e5 + 100;
const int P = 1e9 + 7;
const int i2 = (P + 1) / 2;
ll m;
int n, k, h[M][C];
int tot, vs[N], pr[N], mxp[N], mxc[N], rs[N], f[N];
vector<int> F[N];
vector<ull> pn;
gp_hash_table<ull, int> tF;
void Add(int &x, int y) { x += y, (x >= P) && (x -= P); }
int qpow(int p, int q) {
int res = 1;
for (; q; q >>= 1, p = 1ll * p * p % P)
if (q & 1) res = 1ll * res * p % P;
return res;
}
void dfs(ll x, int y, int i) {
if (!y || i == tot + 1) return;
pn.pb((ull)x * P + y);
for (int j = i; j <= tot; j++) {
if (x > m / pr[j] / pr[j]) break;
ll t = x * pr[j] * pr[j];
for (int l = 2; t <= m; t *= pr[j], l++)
if (h[j][l]) dfs(t, 1ll * y * h[j][l] % P, j + 1);
}
}
void init(int lim) {
for (int i = 2; i <= lim; i++) {
if (!vs[i]) pr[++tot] = i, mxp[i] = i, mxc[i] = rs[i] = 1;
for (int j = 1; j <= tot && i * pr[j] <= lim; j++) {
vs[i * pr[j]] = 1;
mxp[i * pr[j]] = mxp[i];
mxc[i * pr[j]] = mxc[i] + (mxp[i] == pr[j]);
rs[i * pr[j]] = (pr[j] == mxp[i]) ? rs[i] : (rs[i] * pr[j]);
if (i % pr[j] == 0) break;
}
}
f[1] = 1;
for (int i = 2; i <= lim; i++)
f[i] = f[rs[i]] * qpow(mxp[i], __gcd(k, mxc[i]));
for (int i = 1; i <= lim; i++) {
F[i].resize(lim / i + 5);
for (int j = 1; j <= lim / i; j++)
Add(F[i][j] = F[i][j - 1], f[i * j]);
}
for (int i = 1; i <= tot; i++) {
h[i][0] = 1; ll tp = pr[i];
for (int j = 1; tp <= m; tp *= pr[i], j++) {
h[i][j] = qpow(pr[i], __gcd(j, k));
for (int l = 1, t = pr[i]; l <= j; l++, t = 1ll * t * pr[i] % P)
Add(h[i][j], P - 1ll * t * h[i][j - l] % P);
}
}
}
int S(ll x) { x %= P; return 1ll * x * (x + 1) % P * i2 % P; }
int calc(int x, ll y) {
if (x * y <= 1e6) return F[x][y];
if (!y) return 0;
if (tF.find((ull)x * P + y) != tF.end()) return tF[(ull)x * P + y];
int res = 0;
if (x == 1) {
for (ull p : pn) {
if (p / P > y) break;
Add(res, 1ll * (p % P) * S(y / (p / P)) % P);
}
} else {
int p = mxp[x]; ll t = y;
for (int i = 0; t; i++, t /= p)
Add(res, 1ll * (calc(rs[x], t) - calc(rs[x] * p, t / p) + P) % P * qpow(p, __gcd(i + mxc[x], k)) % P);
}
return tF[(ull)x * P + y] = res;
}
void solve() {
cin >> n >> m >> k, init(1e6);
dfs(1, 1, 1), sort(pn.begin(), pn.end(), [](ull &x, ull &y) { return x / P < y / P; });
int res = 0;
for (int i = 1; i <= n; i++)
Add(res, calc(i, m));
cout << res << '\n';
}
bool Mbe;
signed main() {
// ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
cerr << (&Med - &Mbe) / 1048576.0 << " MB\n";
#ifdef FILE
freopen("1.in", "r", stdin);
freopen("1.out", "w", stdout);
#endif
int T = 1;
// cin >> T;
while (T--) solve();
cerr << (int)(1e3 * clock() / CLOCKS_PER_SEC) << " ms\n";
return 0;
}
