Baby Stepsなブログ

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

Bellman-Ford法による最小費用流を求める実装

前提

最小費用流問題のグラフには、最大流問題で与えられる辺の情報に加えて辺のコストが付与される。

以下の実装では、このコストに対してS-T最短路を1つ見つけて目いっぱい流す、その経路が限界になったら次の最短路を見つけて目いっぱい流して、という事を求める流量に達するまで反復して行うことで最小費用流を求めている。

最小費用流問題を求める際にも残余グラフを構築するが、この際、逆辺のコストは元のコストを-1倍した値を設定することになる。

故に残余グラフには負の辺を含むことになるため、そのままではDijkstra法を使うことができないため、最短路を求めるためにBellman-Ford法を利用している。(ポテンシャルを導入することで最短路を求める工程をDijkstra法に置き換える方法もある => 続き

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

参考

P.48 以降参照

http://hos.ac/slides/20150319_flow.pdf

検証用問題

onlinejudge.u-aizu.ac.jp

実装

template <typename capacity_t, typename cost_t>
class MinCostFlow {
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:
  MinCostFlow(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();
    cost_t res = 0;
    while (flow > 0) {
      vector<cost_t> dist(n, INF);
      dist[s] = 0;
      bool update = true;
      while (update) {
        update = false;
        for (int v = 0; v < n; v++) {
          if (dist[v] == INF) 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) {
              dist[edge.to] = dist[v] + edge.cost;
              prev_v[edge.to] = v;
              prev_e[edge.to] = i;
              update = true;
            }
          }
        }
      }
      if (dist[t] == INF) return -1;
      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 * dist[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;
  }
};

最大流問題とDinic法に関する個人的なメモ

前提

Dinic法を理解するには、先にFord-Fulkerson法を知っていると早い。

pione.hatenablog.com

Dinic法の基本的な動作はFord-Fulkerson法と同じで、残余グラフを構築して増加路の辺、逆辺の capacity を更新していき、更新できるパスがなくなるまでそれを繰り返している。

Dinic法ではこれを効率化するためにBFSとDFSを併用している。(DFSとBFSの実行順序は以下の参考1のフローチャートがわかりやすい。)

BFS

始点sから各ノードへの最短距離levelを求める。level配列は-1で初期化し、BFSは移動可能な辺(つまり level[v]<0 && capacity>0 である辺)のみを使って遷移し、level配列を更新していく。 level[t]==-1 であるなら、増加路が存在しないことになるのでその時点で更新を終了する。

DFS

次に、求めたlevelを活用しつつDFSを行う。 グラフに増加路が存在するとき、その増加路の頂点のどの隣り合う頂点間には、level[from]<level[to] && capacity>0 が成り立つはずなので、これを満たす頂点間のみの遷移を見ればよい。

ある1つの増加路の流量を求められたら、Ford-Fulkerson法と同じ要領で残余グラフの辺、逆辺の capacity を更新する。 この際、調べた辺を記録しておく。(以下の実装におけるitr配列)

再度DFSにより増加路を探す際、itr配列を使えば、既にその辺から先に増加路は無いと分かっている辺をスキップすることができるため効率化ができる。(ここがFord-Fulkerson法と比較したDinic法の優位な点だと思う。)

計算量: O(|V|^2|E|)(頂点数をV、辺数をEとする)

実際にはこの計算量よりも早い結果になることが多い。(参考2の記事を参照)

参考

参考1)

以下の記事のフローチャートはとても参考になる。

vartkw.hatenablog.com

参考2)

計算量に関する記事。

topcoder-g-hatena-ne-jp.jag-icpc.org

検証用問題

onlinejudge.u-aizu.ac.jp

実装

template <class T>
class Dinic {
private:
  const int INF = 1e9;
  vector<int> level, itr;

  struct Edge {
    int to, rev; T cap;
    Edge(int to, int rev, T cap) : to(to), rev(rev), cap(cap){};
  };

