题解 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 的题解 或者 yhx 的题解(事实上我的代码里这一部分就是抄的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})$ 的算法了.

#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);
  }
}