Baby Stepsなブログ

競プロとか。間違ったこと書いてあったら@pi0nep1oneにご連絡ください。

行列累乗についてのメモと出題例

最近DPの漸化式を行列式に変換するセンスを磨くために関連問題を解いており、問題や行列累乗全般について雑多なメモを残しておこうと思う。

フィボナッチ数列のN項目を高速に求める

行列累乗の有名な活用例としてよくフィボナッチ数列が挙げられる。

オンラインジャッジで全く同じ問題は見つけられなかったが、蟻本にその記載がある。

フィボナッチ数列の漸化式は以下の様になる。

  • dp[0]=dp[1]=1
  • dp[i+2]=dp[i+1]+dp[i]

この遷移をナイーブに実装すると第N項を求めるためにO(N)の時間が掛かる。

これを高速に求めるため、上記の漸化式を行列積の形に書き換える。(f_iをフィボナッチ数列の第i項目とすると)

{f_1, f_0} = {{1,1}, {1,0}} * {1,0}

更に右辺に {{1,1}, {1,0}} を掛けると {f_2, f_1}、更に掛けると {f_3, f_2} を求めることができる。

つまり第N項目を求めるためには {{1,1}, {1,0}} のN乗を求めればよい。累乗を高速に求める手法に繰り返し二乗法というものがあるが、これが行列にも応用することができる。この高速に行列の累乗を求める手法を競プロ界隈では「行列累乗」と呼称しているらしい。

以下は、行列累乗によりフィボナッチ数列を求めるための実装

フィボナッチ数列

#define mat vector<vector<ll>>

/// 行列積
mat mat_mul(mat &a, mat &b) {
  mat res(a.size(), vector<ll>(b[0].size()));
  for (int i = 0; i < a.size(); i++) {
    for (int j = 0; j < b[0].size(); j++) {
      for (int k = 0; k < b.size(); k++) {
        (res[i][j] += a[i][k] * b[k][j]) %= MOD;
      }
    }
  }
  return res;
}

/// 行列累乗
mat mat_pow(mat a, long long n) {
  mat res(a.size(), vector<ll>(a.size()));
  // 単位行列で初期化
  for (int i = 0; i < a.size(); i++)
    res[i][i] = 1;
  // 繰り返し二乗法
  while (n > 0) {
    if (n & 1) res = mat_mul(a, res);
    a = mat_mul(a, a);
    n >>= 1;
  }
  return res;
}

signed main()
{
  ll n;
  cin >> n;
  
  // フィボナッチ数列の行列累乗
  mat m(2, vc<ll>(2));
  m[0][0] = m[0][1] = m[1][0] = 1;
  m[1][1] = 0;
  m = mat_pow(m, n);
  ll ans = m[1][0];
  cout << ans << endl;
}

例題をいくつか

以下に累積累乗によって解ける問題を載せておく。

下に行くほど難しくなっている。

ABC 009 D - 漸化式

atcoder.jp

この問題で漸化式を立てると以下の様になる。

A[i] = (C[0]&A[i-1]) ^ (C[1]&A[i-2]) ^ ...

この遷移を観察すると、自身より前の項によって自身の項が定まるためフィボナッチ数列DPに似ている。

また、求めたい第M項のMの値が大きい(最大109)かつ、自身を求めるために参照する自身より前の項の数Kが小さい(最大100)であるということから計算量から行列累乗が解法の候補に挙がる。(行列累乗で解けるなら、Kが行列の列幅行幅に対応するので、その計算量はO(K3logM)と十分に高速である。)

ここで、この問題では演算が加法、乗法ではなく、排他的論理和(XOR)、論理積(AND)を使うが半環を満たすため行列累乗が使える。

累乗する行列は以下の様になる。

※ 以下の画像はは解説スライドのP.102の引用です。

この問題では乗法が論理積に置き換わっているので、ここでの1とは、論理積を行った結果が変化しない値である。Aの要素は32ビット符号なし整数なので、32ビット1が立った値、(1<<32)-1 などを扱えばよい。

f:id:pione8888:20210404202641j:plain

解説スライド

www.slideshare.net

解答コード

atcoder.jp

ARC 050 C - LCM 111

atcoder.jp

X桁のレピュニット数をf(X)と表すことにする。

レピュニット数同士の剰余を観察してみると(A>Bとする)

f(A)%f(B) = f(A%B)

という関係があることがわかる。レピュニット数同士の剰余がその桁数の剰余に対応することから、レピュニット数同士のgcdはその桁数のgcdによって求めることができる。つまり、

gcd(f(A), f(B)) = f(gcd(A, B))

という関係が成り立つ。

gcd(A, B)をGとすると求める解は、

f(A)*f(B)/f(G)

で求めることができるので、後は行列累乗でf(A), f(B), f(G)を高速に求めれば良さそう...に思えるが、この問題で与えられる法Mは適当な正整数であり(f(G)と互いに素であるとは限らないため)f(G)の逆元を求められない可能性がある。

ひとまず f(A) を求めることを考える。

f(i+1) = f(i) * 10 + 1

なので、これを行列式で表すと以下の様になる。

{f(2), f(1)} = {{10,1},{0,1}}*{1,1}

よって、f(N)を求めるには、{{10,1},{0,1}}のN乗を行列累乗で求めればよいとわかる。

残る f(B)/f(G) について、これは逆元を使わずに求めることにする。f(B), f(G)はいずれもレピュニット数で f(B)>=f(G) なので、このような関係性の2数の除法について観察してみる。

f(6) / f(2) = 111111 / 11 = 10101

この 10101 という値はレピュニット数の遷移とほぼ同じで求めることができる。つまり、以下の様に表せる。

g(i+1) = g(i) * 10G + 1

これを行列式で表すと以下の様になる。

{g(2), g(1)} = {{10G,1},{0,1}}*{1,1}

B/G = g とすると、f(B)/f(G) は {{10G,1},{0,1}}*{1,1} をg乗することで求めることができる。

解答コード

atcoder.jp

EDPC R - Walk

atcoder.jp

この問題もKという値が大きい(最大1018)かつ、行列の列幅になるであろうNの値が小さい(最大50)ので、行列累乗で解けるならO(N3logK)時間になり間に合いそう。(さらにこの問題は数え上げである)

まず素直に漸化式を考えると、以下の様な遷移が考えられる。

dp[len][pos] := 経路長lenで終点がposであるものの通り数

  • dp[0][i] = 1 (初期化)
  • dp[len+1][to] += dp[len][from]

簡単のため N=2 という小さい場合について考える。有向辺は以下の様に与えられるとする。

0 0
1 0 // 頂点2->頂点1の有向辺

この行列を(i, j)=(頂点i <- 頂点j の有向辺が存在する)というように変換する。

0 1 // 頂点1<-頂点2の有向辺
0 0

この行列を使えば、以下の様に行列式で漸化式を表すことができる。

{dp[1][0], dp[1][1]} = {{0,1}, {0,0}} * {1,1}

経路長2の長さを求めるに右辺に {{0,1}, {0,0}} を掛ければよく、つまり経路長Kの長さを求めるにはこのK乗を求めればよい。

以上を有向グラフ全般に適用させるためには、与えられる行列Aから、(頂点i <- 頂点j の有向辺が存在する)行列を構築してそのK乗を求めればよい。

解答コード

atcoder.jp

yukicoder No.840 ほむほむほむら

yukicoder.me

h(ほ) -> m(む) -> r(ら) の並び順の部分文字列の個数がK倍個であるような文字列の個数を数え上げる問題。

この問題も、文字列の長さとして与えられるNが大きく(最大109)、Kが小さい(最大5)。

まず素直にDPの漸化式を考えてみる。