  vector<vector<Edge>> G;

  Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }

  void bfs(int s) {
    level.assign(G.size(), -1);
    level[s] = 0;
    queue<int> q;
    q.push(s);
    while (!q.empty()) {
      auto v = q.front();
      q.pop();
      for (auto &e : G[v]) {
        if (level[e.to] < 0 and e.cap > 0) {
          level[e.to] = level[v] + 1;
          q.push(e.to);
        }
      }
    }
  }

  T dfs(int v, int t, T flow) {
    if (v == t) return flow;
    for (int &i = itr[v]; i < G[v].size(); i++) {
      auto &edge = G[v][i];
      if (level[v] < level[edge.to] and edge.cap > 0) {
        auto f = dfs(edge.to, t, min(flow, edge.cap));
        if (f > 0) {
          edge.cap -= f;
          get_rev(edge).cap += f;
          return f;
        }
      }
    }
    return 0;
  }

public:
  Dinic(int n) { G.resize(n); }

  void add_edge(int from, int to, T cap) {
    G[from].push_back(Edge(to, (int)G[to].size(), cap));
    G[to].push_back(Edge(from, (int)G[from].size() - 1, 0));
  }

  T get_max_flow(int s, int t) {
    int n = G.size();
    T res = 0;
    while (true) {
      itr.assign(n, 0);
      bfs(s);
      if (level[t] < 0) break;
      while (true) {
        T flow = dfs(s, t, INF);
        if (flow > 0) res += flow;
        else break;
      }
    }
    return res;
  }
};

最大流問題とFord-Fulkerson法に関する個人的なメモ

前提

最大流問題とは、辺に capacity を持つグラフが与えられて、始点から目いっぱい水を流した時に終点に流せる最大量を求める問題。

グラフの最大流は、そのグラフのS-T最小カットに対応する。

計算量: O(|flow||E|) (最大フローをflow、辺数をEとする)

|E| の部分は一度のDFSを意味し、一度のDFSで少なくとも1以上流量が改善され、流量の改善は高々そのグラフの最大流までしかないことから計算量に最大流 flow が登場する。 実際には1ずつの更新が flow 回行われるということはほとんどないので |flow| より計算量は良くなる。

Ford-Fulkerson法では、与えられたグラフに対して以下の様な規則で逆辺をも持つ新しいグラフを構築し、このグラフに対して貪欲的に値を更新していくことで解を求めている。(この逆辺を張ったグラフを残余グラフと呼ぶ。)残余グラフではなく、元のグラフに対して貪欲法で解を求めようとすると正しい解が求まらないことは、以下に載せた 参考2 のスライドに図解がある。

与えられた元のグラフのある辺に対して、同じ方向の辺はその辺に流すことができる残容量、逆方向の辺は既にその辺に流した流量、として双方向に辺を張ったグラフを構築する。すると、逆辺の値を減らすことは今流している流量を減らすこと、元の辺の値を減らすことは残容量を減らす(つまり辺に追加で流す)ことに対応する。場合によっては、ある辺の流量を減らす(逆辺の値を減らす)ことで全体の流量が改善する場合があり(参考1 P23 参照)、そのため流量を増やす、減らすの双方向の操作ができる残余グラフである必要がある。

ここで、同じ辺が流量を増やされたり減らされたり何度も更新される可能性があるため、無限ループしないかという不安がある。 しかし、更新を続けることでいずれS-T最小カットが定まる(つまりS-Tの増加路がなくなる)ため、その時点で更新が終了するため心配はない。

二部グラフの最大マッチング

最大流の有名な応用例として、二部グラフの最大マッチングを求めるものがある。

元の二部グラフに最大流を流すための source と sink 用の2つの頂点を追加することを考える。

二部グラフの頂点部分集合、X、Yのどちらかに対して source からの辺を張り、もう片方から sink への辺を張る。すると s-t最大流が求める最大マッチングの値になる。(参考4に最大流が二部グラフのマッチングになることの図解がある)

