Counting Divisors (Square)

题目描述

Let $\sigma_0(n)$ be the number of positive divisors of $n$.

For example, $\sigma_0(1)=1, \ \sigma_0(2)=2$ and $\sigma_0(6)=4$.

Let $S_2(n)=\sum_{i=1}^n \sigma_0(i^2)$.

Given $N$, find $S_2(N)$.

题意概述

令$\sigma_0(n)$表示$n$的正因数个数，$S_2(n)=\sum_{i=1}^n \sigma_0(i^2)$。给定$N$，求$S_2(N)$。

数据范围：$1 \le N \le 10^{12}$。

算法分析

为了方便，我们令$f(n)=\sigma_0(n^2)$，$\omega(n)$表示$n$的不同质因子个数。

设$n=\prod_{k=1}^{\omega(n)} p_k^{a_k}$，其中$p_k$表示不同质数。

先来推$f(n)$。对于$n$的每个因数$d$，将$d$的某些质因子$p_k$的次数加上$a_k$得到$d’$，有$d’ \nmid n, \ d’ \mid n^2$。这样的方案数有$(2^{\omega(d)}-1)$种，再算上$d$本身，可得

$$
\begin{align}
f(n)=\sum_{d \mid n} 2^{\omega(d)}
\end{align}
$$

化简$2^{\omega(d)}$。它等价于$d$的每个质因数要么选一个，要么不选。那么它就等于$d$的因数中 square-free 的数的个数。因此

$$
\begin{align}
2^{\omega(d)}=\sum_{i \mid d} \mu^2(i)
\end{align}
$$

回到原式

$$
\begin{align}
S_2(n)=\sum_{i=1}^n \sum_{j \mid i} \sum_{k \mid j} \mu^2(k)
\end{align}
$$

还是毫无头绪。我们继续化简$f(n)$。考虑它的卷积形式$f=\mu^2 \ast 1 \ast 1=\mu^2 \ast \sigma_0$。所以

$$
\begin{align}
S_2(n)&=\sum_{i=1}^n \sum_{j \mid i} \mu^2(j) \sigma_0\left({i \over j}\right) \\
&=\sum_{ij \le n} \mu^2(i) \sigma_0(j)
\end{align}
$$

那么只要用反演求出$\mu^2$的前缀和、用分块求出$\sigma_0$的前缀和，就可以利用分块的技巧了。具体如下

$$
\begin{align}
\sum_{i=1}^n \mu^2(i)&=\sum_{i=1}^n [i\text{最大平方因子}=1] \\
&=\sum_{i=1}^n \sum_{j^2 \mid i} \mu(j) \\
&=\sum_{i=1}^{\lfloor \sqrt{n} \rfloor} \mu(i) \left\lfloor {n \over i^2} \right\rfloor \\
\sum_{i=1}^n \sigma_0(i)&=\sum_{i=1}^n \left\lfloor {n \over i} \right\rfloor
\end{align}
$$

代码

#include <cmath>
#include <cstdio>
#include <cstring>

void read(long long &t) {
  char c; while ((c = getchar()) < '0' || c > '9') ; t = c - '0';
  while ((c = getchar()) >= '0' && c <= '9') (t *= 10) += c - '0';
}

int N;
char mui[100000000];
int d[100000000], sum[100000000], prime[10000000];
bool vis[100000000];
long long n, sigma[100000000];

void init() {
  int top = 0; mui[1] = 1, sigma[1] = 1;
  for (int i = 2; i < N; ++ i) {
    if (! vis[i]) prime[top ++] = i, mui[i] = -1, sigma[i] = 2, d[i] = 1;
    for (int j = 0; j < top; ++ j) {
      int k = i * prime[j]; if (k >= N) break;
      vis[k] = true;
      if (i % prime[j]) mui[k] = -mui[i], sigma[k] = sigma[i] << 1, d[k] = 1;
      else { mui[k] = 0, sigma[k] = sigma[i] / (d[i] + 1) * (d[i] + 2), d[k] = d[i] + 1; break; }
    }
  }
  for (int i = 1; i < N; ++ i) sum[i] = sum[i - 1] + mui[i] * mui[i], sigma[i] += sigma[i - 1];
}

inline long long get_mui(long long n) {
  if (n < N) return sum[n];
  int q = sqrt(n); long long ret = 0;
  for (register int i = 1; i <= q; ++ i) if (mui[i]) ret += mui[i] * (n / i / i);
  return ret;
}

inline long long get_sigma(long long n) {
  if (n < N) return sigma[n];
  long long ret = 0;
  for (register long long i = 1, j; i <= n; i = j + 1)
    j = n / (n / i), ret += (j - i + 1) * (n / i);
  return ret;
}

long long calc(long long n) {
  int q = sqrt(n); long long ret = 0;
  for (register int i = 1; i <= q; ++ i) if (mui[i]) ret += get_sigma(n / i);
  for (register long long i = q + 1, j, cur, last = sum[q]; i <= n; i = j + 1, last = cur)
    j = n / (n / i), cur = get_mui(j), ret += (cur - last) * get_sigma(n / i);
  return ret;
}

signed main() {
  long long T; read(T);
  if (T > 800) N = 20000; else N = 100000000;
  init();
  while (T --) read(n), printf("%lld\n", calc(n));
  return 0;
}