LibreOJ 509 动态几何问题

发表于 | 分类于 |

\[ \begin{aligned} &\quad\; \sum_{i=1}^n\sum_{j=1}^m [ij\text{ is perfect square}] \\ &= \sum_{k=1}^{\min(n,m)} \mu^2(k) \left\lfloor\sqrt{\frac nk}\right\rfloor \left\lfloor\sqrt{\frac mk}\right\rfloor \end{aligned} \]

考虑 \(\left\lfloor\sqrt{\frac nk}\right\rfloor\) 有多少种不同的取值。
容易发现以 \(n^{1/3}\) 为分界线可知有 \(O(n^{1/3})\) 种。

接下来处理计算 \(\mu^2\) 的前缀和。易知我们需要所有 \(\left\lfloor\frac n{x^2}\right\rfloor\) 处的 \(\mu^2\) 的前缀和。

注意到 \[ \sum_{i=1}^n \mu^2(i) = \sum_{i=1}^{\lfloor\sqrt n\rfloor} \mu(i) \left\lfloor\frac n{i^2}\right\rfloor \]

先不管怎么求 \(\mu\)。考虑设阈值 \(B\)：对 \(\left\lfloor\frac n{x^2}\right\rfloor \le B\)，线性筛；对 \(\left\lfloor\frac n{x^2}\right\rfloor > B\)，使用 \(O\left(\left\lfloor\frac n{x^2}\right\rfloor^{1/3}\right)\) 的整除分块计算。
后面的复杂度大约是 \[ \begin{aligned} &\quad\; \int_0^{ \sqrt{\frac nB} } \left(\frac n{x^2}\right)^{1/3} \,\mathrm d x \\ &= \left.3 n^{1/3} x^{1/3}\right|_{x=0}^{ \sqrt{\frac nB} } \\ &= O(n^{1/2} B^{-1/6}) \end{aligned} \]

平衡一下，取 \(B = n^{3/7}\) 可得复杂度为 \(O(n^{3/7})\)。

然后考虑求 \(\mu\)。我们需要的位置是 \(\left\lfloor\sqrt{ \frac n{x^2y} }\right\rfloor = \left\lfloor\frac{\left\lfloor\sqrt{\frac ny}\right\rfloor}x\right\rfloor\)。
设阈值 \(T\)，\(\le T\) 的部分用线性筛，\(> T\) 的部分用杜教筛，那么后面部分的复杂度大约是 \[ \begin{aligned} &\quad\; \int_0^{ \frac n{T^2} } \left(\frac ny\right)^{1/3} \mathrm dy \\ &= \left.\frac 32 n^{1/3} y^{2/3} \right|_{y=0}^{ \frac n{T^2} } \\ &= O(n T^{-4/3}) \end{aligned} \]

平衡一下，取 \(T = n^{3/7}\) 可得复杂度为 \(O(n^{3/7})\)。

然后乱写一下就过了，虽然常数贼大。

代码： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include <cmath>
#include <cstdio>
#include <algorithm>
#include <unordered_map>

using ll = long long;

using namespace std;

const int LIM = 9e6;
ll n, m;
int lim;
int vis[LIM + 5], cnt, prime[LIM + 5];
int mu[LIM + 5], mu2[LIM + 5];
unordered_map<ll, ll> memMu, memMu2;
ll sumMu(ll n) {
  if (n <= lim)
    return mu[n];
  if (memMu.count(n))
    return memMu[n];
  ll ret = 1;
  for (ll l = 2, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ret -= (r - l + 1) * sumMu(n / l);
  }
  return memMu[n] = ret;
}
ll sumMu2(ll n) {
  if (n <= lim)
    return mu2[n];
  if (memMu2.count(n))
    return memMu2[n];
  ll ret = 0;
  int nSqrt = sqrt(n);
  for (int l = 1, r; l <= nSqrt; l = r + 1) {
    ll v = n / ((ll)l * l);
    r = sqrt(n / v);
    ret += (sumMu(r) - sumMu(l - 1)) * v;
  }
  return memMu2[n] = ret;
}

ll ans;

int main() {
  scanf("%lld%lld", &n, &m);
  if (n > m)
    swap(n, m);
  lim = pow(m, 3.0 / 7);
  mu[1] = mu2[1] = 1;
  for (int i = 2; i <= lim; ++i) {
    if (!vis[i])
      prime[++cnt] = i, mu[i] = -1;
    for (int j = 1; j <= cnt && i * prime[j] <= lim; ++j) {
      vis[i * prime[j]] = 1;
      if (!(i % prime[j]))
        break;
      mu[i * prime[j]] = -mu[i];
    }
    mu2[i] = mu2[i - 1] + (bool)mu[i],
    mu[i] += mu[i - 1];
  }
  for (ll l = 1, r; l <= n; l = r + 1) {
    int v0 = sqrt(n / l), v1 = sqrt(m / l);
    r = min(n / ((ll)v0 * v0), m / ((ll)v1 * v1));
    ans += (sumMu2(r) - sumMu2(l - 1)) * v0 * v1;
  }
  printf("%lld\n", ans);
}