Official

F - kirinuki Editorial by physics0523


数える部分長方形のうち、下の辺を固定することを考えます。
この時、各列について (そのマスから始めて) 上に何個 . が連続しているかを数えておくことができます。

例えば以下のグリッドで下の辺を最下部に固定した時、各列について上に .\(h=(0,2,1,4,5)\) 個続いています。

...#.
.#...
###..
..#..
#....

この列 \(h\) に対して \(O(M)\) で長方形の数を求められれば、全体で \(O(NM)\) でこの問題を解くことができます。

先程の列に対し、 Cartesian Tree を構築します。 参考1 参考2 Wikipedia(英語版)
但し、何らかの方法で \(h_i\) が等しい要素同士もタイブレークし、その中で順序が付いているものとします。
この Cartesian Tree から、各要素 \(i\) について以下の条件を満たす \(l,r\) を求めることができます。

  • \(j=l-1\) について \(h_j < h_i\)
  • \(j=l,l+1,\dots,i-1,i+1\dots,r\) について \(h_j > h_i\)
  • \(j=r+1\) について \(h_j < h_i\)

この \(l,r\) を使えば、 \(i\) 列目を含むかつ \(h_i\) が最小の序列を持つような長方領域を数え上げることができます。このとき、長方領域の高さは \(h_i\) 以下に制限されます。
横幅が小さい領域では \(h_i\) に高さが制約され、横幅が大きい領域では \(K\) に高さが制約されます。
横幅ごとに、その横幅を満たすように選択する方法は \((1,2,3,\dots,x-1,x,\dots,x,x-1,\dots,2,1)\) 通りという形になります。
以降の詳細は省略しますが、この係数列を満たすように累積和をうまく準備しておくことで、各 \(i\) ごとに \(O(1)\) で数えるべき長方領域の個数を求めることができます。

以上より、この問題を時間計算量 \(O(NM)\) で解くことが出来ました。

実装例 (C++):

#include<bits/stdc++.h>

using namespace std;

// https://nyaannyaan.github.io/library/tree/cartesian-tree.hpp.html
template <typename T>
pair<vector<vector<int>>, int> CartesianTree(vector<T> &a) {
  int N = (int)a.size();
  vector<vector<int>> g(N);
  vector<int> p(N, -1), st;
  st.reserve(N);
  for (int i = 0; i < N; i++) {
    int prv = -1;
    while (!st.empty() && a[i] < a[st.back()]) {
      prv = st.back();
      st.pop_back();
    }
    if (prv != -1) p[prv] = i;
    if (!st.empty()) p[i] = st.back();
    st.push_back(i);
  }
  int root = -1;
  for (int i = 0; i < N; i++) {
    if (p[i] != -1)
      g[p[i]].push_back(i);
    else
      root = i;
  }
  return make_pair(g, root);
}

int main(){
  long long n,m,k;
  cin >> n >> m >> k;
  vector<string> s(n);
  for(auto &nx : s){
    cin >> nx;
  }

  vector<vector<long long>> rw(6,vector<long long>(m+1,0));
  for(long long i=1;i<=m;i++){
    rw[0][i]=rw[0][i-1]+(k/i)*i;
    rw[1][i]=rw[1][i-1]+(k/i);
    rw[2][i]=rw[2][i-1]+(k/i)*(m+1-i);
    rw[3][i]=rw[3][i-1]+1*i;
    rw[4][i]=rw[4][i-1]+1;
    rw[5][i]=rw[5][i-1]+1*(m+1-i);
  }

  long long res=0;
  vector<long long> h(m,0);
  for(long long d=0;d<n;d++){
    for(long long i=0;i<m;i++){
      if(s[d][i]=='#'){h[i]=0;}
      else{h[i]++;}
    }
    auto [g,r]=CartesianTree(h);
    vector<long long> sq;
    queue<long long> q;
    q.push(r);
    while(!q.empty()){
      auto od=q.front(); q.pop();
      sq.push_back(od);
      for(auto &nx : g[od]){
        q.push(nx);
      }
    }
    reverse(sq.begin(),sq.end());
    vector<pair<long long,long long>> vp;
    for(long long i=0;i<m;i++){vp.push_back({i,i});}
    for(auto &nx : sq){
      for(auto &ny : g[nx]){
        vp[nx].first=min(vp[nx].first,vp[ny].first);
        vp[nx].second=max(vp[nx].second,vp[ny].second);
      }
      long long ch=h[nx];
      long long m0=1;
      long long m1=min(nx-vp[nx].first,vp[nx].second-nx)+1;
      long long m2=max(nx-vp[nx].first,vp[nx].second-nx)+1;
      long long m3=(vp[nx].second-vp[nx].first)+1;
      if(ch==0){continue;}
      {
        // unlimited
        long long l=(k/(ch+1))+1,r=1e9;
        {
          long long cl=max(l,m0),cr=min(r,m1);
          l=max(l,cr+1);
          if(cl<=cr){
            res+=rw[0][cr];
            res-=rw[0][cl-1];
          }
        }
        {
          long long cl=max(l,m1+1),cr=min(r,m2);
          l=max(l,cr+1);
          if(cl<=cr){
            res+=rw[1][cr]*m1;
            res-=rw[1][cl-1]*m1;
          }
        }
        {
          long long cl=max(l,m2+1),cr=min(r,m3);
          l=max(l,cr+1);
          if(cl<=cr){
            res+=rw[2][cr];
            res-=rw[2][cl-1];
            res-=rw[1][cr]*(m-m3);
            res+=rw[1][cl-1]*(m-m3);
          }
        }
      }
      {
        // limited
        long long l=1,r=(k/(ch+1));
        {
          long long cl=max(l,m0),cr=min(r,m1);
          l=max(l,cr+1);
          if(cl<=cr){
            res+=rw[3][cr]*ch;
            res-=rw[3][cl-1]*ch;
          }
        }
        {
          long long cl=max(l,m1+1),cr=min(r,m2);
          l=max(l,cr+1);
          if(cl<=cr){
            res+=rw[4][cr]*m1*ch;
            res-=rw[4][cl-1]*m1*ch;
          }
        }
        {
          long long cl=max(l,m2+1),cr=min(r,m3);
          l=max(l,cr+1);
          if(cl<=cr){
            res+=rw[5][cr]*ch;
            res-=rw[5][cl-1]*ch;
            res-=rw[4][cr]*(m-m3)*ch;
            res+=rw[4][cl-1]*(m-m3)*ch;
          }
        }
      }
    }
  }
  cout << res << "\n";
  return 0;
}

posted:
last update: