BZOJ3994 [SDOI2015]约数个数和

Description

 设d(x)为x的约数个数，给定N、M，求$sum_{i=1}^nsum_{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]$$

于是

$$egin{aligned}
 sum_{i=1}^Nsum_{j=1}^Md(ij)
&=sum_{i=1}^Nsum_{j=1}^Msum_{x|i}sum_{y|j} [gcd(x, y) = 1]\
&=sum_{x=1}^Nsum_{y=1}^M [gcd(x, y) = 1]leftlfloorfrac Nx ight floor leftlfloorfrac My ight floor
end{aligned}$$

我们令$f(d) = sum_{x=1}^Nsum_{y=1}^M [gcd(x, y) = d]leftlfloorfrac Nx ight floor leftlfloorfrac My ight floor$，有

$$egin{aligned}
sum_{n|d}f(d)
&= sum_{x=1}^Nsum_{y=1}^M [n|gcd(x, y)]leftlfloorfrac Nx ight floor leftlfloorfrac My ight floor\
&= sum_{i=1}^{leftlfloorfrac Nn ight floor}sum_{j=1}^{leftlfloorfrac Mn ight floor} leftlfloorfrac{leftlfloorfrac Nn ight floor}i ight floorleftlfloorfrac{leftlfloorfrac Mn ight floor}j ight floor
end{aligned}$$

如果我们令$t(n) = sum_{i=1}^n leftlfloorfrac ni ight floor$，那上式就等于$t(leftlfloorfrac Nn ight floor)t(leftlfloorfrac Mn ight floor)$

于是$f(n) = sum_{n|d} mu(frac dn)t(leftlfloorfrac Nd ight floor)t(leftlfloorfrac Md ight floor)$

$ans = f(1) = sum_{n=1}^{min(N, M)}mu(n)t(leftlfloorfrac Nn ight floor)t(leftlfloorfrac Mn ight floor)$

预处理$mu(n)$的前缀和、$O(nsqrt 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
", solve(n, m));
  }
  return 0;
}