题解 SP20173 DIVCNT2
铃悬
2021-08-09 16:40:30
给定 $n$,求]
$$ \sum_{i = 1}^n \sigma_0(i^2). $$
我们发现 $\sigma_0(n^2)$ 这个函数等于 $\mu^2(n)$ 和 $\sigma_0(n)$ 的 Dirichlet 卷积(都是积性函数,观察素数幂处的值就可以发现这一点).
因此利用整除分块,只需要计算
$\mu^2$ 以及 $\sigma_0$ 在每个 $\lfloor n / i \rfloor$ 处的值.
考虑这两者怎么计算:
对于 $\mu^2$,我们有
$$\sum_{i = 1}^n \mu^2(i)
= \sum_{k = 1}^{\lfloor \sqrt{n} \rfloor}
\mu(i) \Bigl\lfloor \frac{n}{d^2} \Bigr\rfloor$$
可以利用类似整除分块的办法($\lfloor n / d^2 \rfloor$ 只有 $O(n^{1/3})$ 种取值)预处理 $\mu$ 在 $1, \dots, \sqrt{n}$ 的前缀和后在 $O(n^{1/3})$ 的时间内算出 $\mu^2$ 的一个前缀和.
对于 $\sigma_0$,用折线拟合双曲线的办法据说是 $O(n^{1/3}\log n)$ 的,
参见 [DIVCNT1 的题解](https://www.luogu.com.cn/problem/solution/SP26073)
或者 [yhx 的题解](https://yhx-12243.github.io/OI-transit/records/spojDIVCNT1.html)(事实上我的代码里这一部分就是抄的x)
如果我们线性预处理出 $\mu^2$ 和 $\sigma_0$ 在 $1, \dots, B$ 的前缀和,而对那些大于 $B$ 的 $n / i$ 都 $O(n^{1/3} \log n)$ 暴力计算,时间复杂度是
$$B + \sum_{i = 1}^{n / B} \Bigl( \frac{n}{i} \Bigr)^{1/3} \log n
= O\Bigl(B + \frac{n\log n}{B^{2/3}}\Bigr)$$
取
$B = (Tn \log n)^{3/5}$ 就可以得到 $O((Tn \log n)^{3/5})$ 的算法了.
```cpp
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
typedef long long LL;
const int B = 20000000;
namespace sum_sigma {
int S[B];
struct P {
LL x, y;
P(LL x = 0, LL y = 0) : x(x), y(y) {}
inline P operator+(const P &a) const { return P(x + a.x, y + a.y); }
} stack[100000];
LL N;
LL solve(LL N) {
if (N < B) return S[N];
sum_sigma::N = N;
LL n13 = std::cbrt(N), n12 = std::sqrt(N);
LL x = N / n12, y = n12 + 1, ans = 0;
int top = 0;
stack[++top] = P(1, 0);
stack[++top] = P(1, 1);
P L, R;
while (true) {
for (L = stack[top--]; (x + L.x) * (y - L.y) > N; x += L.x, y -= L.y)
ans += x * L.y + (L.y + 1) * (L.x - 1) / 2;
if (y <= n13) break;
for (R = stack[top]; (x + R.x) * (y - R.y) <= N; R = stack[--top]) L = R;
while (true) {
P M = L + R;
if ((x + M.x) * (y - M.y) > N) stack[++top] = R = M;
else {
if (N * R.x <= (x + M.x) * (x + M.x) * R.y) break;
L = M;
}
}
}
for (int i = 1; i < y; ++i) ans += N / i;
return ans * 2 - n12 * n12;
}
};
namespace sum_mu {
int mu[B], S[B];
LL solve(LL N) {
if (N < B) return S[N];
LL ans = 0;
int i = 1, j;
while ((LL)i * i <= N) {
j = std::sqrt(N / (N / i / i));
ans += (mu[j] - mu[i - 1]) * (N / i / i);
i = j + 1;
}
return ans;
}
};
void Init() {
int *Ss = sum_sigma::S, *mu = sum_mu::mu, *Su = sum_mu::S;
static bool vis[B];
static int pr[650000], q[B], cnt;
Ss[1] = mu[1] = 1;
for (int i = 2; i < B; ++i) {
if (!vis[i]) {
Ss[i] = 2; mu[i] = -1; q[i] = 1;
pr[cnt++] = i;
}
for (int j = 0, t = (B - 1) / i; pr[j] <= t; ++j) {
int l = i * pr[j]; vis[l] = true;
if (i % pr[j] == 0) {
mu[l] = 0;
Ss[l] = Ss[i] + q[i];
q[l] = q[i];
break;
} else {
mu[l] = -mu[i];
Ss[l] = 2 * Ss[i];
q[l] = Ss[i];
}
}
}
for (int i = 2; i < B; ++i) Ss[i] += Ss[i - 1];
for (int i = 1; i < B; ++i) {
Su[i] = Su[i - 1] + (mu[i] != 0);
mu[i] += mu[i - 1];
}
}
int main() {
Init();
int T; scanf("%d", &T);
while (T--) {
LL n; scanf("%lld", &n);
LL ans = 0, ss = 0;
LL i = 1, j;
while (i <= n) {
j = n / (n / i);
LL t = sum_mu::solve(j);
ans += (t - ss) * sum_sigma::solve(n / i);
ss = t;
i = j + 1;
}
printf("%lld\n", ans);
}
}
```