dp[n][i][j][k] := n文字で、hの個数がi個の倍数、hmの個数がj個の倍数、hmrの個数がkの倍数である場合の通り数(よって配列長はN*K3

  • dp[0][0][0][0] = 1 (初期化)
  • dp[n][(i+1)%K][j][k] += dp[n-1][i][j][k] // hを追加
  • dp[n][i][(j+i)%K][k] += dp[n-1][i][j][k] // mを追加
  • dp[n][i][j][(k+j)%K] += dp[n-1][i][j][k] // rを追加

これを行列式で合わらすために [i][j][k] の3次元配列をK進数によって1次元で表すことを考える。

つまり、hがi個、mがj個、rがk個を iK2+jK+k と表すことにする。この値は高々K3(<=53)と小さく、この行幅列幅の行列を累乗しても十分高速である。

あとは、行列mのm(to, from)を状態toへ状態fromからの遷移があるとして初期化し、このN乗を求めればN文字の各状態の通り数を計算できる。このうち、hmrの個数が0倍の個数の総和を求めればそれが解となる。

解答コード

yukicoder.me

yukicoder No.274 The Wall 解説

問題

yukicoder.me

2-SATについての前提

2-SAT問題の解き方

  1. 論理和論理包含に変換
  2. 1を元にimplication graphを作る
  3. 2のimplication graphを強連結成分分解する
  4. 3の結果、ある要素Xとその否定¬Xが強連結成分であるなら、その2-SATを満たす解は存在しない

式変換について

論理包含P→Qとは、Pが偽またはQが真であることを意味する。論理包含を利用することで論理和は、以下の様に書き換えることが可能。

P∨Q ⇔ ¬P→Q∧¬Q→P

implication graphについて

論理包含をの関係をグラフにおこしたものをimplication graphと言う。

implication graphにおいて強連結である頂点の真偽値は全て等しいという性質がある。(これは参考に載せた、2-SATを解いた時のメモに解説がある)

強連結である頂点は互いに行き来可能なので、(P, Q)が強連結であるなら、P→Q∧Q→P、つまりPとQは等しい。

この性質から、もしある変数Xとその否定¬Xが同じ強連結な集合に属するなら、そのimplication graphは矛盾することになり、その2-SATには解が無い事になる。

解説

1つの区間を2-SATにおける1つの変数と捉える。また、その変数の否定とは、区間を180度回転させたものとする。

すると、ある2つの区間が重なる可能性があるなら、その2つの区間を指す変数をA, Bとして、A∨Bという論理和の関係が成り立つ。(どちらか一方しか取れないため)

全ての2頂点間について、重なる可能性があるかを見ていく。(2頂点間の関係は180度回転させる場合を加味すると4パターンあることに注意する。)

上記の論理和の関係を、論理包含に変換しimplication graphを構築する。後は前提に記述した2-SATの解き方で解ける。

解答コード

struct StronglyConnectedComponents {
  int n;
  vector<vector<int>> G, rG;
  vector<int> order, component;
  vector<bool> used;
  void dfs(int v) {
    used[v] = 1;
    for (auto nv : G[v]) {
      if (!used[nv]) dfs(nv);
    }
    order.push_back(v);
  }
  void rdfs(int v, int k) {
    component[v] = k;
    for (auto nv : rG[v]) {
      if (component[nv] < 0) rdfs(nv, k);
    }
  }

  StronglyConnectedComponents(vector<vector<int>> &_G) {
    n = _G.size();
    G = _G;
    rG.resize(n);
    component.assign(n, -1);
    used.resize(n);
    for (int v = 0; v < n; v++) {
      for (auto nv : G[v]) rG[nv].push_back(v);
    }
    for (int v = 0; v < n; v++) if (!used[v]) dfs(v);
    int k = 0;
    reverse(order.begin(), order.end());
    for (auto v : order) if (component[v] == -1) rdfs(v, k), k++;
  }

  bool is_same(int u, int v) { return component[u] == component[v]; }

  vector<vector<int>> rebuild() {
    int N = *max_element(component.begin(), component.end()) + 1;
    vector<vector<int>> rebuildedG(N);
    set<pair<int, int>> connected;
    for (int v = 0; v < N; v++) {
      for (auto nv : G[v]) {
        if (component[v] != component[nv] and !connected.count(pair(v, nv))) {
          connected.insert(pair(v, nv));
          rebuildedG[component[v]].push_back(component[nv]);
        }
      }
    }
    return rebuildedG;
  }
};

#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i)
#define irep(i, m, n) for (long long i = (long long)(m); i < (long long)(n); ++i)
using ll = long long;
template<class t> using vc=vector<t>;

signed main()
{
  cin.tie( 0 ); ios::sync_with_stdio( false );
  
  int n,m; cin>>n>>m;
  vc<int> l(n), r(n);
  rep(i,n) cin>>l[i]>>r[i];
  
  vc<vc<int>> G(2*n);
  rep(i,n){
    irep(j,i+1,n){
      int l1=l[i], r1=r[i], l2=l[j], r2=r[j];
      int L1=m-r1-1, R1=m-l1-1, L2=m-r2-1, R2=m-l2-1;
      if(l1<=l2 and l2<=r1 or l1<=r2 and r2<=r1
        or l2<=l1 and l1<=r2){
        // A∨B <=> ¬A=>B∧¬B=>A
        G[i+n].push_back(j);
        G[j+n].push_back(i);
      }
      if(l1<=L2 and L2<=r1 or l1<=R2 and R2<=r1
        or L2<=l1 and l1<=R2){
        // A∨¬B <=> ¬A=>¬B∧B=>A
        G[i+n].push_back(j+n);
        G[j].push_back(i);
      }
      if(L1<=l2 and l2<=R1 or L1<=r2 and r2<=R1
        or l2<=L1 and L1<=r2){
        // ¬A∨B <=> A=>B∧¬B=>¬A
        G[i].push_back(j);
        G[j+n].push_back(i+n);
      }
      if(L1<=L2 and L2<=R1 or L1<=R2 and R2<=R1
        or L2<=L1 and L1<=R2){
        // ¬A∨¬B <=> A=>¬B∧B=>¬A
        G[i].push_back(j+n);
        G[j].push_back(i+n);
      }
    }
  }
  
  auto scc = StronglyConnectedComponents(G);
  rep(i,n){
    if(scc.is_same(i,i+n)){
      cout<<"NO"<<endl;
      return 0;
    }
  }
  cout<<"YES"<<endl;
}

参考

japlj.hatenadiary.org

論理包含に関する解説

applelife100.blogspot.com

強連結成分分解の実装と個人的なメモ

前提

有向グラフGの頂点の部分集合をSとする。

  • Sが強連結であるとは、Sの任意の2頂点間が行き来可能であること
  • Sが強連結成分であるとは、Sに他のどの頂点を追加してもそれ以上強連結にならないこと

強連結成分分解とは、強連結成分を1つの頂点にまとめること。こうしてまとめられた頂点でグラフを再構築すると、その有向グラフはDAGになる。

  • グラフをDAGに再構築し、DAGに対してDPするのは典型
  • 強連結成分分解を使って解くことができる有名問題として2-SATがある

強連結成分分解のアルゴリズムは単純で2回DFSを行うだけ。

1回目のDFSで各頂点に帰りがけの順を記録する。

2回目のDFSは元のグラフの辺を逆辺に置き換えたものに対して行う。1回目で採番した帰りがけの遅い順にDFSする。既に訪れた頂点を避けながらDFSすることで、その頂点から行き来できる頂点集合が1つの強連結成分となる。

2回のDFSで何をしているのか

これは参考に載せた、目指せグラフマスターのスライドP.11によくまとめられている。

  1. 強連結な集合は辺を逆辺に置き換えても強連結のまま
  2. DFSの帰りがけ順が遅い => 完成したDAGのトポロジカル順序において早い
  3. 逆辺のDFSなら今見ている頂点から(2のトポロジカル順序における)後の強連結成分を訪れない

3が重要で、このトポロジカル順序を無視して逆辺のグラフにDFSすると、ある連結成分から、自身よりトポロジカル順序が前の連結成分へ移動してしまう可能性があり、本来連結成分ではない集合を連結成分としてしまう可能性がある。

実装

struct StronglyConnectedComponents {
  int n;
  vector<vector<int>> G, rG;
  vector<int> order, component;
  vector<bool> used;
  void dfs(int v) {
    used[v] = 1;
    for (auto nv : G[v]) {
      if (!used[nv]) dfs(nv);
    }
    order.push_back(v);
  }
  void rdfs(int v, int k) {
    component[v] = k;
    for (auto nv : rG[v]) {
      if (component[nv] < 0) rdfs(nv, k);
    }
  }

  StronglyConnectedComponents(vector<vector<int>> &_G) {
    n = _G.size();
    G = _G;
    rG.resize(n);
    component.assign(n, -1);
    used.resize(n);
    for (int v = 0; v < n; v++) {
      for (auto nv : G[v]) rG[nv].push_back(v);
    }
    for (int v = 0; v < n; v++) if (!used[v]) dfs(v);
    int k = 0;
    reverse(order.begin(), order.end());
    for (auto v : order) if (component[v] == -1) rdfs(v, k), k++;
  }

  /// 頂点(u, v)が同じ強連結成分に含まれるか
  bool is_same(int u, int v) { return component[u] == component[v]; }

  /// 強連結成分を1つのノードに潰したグラフを再構築する
  vector<vector<int>> rebuild() {
    int N = *max_element(component.begin(), component.end()) + 1;
    vector<vector<int>> rebuildedG(N);
    set<pair<int, int>> connected;
    for (int v = 0; v < N; v++) {
      for (auto nv : G[v]) {
        if (component[v] != component[nv] and !connected.count(pair(v, nv))) {
          connected.insert(pair(v, nv));
          rebuildedG[component[v]].push_back(component[nv]);
        }
      }
    }
    return rebuildedG;
  }
};

検証用問題

onlinejudge.u-aizu.ac.jp

参考

p.82から強連結成分分解の解説

https://hcpc-hokudai.github.io/archive/graph_scc_001.pdf

www.slideshare.net

ダブリングによる接尾辞配列(Suffix-Array)の実装と個人的なメモ

AC Libraryに、より高速なライブラリが存在する。

以下に載せた実装は、自分の理解を深めるために実装したもの。(蟻本のものにアレンジを加えclass化した。計算量は変わらず、Manber & Myers の O(n log2 n) 時間のアルゴリズム

前提

接尾辞配列(Suffix-Array)とは、ある文字列の接尾辞(Suffix)を辞書順にソートしたもの。

ナイーブにSuffix-Arrayを構築しようとすると、文字列をsとして、空間計算量 |s|^2、時間計算 |s|^2 log |s| となる。後述するが、Suffix-Arrayを構築するのは文字列中に含まれる部分文字列の検索を高速に行うためであり、ナイーブなSuffix-Arrayの構築の時間計算量は、ナイーブな検索の時間計算量(|s|^2)を上回ってしまうので、このままでは構築する意味があまりない。また空間計算量も、扱う文字列の長さが105以上の問題ではMLEになってしまう。

そのため、Suffix-Array構築時の計算量を改善する必要がある。

Suffix-Arrayでもつ値は部分文字列ではなくSuffixの始点のindexにする。(これで空間計算量は|s|になる)

各indexからのSuffixの辞書順の順位が求められれば、その順位を使ってSuffix-Arrayのindexをソートすることができるので、rank配列を作り、これをダブリングによって段階的に更新していく。

以下は、rank配列を更新する実装部分についての簡単な説明。

rank配列の初期化

rankはそのindexの文字コードによって初期化している。rank[n](空文字)だけは、-1で初期化。

    for (int i = 0; i <= n; i++) {
      suffix_array[i] = i;
      rank[i] = i < n ? (int)s[i] : -1;
    }
確定した2*k文字まで長さのrankを元にSuffix-Arrayをソートする

k=1から始めて、各indexから2*k文字まで見た時の順位をrank配列に格納していく。

nを対象の文字列の長さとして、2*k>=nまで更新を行えば全てのSuffixのrankを求めたことになる。

ループ回数が log n 回、各シーケンス内のソートが n log n 時間なので、時間計算量は n log2 n

    vector<int> tmp(n + 1); // 一度tmp配列に
    for (k = 1; k <= n; k *= 2) {
      sort(suffix_array.begin(), suffix_array.end(), compare);
      tmp[suffix_array[0]] = 0;
      for (int i = 1; i <= n; i++) {
        tmp[suffix_array[i]] = tmp[suffix_array[i - 1]] + (compare(suffix_array[i - 1], suffix_array[i]) ? 1 : 0);
      }
      for (int i = 0; i <= n; i++) {
        rank[i] = tmp[i];
      }
    }
ソートに使う比較関数

k=1から始めて、2*kまでの長さの部分文字列の大小比較を行う。

    auto compare = [&](int a, int b) {
      if (rank[a] != rank[b]) return rank[a] < rank[b]; // 前半k文字のrankが異なるなら、後半k文字を比較するまでもなく大小が確定する
      else {
        // 前半k文字のrankが同等なら、後半のk文字のrankの大小で2*k文字の大小が決まる
        // +k文字目が、文字列をはみ出すなら空文字としてrank=-1とする
        auto rank_ak = a + k <= n ? rank[a + k] : -1;
        auto rank_bk = b + k <= n ? rank[b + k] : -1;
        return rank_ak < rank_bk;
      }
    };

実装

class SuffixArray {
private:
  string s;
  int n, k;
  vector<int> suffix_array;
  vector<int> rank;
  void create_suffix_array() {
    for (int i = 0; i <= n; i++) {
      suffix_array[i] = i;
      rank[i] = i < n ? (int)s[i] : -1;
    }
    auto compare = [&](int a, int b) {
      if (rank[a] != rank[b]) return rank[a] < rank[b];
      else {
        auto rank_ak = a + k <= n ? rank[a + k] : -1;
        auto rank_bk = b + k <= n ? rank[b + k] : -1;
        return rank_ak < rank_bk;
      }
    };
    vector<int> tmp(n + 1);
    for (k = 1; k <= n; k *= 2) {
      sort(suffix_array.begin(), suffix_array.end(), compare);
      tmp[suffix_array[0]] = 0;
      for (int i = 1; i <= n; i++) {
        tmp[suffix_array[i]] = tmp[suffix_array[i - 1]] + (compare(suffix_array[i - 1], suffix_array[i]) ? 1 : 0);
      }
      for (int i = 0; i <= n; i++) {
        rank[i] = tmp[i];
      }
    }
  }

public:
  SuffixArray(const string &_s) {
    s = _s;
    n = _s.size();
    suffix_array.resize(n + 1);
    rank.resize(n + 1);
    create_suffix_array();
  }

  // 文字列sの部分文字列として、パラメータで渡した文字列tが現れるかを返す
  bool contains(const string &t) {
    int l = 0, r = n;
    while (r - l > 1) {
      int mid = (l + r) / 2;
      int index = suffix_array[mid];
      if (s.compare(index, t.size(), t) < 0) l = mid;
      else r = mid;
    }
    return s.compare(suffix_array[r], t.size(), t) == 0;
  }
};
Suffix-Arrayを使った文字列検索

上記の実装のcontainsという関数は、文字列sの中にtが現れるかを判定している。

Suffix-Arrayがあることで二分探索によって部分文字列の検索を行えるので、O(|t| log |s|) 時間で判定を行える。

検証用問題

Suffix-Arrayを使った高速文字列検索によって解くことができる。(上記のcontainsにより検証可能)

onlinejudge.u-aizu.ac.jp

ABC168 F - . (Single Dot) 解説

制約から座標圧縮をすべきであることは想像しやすいが、この問題はかなり難しいと思った。

atcoder.jp

前提

  • 座標圧縮

pione.hatenablog.com

解説

座標圧縮を行う問題でよく見かけるパターンは以下の2つである。

  1. いくつかの長方形が与えられて、長方形が重なる領域の面積の総和はいくつか?
  2. いくつかの線分が与えられて、線分で区切られた領域の数はいくつか?

今回の問題は、線分が与えられてそれにより区切られた領域の面積を求めるというもので、上記の2点とは微妙に異なる。

1に類する問題は、座標圧縮後の状態をグリッドにマッピングすると、各マスと面積が対応するので対象のマスの面積をそれぞれ求めてその総和を求めればよい。

2に類する問題は、マスと線分を対応させてグリッドにマッピングし、BFSすることで求めることができる。(この時、グリッドの1マスは面積に対応しない)

今回の問題では、マスを線分に対応させつつ、移動できるマスからその面積を求める必要がある。

線分をグリッドにマッピングしたとき、線分間に空間があるかをわかるようにするため、与えられる座標を全て2倍にする。

実際にサンプル1で試すとわかるが、座標圧縮しなくても、与えられる線分の情報をそのままグリッドにマッピングすると、線分間に空間があるのかわからない。例えば、x軸方向に延びる線分で、x=0, x=1, x=2 の3つの線分をグリッドにマッピングすると、当然3行を連続して塗りつぶしてしまい、この行間に空間が存在するのかがわからない。

そのため、各座標を2倍することを考える。x=0, x=1, x=2 → x=0, x=2, x=4 このようにすると、各線分間には必ず1以上の差が生まれるため、適切に座標圧縮する(線分の座標の値+1を空白として登録する。前提参照)ことでグリッドに線分間の空間をマッピングすることができる。

グリッドのマスから面積を求める方法については結論から書くと、マス (x, y) が空白でかつ x, y が共に奇数であるマスだけを数えるようにすると良い。

あるマスが空白であるなら、そのマスの面積は (隣のY座標 - 今のY座標) * (隣のX座標 - 今のX座標) で求めることができる。

これは、前工程で線分の各座標値を2倍にし、その隣の座標を空白としたため、(x, y) のいずれかが偶数のマスは線分間の空白をマッピングするためのバッファであり、そのマスの面積は線分間の面積に対応しないため。

実装

#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i)
#define all(v) v.begin(), v.end()
using ll = long long;
template<class t> using vc=vector<t>;

template <typename T>
vector<T> compress(vector<T> &v1, vector<T> &v2, vector<T> &v3) {
  vector<T> vals;
  int n = v1.size(), m = v3.size();
  for (int i = 0; i < n; i++) {
    vals.push_back(v1[i]);
    vals.push_back(v2[i]);
    vals.push_back(v1[i] + 1);
    vals.push_back(v2[i] + 1);
  }
  for (int i = 0; i < m; i++) {
    vals.push_back(v3[i]);
    vals.push_back(v3[i] + 1);
  }
  vals.push_back(0);
  vals.push_back(1);
  sort(vals.begin(), vals.end());
  vals.erase(unique(vals.begin(), vals.end()), vals.end());
  for (int i = 0; i < n; i++) {
    v1[i] = lower_bound(vals.begin(), vals.end(), v1[i]) - vals.begin();
    v2[i] = lower_bound(vals.begin(), vals.end(), v2[i]) - vals.begin();
  }
  for (int i = 0; i < m; i++) {
    v3[i] = lower_bound(vals.begin(), vals.end(), v3[i]) - vals.begin();
  }
  return vals;
}

void inf() {
  cout << "INF" << endl;
  exit(0);
}

const ll dy[] = {0, 1, 0, -1, -1, 1, 1, -1};
const ll dx[] = {1, 0, -1, 0, 1, 1, -1, -1};

inline bool inside(ll y, ll x, ll H, ll W) {
    return (y >= 0 && x >= 0 && y < H && x < W);
}

signed main()
{
  cin.tie( 0 ); ios::sync_with_stdio( false );

  int n, m;
  cin >> n >> m;
  vc<ll> x1(n), x2(n), y1(n), x3(m), y2(m), y3(m);
  rep(i, n) cin >> x1[i] >> x2[i] >> y1[i];
  rep(i, m) cin >> x3[i] >> y2[i] >> y3[i];

  // 全ての座標の値を2倍
  rep(i, n) {
    x1[i] *= 2;
    x2[i] *= 2;
    y1[i] *= 2;
  }
  rep(i, m) {
    x3[i] *= 2;
    y2[i] *= 2;
    y3[i] *= 2;
  }

  // 座標圧縮
  auto X = compress(x1, x2, x3);
  auto Y = compress(y2, y3, y1);

  int h = X.size(), w = Y.size();

  // X軸方向、Y軸方向、別々にimos法を行う
  vc<vc<int>> g1(h, vc<int>(w)); // X軸方向
  vc<vc<int>> g2(h, vc<int>(w)); // Y軸方向

  rep(i, n) {
    g1[x1[i]][y1[i]]++;
    g1[x2[i] + 1][y1[i]]--;
  }
  rep(j, w) {
    rep(i, h - 1) { g1[i + 1][j] += g1[i][j]; }
  }
  rep(i, m) {
    g2[x3[i]][y2[i]]++;
    g2[x3[i]][y3[i] + 1]--;
  }
  rep(i, h) {
    rep(j, w - 1) { g2[i][j + 1] += g2[i][j]; }
  }

  // BFS で移動可能な領域の面積の総和を求める。h * w の範囲外に出られるなら INF
  vc<vc<bool>> visited(h, vc<bool>(w));
  queue<pair<int, int>> q;
  int sx = lower_bound(X.begin(), X.end(), 0) - X.begin();
  int sy = lower_bound(Y.begin(), Y.end(), 0) - Y.begin();
  q.emplace(sx, sy);
  visited[sx][sy] = 1;
  ll ans = 0;
  while (!q.empty()) {
    int x, y;
    tie(x, y) = q.front();
    q.pop();
    if (x % 2 == 1 and y % 2 == 1) ans += (X[x + 1] / 2 - (X[x] - 1) / 2) * (Y[y + 1] / 2 - (Y[y] - 1) / 2);
    rep(i, 4) {
      int nx = x + dx[i], ny = y + dy[i];
      if (inside(nx, ny, h, w) and g1[nx][ny] == 0 and g2[nx][ny] == 0 and visited[nx][ny] == 0) {
        visited[nx][ny] = 1;
        q.emplace(nx, ny);
      } else if (!inside(nx, ny, h, w)) {
        inf();
      }
    }
  }
  cout << ans << endl;
}

座標圧縮に関する解説と練習問題

前提

一次元座標圧縮はBITによる転倒数の数え上げでお馴染み。

二次元の場合は、x軸、y軸を独立に捉えることでほぼ同様に求めることができる。

二次元座標圧縮における注意点は、登場する座標の1つ隣の座標も登録する点にある。これは登場する座標だけで圧縮を行うと、元のグリッドの線分(や長方形)間に存在した空白を潰してしまう可能性があるため。(この図解が参考1のスライドにある)

座標圧縮の基本的な実装

  • パラメータに渡した座標をmutableに扱い、圧縮後の座標に置き換える
    • 座標圧縮が必要になる問題は、線分の始点と終点や長方形(3Dなら直方体)の対角など、同数の2つの座標が与えられる形式が多く、そのため引数に2つのリストを取るようにすると便利なことが多い(練習問題での compress の呼び出し元を参照)
  • 戻り値として、引数に渡した座標(とその1つ隣の座標)の圧縮前の座標をソートし重複を排除したリストvalsを返す

上記の様にすることで、圧縮後の値xから圧縮前の値を復元したい場合、vals[x] と記述するだけで簡単に行うことができる。

実際の使い方は以下の練習問題の解法を参照。

template<typename T>
vector<T> compress(vector<T> &v1, vector<T> &v2){
  vector<T> vals;
  int n = v1.size();
  for(int i=0; i<n; i++){
    for(int j=0; j<=1; j++){
      vals.push_back(v1[i]+j);
      vals.push_back(v2[i]+j);
    }
  }
  sort(vals.begin(), vals.end());
  vals.erase(unique(vals.begin(), vals.end()), vals.end());
  for(int i=0; i<n; i++){
    v1[i] = lower_bound(vals.begin(), vals.end(), v1[i]) - vals.begin();
    v2[i] = lower_bound(vals.begin(), vals.end(), v2[i]) - vals.begin();
  }
  return vals;
}

練習問題

座標圧縮が必要になる問題では、座標圧縮を行い、imos法により圧縮後の座標をグリッドにマッピングするという処理を書くことが多い。(そのため、ここでは多次元imos法の書き方を知っていることを前提とする)

そして以下の様なことを問う問題が多い。

  1. 座標平面と長方形が与えられて、長方形が重なる領域の面積の総和はいくつか?(三次元の問題もある)
  2. 座標平面と線分が与えられて、線分で区切られた領域の数はいくつか?
二次元座標圧縮

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DSL_4_A&lang=ja

問題形式1の問題で、圧縮した座標をグリッドにマッピングし、そのグリッドの1マス1マスの面積を求め、その総和を求めることで解を求めている。

signed main()
{
  int n; cin>>n;
  vector<long long> x1(n), y1(n), x2(n), y2(n);
  for(int i=0; i<n; i++) cin>>x1[i]>>y1[i]>>x2[i]>>y2[i];
  
  vector<long long> X = compress(x1, x2);
  vector<long long> Y = compress(y1, y2);
  
  int h = Y.size(), w = X.size();
  vector<vector<int>> G(h, vector<int>(w));
  
  // 二次元imos
  for(int i=0; i<n; i++){
    G[y1[i]][x1[i]]++;
    G[y2[i]][x2[i]]++;
    G[y1[i]][x2[i]]--;
    G[y2[i]][x1[i]]--;
  }
  for(int i=0; i<h-1; i++){
    for(int j=0; j<w; j++){
      G[i+1][j] += G[i][j];
    }
  }
  for(int i=0; i<h; i++){
    for(int j=0; j<w-1; j++){
      G[i][j+1] += G[i][j];
    }
  }
  
  // グリッドにおける各セルの面積の総和を求めることで解を得る
  long long ans=0;
  for(int i=0; i<h; i++){
    for(int j=0; j<w; j++){
      if(G[i][j]>0){
        ans += (Y[i+1]-Y[i]) * (X[j+1]-X[j]);
      }
    }
  }
  cout<<ans<<endl;
}

atcoder.jp

問題形式2の問題。座標圧縮したグリッドに対してBFSを行い、その回数が解となる。

この問題では圧縮する座標にグリッドの両端を加えるため、引数に head, tail を追加している。

#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i)
template<class t> using vc=vector<t>;

template <typename T>
vector<T> compress(vector<T> &v1, vector<T> &v2, T head, T tail) {
  vector<T> vals;
  int n = v1.size();
  for (int i = 0; i < n; i++) {
    vals.push_back(v1[i]);
    vals.push_back(v2[i]);
    if(v2[i] + 1 <= tail) vals.push_back(v2[i] + 1);
  }
  vals.push_back(head);
  vals.push_back(tail);
  sort(vals.begin(), vals.end());
  vals.erase(unique(vals.begin(), vals.end()), vals.end());
  for (int i = 0; i < n; i++) {
    v1[i] = lower_bound(vals.begin(), vals.end(), v1[i]) - vals.begin();
    v2[i] = lower_bound(vals.begin(), vals.end(), v2[i]) - vals.begin();
  }
  return vals;
}

const int dy[] = {0, 1, 0, -1};
const int dx[] = {1, 0, -1, 0};

inline bool inside(int y, int x, int H, int W) {
    return (y >= 0 && x >= 0 && y < H && x < W);
}

signed main()
{
  cin.tie( 0 ); ios::sync_with_stdio( false );
  
  int H,W; cin>>H>>W;
  int n; cin>>n;
  vc<int> x1(n),y1(n),x2(n),y2(n);
  rep(i,n) cin>>x1[i]>>y1[i]>>x2[i]>>y2[i];
  
  auto X=compress(x1,x2,0,H);
  auto Y=compress(y1,y2,0,W);
  
  int h=X.size(), w=Y.size();
  vc<vc<int>> G(h, vc<int>(w));
  
  // 二次元imos
  rep(i,n){
    G[x1[i]][y1[i]]++;
    G[x1[i]][y2[i]]--;
    G[x2[i]][y1[i]]--;
    G[x2[i]][y2[i]]++;
  }
  rep(i,h){
    rep(j,w-1){
      G[i][j+1]+=G[i][j];
    }
  }
  rep(j,w){
    rep(i,h-1){
      G[i+1][j]+=G[i][j];
    }
  }

  h--, w--;
  
  // BFS
  int ans=0;
  rep(i,h){
    rep(j,w){
      if(G[i][j]==0){
        ans++;
        queue<pair<int,int>> q;
        q.emplace(i,j);
        G[i][j]=1;
        while(!q.empty()){
          int x,y; tie(x,y)=q.front(); q.pop();
          rep(k,4){
            int nx=x+dx[k], ny=y+dy[k];
            if(inside(nx,ny,h,w) and G[nx][ny]==0){
              G[nx][ny]=1;
              q.emplace(nx,ny);
            }
          }
        }
      }
    }
  }
  cout<<ans<<endl;
}
三次元座標圧縮

atcoder.jp

これも問題形式1の問題で三次元に拡張したもの。(三次元imos法の書き方の練習にもなる)

#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i)
using ll = long long;
template<class t> using vc=vector<t>;

signed main()
{
  cin.tie( 0 ); ios::sync_with_stdio( false );
  
  int n,K; cin>>n>>K;
  vc<ll> x1(n), y1(n), z1(n), x2(n), y2(n), z2(n);
  rep(i,n){
    cin>>x1[i]>>y1[i]>>z1[i]>>x2[i]>>y2[i]>>z2[i];
  }
  vc<ll> X=compress(x1,x2);
  vc<ll> Y=compress(y1,y2);
  vc<ll> Z=compress(z1,z2);
  int h=X.size(), w=Y.size(), d=Z.size();
  vc<vc<vc<ll>>> G(h, vc<vc<ll>>(w, vc<ll>(d)));
  
  // 三次元imos
  rep(i,n){
    G[x1[i]][y1[i]][z1[i]]++;
    G[x1[i]][y2[i]][z1[i]]--;
    G[x1[i]][y1[i]][z2[i]]--;
    G[x2[i]][y1[i]][z1[i]]--;
    G[x1[i]][y2[i]][z2[i]]++;
    G[x2[i]][y1[i]][z2[i]]++;
    G[x2[i]][y2[i]][z1[i]]++;
    G[x2[i]][y2[i]][z2[i]]--;
  }
  rep(i,h){
    rep(j,w){
      rep(k,d-1){
        G[i][j][k+1]+=G[i][j][k];
      }
    }
  }
  rep(i,h){
    rep(k,d){
      rep(j,w-1){
        G[i][j+1][k]+=G[i][j][k];
      }
    }
  }
  rep(k,d){
    rep(j,w){
      rep(i,h-1){
        G[i+1][j][k]+=G[i][j][k];
      }
    }
  }
  
  ll ans=0;
  rep(i,h){
    rep(j,w){
      rep(k,d){
        if(G[i][j][k]>=K){
          ans+=(X[i+1]-X[i])*(Y[j+1]-Y[j])*(Z[k+1]-Z[k]);
        }
      }
    }
  }
  cout<<ans<<endl;
}

参考

参考1)

www.slideshare.net

参考2)実装を参考にさせていただきました。

