题解 P7277 【平凡点滴】
ForgotMe
·
·
题解
写个比较好懂的题解，看官方题解理解了半天。。。。
先写一下怎么推的柿子。（不想看的可以跳过）
\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
问题传化为怎么快速求
$$
\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$ 的数我们会减去两次，于是需要加回来。
如果还不懂的话可以去搜一搜莫比乌斯函数与容斥，~~包看包懂~~。
代码
```cpp
#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;
}
```