参考

参考1)残余グラフの更新の様子が図解してあり非常にわかりやすいです。

http://hos.ac/slides/20150319_flow.pdf

参考2)P.13 に間違った貪欲法では求まらない例の図が参考になります。

www.slideshare.net

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

algo-logic.info

参考4)二部グラフの最大マッチングについて

qiita.com

検証用問題

最大流

onlinejudge.u-aizu.ac.jp

二部グラフの最大マッチング

onlinejudge.u-aizu.ac.jp

実装

template <class T>
class FordFulkerson {
private:
  const T INF = 1e9;
  
  struct Edge {
    int to, from, rev; T cap;
    Edge(int from, int to, int rev, T cap) : to(to), from(from), rev(rev), cap(cap) {}
  };
  
  vector<vector<Edge>> G;
  
  Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }
  
  vector<bool> used;
  
  T dfs(int v, int t, T flow) {
    if (v == t) return flow;
    used[v] = true;
    for (auto &e : G[v]) {
      if (used[e.to] or e.cap <= 0) continue;
      T d = dfs(e.to, t, min(flow, e.cap));
      if (d > 0) {
        e.cap -= d;
        get_rev(e).cap += d;
        return d;
      }
    }
    return 0;
  }

public:
  FordFulkerson(int n) { G.resize(n); }
  
  void add_edge(int from, int to, T cap) {
    G[from].push_back(Edge(from, to, (int)G[to].size(), cap));
    G[to].push_back(Edge(to, from, (int)G[from].size() - 1, 0));
  }
  
  T get_max_flow(int s, int t) {
    int n = G.size();
    T res = 0;
    while (true) {
      used.assign(n, false);
      T flow = dfs(s, t, INF);
      if (flow > 0) {
        res += flow;
      } else {
        return res;
      }
    }
    return 0;
  }
};

Codeforces Round #703 (Div. 2) 参加記

A(500) C1(750) E(2250) の3完。良問が多い回だった。

A - Shifting Stacks (500点)

英文読解の問題。

前の要素を自身より後ろの任意の場所に移動させられるという条件で、狭義単調増加を作れるかという問題。

B - Eastern Exhibition (1000点)

本番中にこれが解けないのはだめ。

いつも通り x, y を独立に考える。

すると数直線上にある1点を取って、他点との差の絶対値の和を最小にする問題になる。(中央値の性質の有名問題)ここまでこれば奇数なら中央値、偶数なら2つの中央値の間の任意の整数値が1方向の取り得る座標の値だとわかる。

C1 - Guessing the Greatest (easy version) (750点)

平方分割。

まず区間全体のsecond maxのindexを取得する。次にそのindexを含むように区間を狭めてクエリを投げてみて、second maxのindexが同じなら狭めた区間に最大値が存在し、異なるのならその区間に最大値はないことがわかるので対象区間が絞れる。これを繰り返せば解が求まる。

105 < 220 なので、同じ実装でC2も通ると考えたがWAだった。(どうしてこんなことに)

E - Paired Payment (2250点)

単純にグラフを再構築してダイクストラをする。これが初期値2250点はお得。

まとめ

#include <atcoder/all>を消し忘れて開幕CEにならないようにする。

JOI 本選 2012 E - JOI 国のお祭り事情(Festivals in JOI Kingdom) 解説

問題

onlinejudge.u-aizu.ac.jp

前提

解説

ダイクストラ法、クラスカル法(的な発想)、Lowest Common Ancestor、3つ基本アルゴリズムの合わせ技で解ける良問。

ざっと他の人の解説を見てみたところ、並行二分探索による解説が多かったので、3つのアルゴリズムを使う解説を載せておく。

(1)ダイクストラ

まず最初に各頂点に、自身から最も近いお祭りまでの距離を記録する => ダイクストラ法。(辺にコストがあるためBFSだとTLEになることに注意する)

