首先显然地，有答案为 \[ \sum\limits_{d=1}^n f(d) \left(2\sum\limits_{i=1}^{\left\lfloor\frac nd\right\rfloor} \varphi(i) - 1\right) \]
杜教筛 / Min_25 筛一下，那么问题变为求 \(f(i)\) 的前缀和。
设 \(s(n) = \sum\limits_{i=1}^n f(i)\)。
观察到 \(f(n)\) 实际上并不好处理，考虑把那个 \(\max(m-g(n)+1,0)\) 的系数处理一下： \[ f(n) = \sum\limits_{u=1}^m [g(n) \le u] n^k \]
但是 \(m\) 可能很大。
然而注意到 \(g(n) \le \log_2 n\)，所以可以把 \(m\) 降到 \(O(\log n)\) 级别，超出部分直接算一下自然数幂和即可。
则 \[ \begin{aligned} s(n) &= \sum\limits_{i=1}^n f(i) \\ &= \sum\limits_{i=1}^n \sum\limits_{u=1}^m [g(i) \le u] i^k \\ &= \sum\limits_{u=1}^m \sum\limits_{i=1}^n [g(i) \le u] i^k \end{aligned} \]
考虑求 \(\sum\limits_{i=1}^n [g(i) \le u]i^k\)。
这可以考虑容斥，即枚举某些质因子，使它们的指数强制超过 \(u\)。
则 \[
\begin{aligned}
\sum\limits_{i=1}^n [g(i) \le u]i^k
&= \sum\limits_{p=1}^{\lfloor n^{1/(u+1)}\rfloor} \mu(p) \sum\limits_{i=1}^{\left\lfloor\frac n{p^{u+1}}\right\rfloor} (ip^{u+1})^k \\
&= \sum\limits_{p=1}^{\lfloor n^{1/(u+1)}\rfloor} \mu(p) p^{(u+1)k} \sum\limits_{i=1}^{\left\lfloor\frac n{p^{u+1}}\right\rfloor} i^k
\end{aligned}
\]
然后注意到 \(\sum\limits_{u=1}^m O(n^{1/(u+1)}) = O(\sqrt n)\)，所以对于 \(u=1,\dots, m\) 如此计算一次的复杂度是 \(O(\sqrt n)\) 的。
另外可以通过线筛 \(g(n)\) 预处理出一定范围内的 \(s(n)\)。
于是根据杜教筛的思想，平衡至 \(O(n^{2/3})\) 即可。
不过还有一个问题，那就是计算自然数 \(k\) 次幂和。
这是一个经典问题，容易发现需要计算的位置只有 \(O(\sqrt n)\) 种，有很多方法可以做到 \(O(k^2)\) 预处理后 \(O(k\sqrt n)\) 解决。
标程使用的是第二类斯特林数。
总复杂度 \(O(k^2 + k\sqrt n + n^{2/3})\)。
可能略微卡常。
代码： 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
using namespace std;
const long long N = 1e10;
const int MX = 5e6;
const int LIM = 1e5;
const int M = 33;
const int K = 100;
const int mod = 998244353;
inline int fpow(int a,int b)
{
int ret = 1;
for(;b;b >>= 1)
(b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
return ret;
}
long long n;
int lim,m,k;
int inv[K + 5],S[K + 5][K + 5];
int vis[MX + 5],cnt,prime[MX + 5];
int mn[MX + 5],mx[MX + 5],f[MX + 5],phi[MX + 5],mu[MX + 5],pw[MX + 5];
int le[LIM + 5],ge[LIM + 5],tot;
int sk[2 * LIM + 5];
inline int &id(long long x)
{
return x <= lim ? le[x] : ge[n / x];
}
int w[2][MX + 5];
inline int &mem_sphi(long long x)
{
return w[0][n / x];
}
inline int &mem_sf(long long x)
{
return w[1][n / x];
}
int sphi(long long n)
{
if(n <= MX)
return phi[n];
if(~mem_sphi(n))
return mem_sphi(n);
int ret = n & 1 ? ((n + 1 >> 1) % mod) * (n % mod) % mod : ((n >> 1) % mod) * ((n + 1) % mod) % mod;
for(register long long l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret = (ret - (r - l + 1) % mod * sphi(n / l) % mod + mod) % mod;
}
return mem_sphi(n) = ret;
}
int sf(long long n)
{
if(n <= MX)
return f[n];
if(~mem_sf(n))
return mem_sf(n);
int ret = m > M ? (long long)(m - M) * sk[id(n)] % mod : 0;
for(register int i = 1;(long long)i * i <= n;++i)
if(mu[i])
{
long long j = n / i / i;
int s = (long long)pw[i] * pw[i] % mod;
for(register int u = 1;u <= min(m,M) && j;++u,j /= i,s = (long long)s * pw[i] % mod)
ret = (ret + (long long)s * mu[i] % mod * sk[id(j)]) % mod;
}
return mem_sf(n) = ret;
}
int ans;
int main()
{
memset(w,-1,sizeof w);
scanf("%lld%d%d",&n,&m,&k),lim = sqrt(n);
inv[1] = 1;
for(register int i = 2;i <= k + 1;++i)
inv[i] = (long long)(mod - mod / i) * inv[mod % i] % mod;
S[0][0] = 1;
for(register int i = 1;i <= k;++i)
for(register int j = 1;j <= k;++j)
S[i][j] = (S[i - 1][j - 1] + (long long)j * S[i - 1][j]) % mod;
phi[1] = mu[1] = pw[1] = 1,f[1] = m;
for(register int i = 2;i <= MX;++i)
{
if(!vis[i])
mn[prime[++cnt] = i] = mx[i] = 1,phi[i] = i - 1,mu[i] = mod - 1,pw[i] = fpow(i,k);
for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
{
vis[i * prime[j]] = 1,pw[i * prime[j]] = (long long)pw[i] * pw[prime[j]] % mod;
if(!(i % prime[j]))
{
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),mu[i * prime[j]] = (mod - mu[i]) % mod;
}
f[i] = (f[i - 1] + (long long)pw[i] * max(m - mx[i] + 1,0)) % mod,
phi[i] = (phi[i] + phi[i - 1]) % mod;
}
for(register long long l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
id(n / l) = ++tot;
int prod = 1;
for(register int i = 0;i <= k;++i)
prod = (n / l + 1 - i) % mod * prod % mod,
sk[tot] = (sk[tot] + (long long)S[k][i] * prod % mod * inv[i + 1]) % mod;
}
for(register long long l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
ans = (ans + (2LL * sphi(n / l) - 1 + mod) * (sf(r) - sf(l - 1) + mod)) % mod;
}
printf("%d\n",ans);
}