algo-logic.info

Dijkstra法による最小費用流を求める実装(Primal-Dual)とポテンシャルに関する個人的なメモ

前提

基本的な考えは、Bellman-Ford法による最小費用流の求め方と同じである。

以下の実装では高速化のためBellman-Ford法の実装個所をDijkstra法に置き換えており、同時に負の辺が存在するグラフに対してDijkstra法を適用するため、ポテンシャルを導入している。

計算量: O(flow|E|log|V|)(最大流flow、頂点数V、辺数Eとする)

ポテンシャルについて

※ 以下の説明は、niuezさんのブログ(参考1)を参考に自分の理解を整理するために書いたものです。

最短距離の数式化

単一始点最短経路を求め済みのグラフにおいて、頂点iの最短距離をd[i]、頂点(i,j)の辺のコストをdist(i,j)とすると、辺が存在する頂点(i, j)間には以下のことが成り立つ。(これが最短距離を数式として落とし込んだものになる)

  • d[s] = 0 とする。(sは始点)
  • d[i] + dist(i, j) >= d[j] (iが始点に近く、jが終点に近い頂点とする)

ポテンシャルを導入し、グラフ(辺のコスト)を再構築する

各頂点にポテンシャルhを決めて、辺の距離を dist'(i, j) = dist(i, j) + h[i] - h[j] に置き換えたグラフG'を考える。 G'における各頂点の最短距離をd'とする。G'にも、最短距離の数式が成り立ち、更に2行目の様に式変換ができる。