(2)クラスカル

(1)を元に、辺のコストを両端の2頂点のお祭りまでの距離のminとするグラフを作る。このグラフから、ある経路を通った時の祭りまでの距離の最小は通った辺のコストの最小に等しくなる。

このグラフを更に辺のコストの大きい順に取っていくことで全域木を作ることを考える。(ここがクラスカル法的な発想)

全域木になった時点で任意の2頂点の行き来は可能である。そして、この過程で全域木を作り残った辺を追加しても、ある頂点間のお祭りまでの距離が改善すること(大きくなること)がない。よって、この全域木における2頂点間のコストの最小がお祭りまでの距離になる。

(3)Lowest Common Ancestor

最後に、2頂点間の経路は(s, t)のLCAを求めて、s, t それぞれからLCAに向けて頂点を見ていけば無駄なく探索を行える。これはダブリングのテーブルがあれば行える。

(3)における2頂点間のお祭りまでの距離の最小値を求める処理がボトルネックで、O(n * q) となるが、制約から N ≤ 5000, Q ≤ 5000 のいずれかを満たすため間に合う。

実装

#include <bits/stdc++.h>
using namespace std;

#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i)

template <class T> inline bool chmax(T &a, T b) {
  if (a < b) {
    a = b;
    return 1;
  }
  return 0;
}
template <class T> inline bool chmin(T &a, T b) {
  if (a > b) {
    a = b;
    return 1;
  }
  return 0;
}

const int INF = 1e9 + 10;

class UnionFind {
private:
  vector<int> parent;

public:
  UnionFind(int n) { parent = vector<int>(n, -1); }

  int root(int a) {
    if (parent[a] < 0)
      return a;
    return parent[a] = root(parent[a]);
  }

  int size(int a) { return -parent[root(a)]; }

  bool connect(int a, int b) {
    a = root(a);
    b = root(b);
    if (a == b) {
      return false;
    }
    if (size(a) < size(b))
      swap(a, b);
    parent[a] += parent[b];
    parent[b] = a;
    return true;
  }

  bool is_same(int a, int b) { return root(a) == root(b); }
};

class LCA {
private:
  int root;
  int k;
  vector<vector<int>> dp;
  vector<int> depth;
  vector<int> dist;

public:
  LCA(const vector<vector<int>> &_G, vector<int> _dist, const int _root = 0) {
    int n = _G.size();
    dist = _dist;
    root = _root;
    k = 1;
    int nibeki = 2;
    while (nibeki < n) {
      nibeki <<= 1;
      k++;
    }
    dp = vector<vector<int>>(k + 1, vector<int>(n, -1));
    depth.resize(n);
    function<void(int, int)> _dfs = [&](int v, int p) {
      dp[0][v] = p;
      for (auto nv : _G[v]) {
        if (nv == p)
          continue;
        depth[nv] = depth[v] + 1;
        _dfs(nv, v);
      }
    };
    _dfs(root, -1);
    for (int i = 0; i < k; i++) {
      for (int j = 0; j < n; j++) {
        if (dp[i][j] == -1)
          continue;
        dp[i + 1][j] = dp[i][dp[i][j]];
      }
    }
  }

  /// get LCA
  int get(int u, int v) {
    if (depth[u] < depth[v])
      swap(u, v);
    if (depth[u] != depth[v]) {
      long long d = depth[u] - depth[v];
      for (int i = 0; i < k; i++)
        if ((d >> i) & 1)
          u = dp[i][u];
    }
    if (u == v)
      return u;

    for (int i = k; i >= 0; i--) {
      if (dp[i][u] != dp[i][v]) {
        u = dp[i][u], v = dp[i][v];
      }
    }
    return dp[0][u];
  }

