LibreOJ 509 动态几何问题

n∑i=1m∑j=1[ij is perfect square]=min(n,m)∑k=1μ2(k)⌊√nk⌋⌊√mk⌋

考虑 ⌊√nk⌋ 有多少种不同的取值。
容易发现以 n1/3 为分界线可知有 O(n1/3) 种。

接下来处理计算 μ2 的前缀和。易知我们需要所有 ⌊nx2⌋ 处的 μ2 的前缀和。

注意到 n∑i=1μ2(i)=⌊√n⌋∑i=1μ(i)⌊ni2⌋

先不管怎么求 μ。考虑设阈值 B：对 ⌊nx2⌋≤B，线性筛；对 ⌊nx2⌋>B，使用 O(⌊nx2⌋1/3) 的整除分块计算。
后面的复杂度大约是 ∫√nB0(nx2)1/3dx=3n1/3x1/3|√nBx=0=O(n1/2B−1/6)

平衡一下，取 B=n3/7 可得复杂度为 O(n3/7)。

然后考虑求 μ。我们需要的位置是 ⌊√nx2y⌋=⌊⌊√ny⌋x⌋。
设阈值 T，≤T 的部分用线性筛，>T 的部分用杜教筛，那么后面部分的复杂度大约是 ∫nT20(ny)1/3dy=32n1/3y2/3|nT2y=0=O(nT−4/3)

平衡一下，取 T=n3/7 可得复杂度为 O(n3/7)。

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

代码：

#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);
}