d'[i] + dist'(i, j) >= d'[j]
⇔ d'[i] + h[i] + dist(i, j) >= d'[j] + h[j]

2行目から、d[i] = d'[i] + h[i] と見れば、d[i]は最短距離の式の形をしており、d[s] = d'[s] + h[s] より各頂点iについて d'[i] + h[i] - h[s] を計算すれば元の最短距離が求められるとわかる。

各辺を非負とするポテンシャルの求め方

ここまでで、ポテンシャルを導入して元のグラフの最短距離を求められることが分かった。次に、Dijkstra法を適用できるようにするためdist'を非負にしたいので、ポテンシャルに適切な値を設定することを考える。

dist(i, j) + h[i] - h[j] >=0
⇔ h[i] + dist(i, j) >= h[j]

2行目から最短距離の式の形をしているため、始点からの最短距離がdist'を非負にするためのポテンシャルであるとわかる。つまり1度Bellman-Ford法等により、各頂点の最短距離を求めておけば、以降その値をポテンシャルとすることでDijkstra法によっても最短距離を求めることができるようになる。

まとめ

  1. 一度Bellman-Ford法等によって最短距離を求められればポテンシャルも求まる
  2. ポテンシャルから非負の辺のみからなるグラフG'を構築でき、これに対してDijkstra法を適用できる
  3. 2で求めたG'における各頂点の最短距離d'から始点sのポテンシャルを減算すれば元のグラフにおける最短距離を求められる