  int get_min_cost(int u, int v) {
    int lca = get(u, v);
    int min_cost = dist[lca];
    while (u != lca) {
      chmin(min_cost, dist[u]);
      u = dp[0][u];
    }
    while (v != lca) {
      chmin(min_cost, dist[v]);
      v = dp[0][v];
    }
    return min_cost;
  }
};

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

  // input
  int n, m, k, q;
  cin >> n >> m >> k >> q;
  vector<int> a(k);
  vector<vector<pair<int, int>>> G(n);
  rep(i, m) {
    int a, b, l;
    cin >> a >> b >> l;
    a--, b--;
    G[a].emplace_back(b, l);
    G[b].emplace_back(a, l);
  }

  // dijkstra
  vector<int> dist(n, INF);
  priority_queue<pair<int, int>, vector<pair<int, int>>, greater<pair<int, int>>> que;
  rep(i, k) {
    int f;
    cin >> f;
    f--;
    dist[f] = 0;
    que.emplace(0, f);
  }

  while (!que.empty()) {
    int c, v;
    tie(c, v) = que.top();
    que.pop();
    if (c > dist[v])
      continue;
    for (auto p : G[v]) {
      int nv, nc;
      tie(nv, nc) = p;
      if (chmin(dist[nv], dist[v] + nc)) {
        que.emplace(dist[nv], nv);
      }
    }
  }

  // kruskal
  vector<tuple<int, int, int>> edge;
  rep(v, n) for (auto p : G[v]) {
    int nv, c;
    tie(nv, c) = p;
    edge.emplace_back(min(dist[v], dist[nv]), v, nv);
  }
  sort(edge.rbegin(), edge.rend());

  auto uf = UnionFind(n);
  vector<vector<int>> tree(n);
  for (auto e : edge) {
    int c, v, nv;
    tie(c, v, nv) = e;
    if (uf.is_same(v, nv))
      continue;
    uf.connect(v, nv);
    tree[v].push_back(nv);
    tree[nv].push_back(v);
  }

  // lca
  auto lca = LCA(tree, dist);
  while (q--) {
    int s, t;
    cin >> s >> t;
    s--, t--;
    cout << lca.get_min_cost(s, t) << endl;
  }
}

ABC 185 E - Sequence Matching 解説

問題

atcoder.jp

前提

解説

最長共通部分列問題(Longest Common Subsequence)が解法の基盤になっていると思った。

LCSのDPでは、複数の列における部分列の長さの最大値を管理するが、この問題では不一致の部分を残した方がスコアが良い場合があるので、それを考える必要がある。

例) A={1,2,9,9,3}, B={1,2,8,8,3} LCSは{1,2,3} で、LCS以外を全て除外する場合のスコアは4だが、そのまま残しておけばスコアは2になる

よって以下の様にDPテーブルを定義して、ある要素を除外する場合と残す場合それぞれの遷移を考えて更新していく。

dp[i][j] := 数列Aの先頭i、数列Bの先頭jまで見た時のスコアの最小値

まずi=0 or j=0 という場合は、長さを一致させるために片方の文字列は全消しするしかないので、この場合は、dp[i][j] = max(i, j)

これがテーブルの初期化に対応する。

i, jがいずれも1文字以上の場合の遷移は2通りに分けて考える。

  • A[i] == B[j] の場合

dp[i+1][j+1] = dp[i][j] + 0 // i-1, j-1までの部分列の末尾に同じ文字を追加してもスコアは変わらない

  • A[i] != B[j] の場合

ここで更に末尾の要素を残すか、残さないか2通りの遷移を考える必要がある。

  1. dp[i+1][j+1] = dp[i][j] + 1 // i-1, j-1 までの部分列の末尾にA[i], B[i]を残す。スコア+1
  2. dp[i+1][j+1] = min(dp[i][j + 1], dp[i + 1][j]) + 1 // 少なくとも片方の末尾を除外する場合。A側を除外する場合、B側を除外する場合の小さい方のスコアに+1

実装

#include <bits/stdc++.h>
using namespace std;

#define rep(i, n) for (long long i = (long long)(0); i < (long long)(n); ++i)
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }

