[HAOI2011]Problem b

Description

对于给出的 (n) 个询问，每次求有多少个数对 ((x,y))，满足 (a le x le b)，(c le y le d)，且 (gcd(x,y) = k)，(gcd(x,y)) 函数为 (x) 和 (y) 的最大公约数。

其中 (1 le n,k le 5 imes 10^4)，(1 le a le b le 5 imes 10^4)，(1 le c le d le 5 imes 10^4)。

Solution

首先肯定是自然而然地想到分成四块容斥，即

[egin{split} sum limits_{i=a}^{b} sum limits_{j=c}^{d} left[ gcd(i,j) = k ight] &= sum limits_{i=1}^{b} sum limits_{j=1}^{d} left[ gcd(i,j) =k ight] \ &- sum limits_{i=1}^{a} sum limits_{j=1}^{d} left[ gcd(i,j) =k ight]\ &- sum limits_{i=1}^{b} sum limits_{j=1}^{c} left[ gcd(i,j) =k ight]\ &+ sum limits_{i=1}^{a} sum limits_{j=1}^{c} left[ gcd(i,j) =k ight] end{split} ]

这样，每一部分都是 (sum limits_{i=1}^{m} sum limits_{j=1}^{m} left[ gcd(i,j) = k ight]) 的形式。

接下来考虑如何将式子进行变换：

[sum limits_{i=1}^{m} sum limits_{j=1}^{m} left[ gcd(i,j) = k ight] = sum limits_{i=1}^{lfloor frac{n}{k} floor} sum limits_{j=1}^{lfloor frac{m}{k} floor} left[ gcd(i,j) = 1 ight] ]

看到右边的 (left[ gcd(i,j) = 1 ight]) 是不是感觉很熟悉？根据在莫比乌斯反演学习笔记中提到的反演结论：

[left[ gcd(i,j) = 1 ight] iff sum limits_{d mid gcd(i,j)} mu(d) ]

可以继续化简：

[sum limits_{i=1}^{lfloor frac{n}{k} floor} sum limits_{j=1}^{lfloor frac{m}{k} floor} sum limits_{d mid gcd(i,j)} mu(d) ]

考虑变换求和顺序：

[sum limits_{d=1}^{min(n,m)} mu(d) sum limits_{i=1}^{lfloor frac{n}{k} floor} left[ d mid i ight] sum limits_{j=1}^{lfloor frac{m}{k} floor} left[ d mid j ight] ]

易知在 (left[ 1, n ight]) 中约数包含 (d) 的数有 (lfloor frac{n}{d} floor) 个，于是式子可以化简为：

[sum limits_{d=1}^{min(n,m)} mu(d) lfloor frac{lfloor frac{n}{k} floor}{d} floor lfloor frac{lfloor frac{m}{k} floor}{d} floor ]

再由莫比乌斯反演学习笔记中提到的引理1，可得：

[sum limits_{d=1}^{min(n,m)} mu(d) lfloor frac{n}{kd} floor lfloor frac{m}{kd} floor ]

所以最终，就得到了

[sum limits_{i=a}^{b} sum limits_{j=c}^{d} left[ gcd(i,j) = k ight] = sum limits_{d=1}^{min(n,m)} mu(d) lfloor frac{n}{kd} floor lfloor frac{m}{kd} floor ]

这个式子就可以利用数论分块进行求解了。

Code

Talk is cheap. Show me the code.

#include <bits/stdc++.h>
using namespace std;

int ty() {
  int x = 0, f = 1; char ch = getchar();
  while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
  while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
  return x * f;
}

const int _ = 5e4 + 10;
int tot, vis[_], p[_], mu[_], sum[_];

void getMu(int lim = 5e4) {
  mu[1] = 1;
  for (int i = 2; i <= lim; ++i) {
    if (!vis[i]) p[++tot] = i, mu[i] = -1;
    for (int j = 1; j <= tot && p[j] * i <= lim; ++j) {
      vis[i * p[j]] = 1;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
  for (int i = 1; i <= lim; ++i) sum[i] = sum[i - 1] + mu[i];
}

int calc(int n, int m, int k) {
  n /= k, m /= k;
  int ret = 0;
  for (int l = 1, r; l <= min(n, m); l = r + 1) {
    r = min(n / (n / l), m / (m / l));
    ret += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
  }
  return ret;
}

int main() {
  getMu();
  int T = ty(), a, b, c, d, k;
  while (T--) {
    a = ty(), b = ty(), c = ty(), d = ty(), k = ty();
    printf("%d
", calc(b, d, k) - calc(a - 1, d, k) - calc(b, c - 1, k) +
                       calc(a - 1, c - 1, k));
  }
  return 0;
}