BZOJ3994 [SDOI2015]约数个数和

Description

 设d(x)为x的约数个数，给定N、M，求$\sum_{i=1}^n\sum_{j=1}^md(ij)$

Input

输入文件包含多组测试数据。

第一行，一个整数T，表示测试数据的组数。

接下来的T行，每行两个整数N、M。

Output

 T行，每行一个整数，表示你所求的答案。

Sample Input

2
7 4
5 6

Sample Output

110
121

HINT

 1<=N, M<=50000

1<=T<=50000

Solution

计算$d(ij)$时，我们把$ij$的每个约数$d$映射到$(gcd(d, i), \frac{d}{gcd(d, i)})$，那么这两个数一定分别是$i, j$的因数，且$(a, b)$对应一个因数当且仅当$gcd(\frac ia, b) = 1$，所以

$$d(ij) = \sum_{x|i}\sum_{y|j} [gcd(\frac ix, y) = 1] = \sum_{x'|i}\sum_{y|j} [gcd(x', y) = 1]$$

于是

$$\begin{aligned}
 \sum_{i=1}^N\sum_{j=1}^Md(ij)
&=\sum_{i=1}^N\sum_{j=1}^M\sum_{x|i}\sum_{y|j} [gcd(x, y) = 1]\\
&=\sum_{x=1}^N\sum_{y=1}^M [gcd(x, y) = 1]\left\lfloor\frac Nx\right\rfloor \left\lfloor\frac My\right\rfloor
\end{aligned}$$

我们令$f(d) = \sum_{x=1}^N\sum_{y=1}^M [gcd(x, y) = d]\left\lfloor\frac Nx\right\rfloor \left\lfloor\frac My\right\rfloor$，有

$$\begin{aligned}
\sum_{n|d}f(d)
&= \sum_{x=1}^N\sum_{y=1}^M [n|gcd(x, y)]\left\lfloor\frac Nx\right\rfloor \left\lfloor\frac My\right\rfloor\\
&= \sum_{i=1}^{\left\lfloor\frac Nn\right\rfloor}\sum_{j=1}^{\left\lfloor\frac Mn\right\rfloor} \left\lfloor\frac{\left\lfloor\frac Nn\right\rfloor}i \right\rfloor\left\lfloor\frac{\left\lfloor\frac Mn\right\rfloor}j\right\rfloor
\end{aligned}$$

如果我们令$t(n) = \sum_{i=1}^n \left\lfloor\frac ni\right\rfloor$，那上式就等于$t(\left\lfloor\frac Nn\right\rfloor)t(\left\lfloor\frac Mn\right\rfloor)$

于是$f(n) = \sum_{n|d} \mu(\frac dn)t(\left\lfloor\frac Nd\right\rfloor)t(\left\lfloor\frac Md\right\rfloor)$

$ans = f(1) = \sum_{n=1}^{min(N, M)}\mu(n)t(\left\lfloor\frac Nn\right\rfloor)t(\left\lfloor\frac Mn\right\rfloor)$

预处理$\mu(n)$的前缀和、$O(n\sqrt n)$预处理所有$t(n)$，查询时$O(\sqrt n)$即可。

代码：

#include <algorithm>
#include <cstdio>
typedef long long LL;
const int N = 50050;
LL t[N];
LL calcT(int n) {
  LL ans = 0;
  for (int i = 1, last; i <= n; i = last + 1) {
    last = n / (n / i);
    ans += n / i * (last - i + 1);
  }
  return ans;
}
bool mark[N];
int pr[N], pcnt = 0, mu[N], S[N];
void getMu() {
  mu[1] = 1;
  for (int i = 2; i < N; ++i) {
    if (!mark[i]) mu[pr[pcnt++] = i] = -1;
    for (int j = 0; j < pcnt && (LL)i * pr[j] < N; ++j) {
      mark[i * pr[j]] = 1;
      if (!(i % pr[j])) {
        mu[i * pr[j]] = 0;
        break;
      }
      mu[i * pr[j]] = -mu[i];
    }
  }
  for (int i = 1; i < N; ++i) S[i] = S[i - 1] + mu[i];
}
LL solve(int n, int m) {
  LL ans = 0;
  for (int i = 1, last; i <= n && i <= m; i = last + 1) {
    last = std::min(n / (n / i), m / (m / i));
    ans += t[n / i] * t[m / i] * (S[last] - S[i - 1]);
  }
  return ans;
}
int main() {
  for (int i = 1; i < N; ++i) t[i] = calcT(i);
  getMu();
  int T, n, m;
  scanf("%d", &T);
  while (T--) {
    scanf("%d%d", &n, &m);
    printf("%lld\n", solve(n, m));
  }
  return 0;
}