const int INF = 1e9;

signed main() {
  int n, m;
  cin >> n >> m;
  vector<int> a(n), b(m);
  rep(i, n) cin >> a[i];
  rep(i, m) cin >> b[i];

  vector<vector<int>> dp(n + 1, vector<int>(m + 1, INF));
  rep(i, n + 1) dp[i][0] = i;
  rep(i, m + 1) dp[0][i] = i;
  rep(i, n) {
    rep(j, m) {
      if (a[i] == b[j]) {
        chmin(dp[i + 1][j + 1], dp[i][j]);
      } else {
        chmin(dp[i + 1][j + 1], dp[i][j] + 1);
        chmin(dp[i + 1][j + 1], min(dp[i][j + 1], dp[i + 1][j]) + 1);
      }
    }
  }
  cout << dp[n][m] << endl;
}

yukicoder バレンタインコンテスト2021 B Giri 解説

yukicoder.me

問題

1~nからn-1個を選んで作れる最小公倍数のうち、その最小値を求め、998244353で割った余りを求めよ。

前提

  • エラトステネスの篩

pione.hatenablog.com

解説

まず除外すべき1つの数は自明に、1~nの中の最大の素数である。

次にn-1個の数の最小公倍数を求め方について、結果が非常に大きな値になる可能性があるためlcm関数を使って求めることはできない。

この手の、n個全ての数の最小公倍数を求めたいが、結果が非常に大きくなる可能性があるからmodを取る問題の求め方は決まっていて、各要素を素因数分解すればよい。

各要素を素因数分解し、各素因数の種類とその最大個数が得られれば、(素因数)^(その個数)の全ての積が最小公倍数になる。

後は、上記をmodを取りながら計算すればよい。

高速素因数分解を行うためにエラトステネスの篩を行うため、計算量はO(n loglog n)となる。

実装

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 998244353;

/// エラトステネスの篩
struct Eratosthenes {
  vector<bool> is_prime;
  /// nまでの素数一覧
  vector<int> primes;
  /// 整数xの最小の素因数
  vector<int> min_prime_factor;

  Eratosthenes(const int n) {
    min_prime_factor.resize(n + 1);
    /// nまでの素数一覧
    is_prime.assign(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    min_prime_factor[0] = min_prime_factor[1] = 1;
    for (int i = 2; i <= n; i++) {
      if (is_prime[i]) {
        for (int j = 2 * i; j <= n; j += i) {
          is_prime[j] = false;
          if (min_prime_factor[j] == 0)
            min_prime_factor[j] = i;
        }
        primes.push_back(i);
        min_prime_factor[i] = i;
      }
    }
  }

  /**
   * xの素因数の種類を取得
   * @param x n以下の数値を指定する必要がある
   */
  set<int> get_prime_factors(int x) {
    set<int> s;
    while (x > 1) {
      s.insert(min_prime_factor[x]);
      x /= min_prime_factor[x];
    }
    return s;
  }

  /**
   * xの素因数とその数を取得
   * @param x n以下の数値を指定する必要がある
   */
  map<int, int> get_prime_factor_and_cnt(int x) {
    map<int, int> mp;
    while (x > 1) {
      mp[min_prime_factor[x]]++;
      x /= min_prime_factor[x];
    }
    return mp;
  }
};

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

  int n;
  cin >> n;
  auto e = Eratosthenes(n);
  auto ma = e.primes.back();

  map<int, int> mp;
  for (int i = 1; i <= n; i++) {
    if (ma == i)
      continue;
    auto tmp = e.get_prime_factor_and_cnt(i);
    for (auto e : tmp) {
      mp[e.first] = max(mp[e.first], e.second);
    }
  }
  long long ans = 1;
  for (auto e : mp) {
    for (int i = 0; i < e.second; i++) {
      (ans *= e.first) %= MOD;
    }
  }
  cout << ans << endl;
}