题解 SP20173 DIVCNT2

铃悬

2021-08-09 16:40:30

Solution

给定 $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); } } ```