题解 Matrix and GCD

本题解仅提供F题的另一种求解全一子矩阵个数的方法

众所周知， 求解全一子矩阵的问题可以用 单调栈 解决， 例如最大子矩阵， 子矩阵个数

而有另一种工具也可以有效而优雅地解决该类问题， 那便是 笛卡尔树

(上面的超链接里有介绍笛卡尔树，并且提供了求解 最大子矩阵 的做法和例题)

如果你熟悉了 笛卡尔树求解最大子矩阵 的做法，相信 求解子矩阵的数量 也很好理解。

假设以当前行为底列连续全1列的高分别为 {}

那么对这k个数字建立小根堆笛卡尔树之后，对于当前分治中心u，区间假设为

计算经过u这一点的区间，满足

那么区间数量即为

因为是小根堆，所以目前任何区间 ，都可以向上延展最多至

从而可以得出 经过当前分治中心u的矩阵数量为 =

// 分治代码
long long sum;
int dfs(int u) {
  if (!u) return 0;
  int L = dfs(t[u].ls) + 1;
  int R = dfs(t[u].rs) + 1;
  sum += 1LL * L * R * t[u].w;
  return L + R - 1;
}

// F题完整代码
#include <bits/stdc++.h>
using namespace std;
const int N = 1005;
const int M = 1000005;
#define For(i,j,k) for(int i=(int)(j);i<=(int)(k);i++)
namespace Sieve {
  int cnt[M];
  void sieve(int n) {
    iota(cnt, cnt + n + 1, 0);
    for (int i = 1; i <= n; i++)
      for (int j = i << 1; j <= n; j += i)
        cnt[j] -= cnt[i];
  }
}
struct Node { int w, f, ls, rs; } t[N];
int build(int n) { // 建树
  for (int i = 1; i <= n; i++) {
    int k = i - 1;
    while (t[k].w > t[i].w) k = t[k].f;
    t[i].ls = t[k].rs, t[t[i].ls].f = i;
    t[i].f = k, t[k].rs = i;
  }
  return t[0].rs;
}
long long sum;
int dfs(int u) { // 分治
  if (!u) return 0;
  int L = dfs(t[u].ls) + 1;
  int R = dfs(t[u].rs) + 1;
  sum += 1LL * L * R * t[u].w;
  return L + R - 1;
}
int n, m, k, a[N][N];
short X[M], Y[M];
short f[N][N];
long long solve(int s) {
  sum = 0;
  for (int i = s; i <= n * m; i += s) {
    f[X[i]][Y[i]] = -1;
  }
  short x, y, cur;
  for (int i = s; i <= n * m; i += s) {
    x = X[i], y = Y[i], cur = 0;
    if (f[x][y] == -1 && f[x-1][y] == 0)
      while (f[x][y] == -1) f[x][y] = ++cur, ++x;
  }
  for (int i = s; i <= n * m; i += s) {
    x = X[i], y = Y[i];
    if (f[x][y] > 0 && f[x][y-1] == 0) {
      t[k = 0] = { 0, 0, 0, 0 };
      for (int j = y; f[x][j] > 0; j++) {
        t[++k] = { f[x][j], 0, 0, 0 };
        f[x][j] = 0;
      }
      dfs(build(k));
    }
  }
  return sum;
}
int main() {
  ios_base::sync_with_stdio(false);
  cin >> n >> m;
  Sieve::sieve(n * m);
  For(i, 1, n) For(j, 1, m)
    cin >> k, X[k] = i, Y[k] = j;
  long long res = 0;
  For(i, 1, n * m) res += solve(i) * Sieve::cnt[i];
  cout << res;
  return 0;
}