ちなみに以下の最小費用流を求める実装においては、Dijkstra法によって1度目の最短距離を求めている。これは、初期状態のグラフにおいて移動可能な経路上に負の辺が存在しないことを前提としているため。

参考

参考1)ポテンシャルについて

niuez.hatenablog.com

参考2)負閉路が無い事とポテンシャルが存在することは同値であることの説明

検証用問題

onlinejudge.u-aizu.ac.jp

実装

template <typename capacity_t, typename cost_t>
class PrimalDual {
private:
  const cost_t INF = 1e9;
  struct Edge {
    int to, rev;
    capacity_t cap;
    cost_t cost;
    Edge(int to, int rev, capacity_t cap, cost_t cost) : to(to), rev(rev), cap(cap), cost(cost) {}
  };
  vector<vector<Edge>> G;
  vector<int> prev_v, prev_e;
  Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }

public:
  PrimalDual(int n) {
    G.resize(n);
    prev_v.resize(n);
    prev_e.resize(n);
  }
  void add_edge(int from, int to, capacity_t cap, cost_t cost) {
    G[from].push_back(Edge(to, (int)G[to].size(), cap, cost));
    G[to].push_back(Edge(from, (int)G[from].size() - 1, 0, -cost));
  }

  cost_t get_min_cost_flow(int s, int t, capacity_t flow) {
    int n = G.size();
    vector<cost_t> h(n); // ポテンシャル
    cost_t res = 0;
    while (flow > 0) {
      vector<cost_t> dist(n, INF);
      dist[s] = 0;
      priority_queue<pair<cost_t, int>, vector<pair<cost_t, int>>, greater<pair<cost_t, int>>> que;
      que.emplace(0, s);
      while (!que.empty()) {
        cost_t c;
        int v;
        tie(c, v) = que.top();
        que.pop();
        if (c > dist[v]) continue;
        for (int i = 0; i < G[v].size(); i++) {
          auto &edge = G[v][i];
          if (edge.cap > 0 and dist[edge.to] > dist[v] + edge.cost + h[v] - h[edge.to]) {
            dist[edge.to] = dist[v] + edge.cost + h[v] - h[edge.to];
            prev_v[edge.to] = v;
            prev_e[edge.to] = i;
            que.emplace(dist[edge.to], edge.to);
          }
        }
      }
      if (dist[t] == INF) return -1;
      for (int v = 0; v < n; v++) h[v] += dist[v];
      capacity_t d = flow;
      for (int v = t; v != s; v = prev_v[v]) {
        d = min(d, G[prev_v[v]][prev_e[v]].cap);
      }
      flow -= d;
      res += d * h[t];
      for (int v = t; v != s; v = prev_v[v]) {
        G[prev_v[v]][prev_e[v]].cap -= d;
        get_rev(G[prev_v[v]][prev_e[v]]).cap += d;
      }
    }
    return res;
  }
};