Baby Stepsなブログ

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

DFSによるオイラーツアー順序を求める実装

高難度帯の部分問題として出てくるので実装を置いておく。

実装

以下の実装は、訪れた順に頂点を採番していく。

/// オイラーツアー
struct EulerianTrail{
  vector<int> order;
  
  EulerianTrail(vector<vector<int>>& G, int root=0){
    int n=G.size();
    order.resize(n);
    int cur=0;
    function<void(int,int)> dfs=[&](int v, int p){
      order[v]=cur;
      cur++;
      for(auto nv: G[v]){
        if(nv==p) continue;
        dfs(nv, v);
      }
    };
    dfs(root, -1);
  }
};

トポロジカルソートの実装をライブラリ化した

過去問埋めしてて、トポロジカルソートで解く問題に当たって面白かった↓。

atcoder.jp

良い機会なので、改めて実装についてまとめておこうと思う。

前提

トポロジカルソートとは、閉路無し有向グラフ(DAG)において、どの頂点もその出力辺の先の頂点より前にくるように並べること。

ようするに、そのグラフの頂点を「実行順序があるタスク」と見た時、先に実行すべきタスクの順に並べることを言う。

  1. 各頂点の入次数をカウントしておく
  2. 閉路無し有向グラフなら必ず入次数0の頂点があるので、queue に push する
  3. 以下の処理を queue が空になるまで繰り返す
    1. queue に入っている頂点を順に取り出し(pop)、今見ている頂点から移動可能な頂点の入次数を -1 する
    2. 入次数が0になった頂点は都度 queue に push する

こうして queue に push された頂点の順序がトポロジカル順序になる

検証用問題

onlinejudge.u-aizu.ac.jp

実装

/**
 * トポロジカルソート
 * 頂点数をV、辺数をEとして O(V+E)
 */
vector<int> topological_sort(const vector<vector<int>>& G){
  const int n=G.size();
  vector<int> in_cnt(n); // 各頂点の入次数
  // 各頂点の入次数をカウント
  for(int v=0; v<n; v++){
    for(auto nv: G[v]){
      in_cnt[nv]++;
    }
  }
  // 入次数が0の頂点を始点としてqueueにpush(同時にresにもpush)
  vector<int> res;
  queue<int> q;
  for(int v=0; v<n; v++){
    if(in_cnt[v]==0){
      q.push(v);
      res.push_back(v);
    }
  }
  // queueに入っている頂点から移動可能な頂点の入次数をデクリメントし、入次数が0になったものをqueueにpush
  while(!q.empty()){
    int v=q.front(); q.pop();
    for(auto nv: G[v]){
      in_cnt[nv]--;
      if(in_cnt[nv]==0){
        q.push(nv);
        res.push_back(nv);
      }
    }
  }
  // reverse(res.begin(), res.end()); // 入次数0の頂点を先頭にしたい場合
  return res;
}

個人的なメモ

頂点 i => j (i < j)という様に辺が与えられて、2点間から計算されるスコアの最高値を求めよ、みたいな問題(上に載せたABC 188 E Peddler がまさにそう)はトポロジカルソートで解く。

この手の問題を dfs で解こうとすると大概TLEになる。(たまにやらかす)それはそうで、1つ頂点が複数の親を持つ可能性があるので、各2点間におけるスコアを求めるには少なくとも入次数0である全ての頂点を始点として dfs する必要がある。つまり頂点数をV、辺数をEとすると計算量は、O( V * (V + E) )になる。

与えられるグラフが木であるなら dfs で問題ないが、上の様な条件では木にならないことに注意する。(木なら各頂点の親の個数は1個以下になる。また辺の数が V-1 個になるので、まず辺数の制約を確認する)

エラトステネスの篩の実装

エラトステネスの篩と、それを利用した素因数の高速列挙の実装をライブラリ化した。(どちらかというと後者の実装を残す目的で作った)

前提

高速素因数分解

  • エラトステネスの篩でnまでの整数を素数判定する過程で、各値の最小の素因数をテーブルmin_prime_factorに保存しておく
  • x / min_prime_factor[x]でxの約数が得られるので、値xの最小の素因数を確認すると同時に、xをその約数に上書きする
  • 上記をxが1になるまでの間繰り返すことで素因数を列挙できる
  • ループ回数は高々log x回なので、O(log x)

検証用問題

atcoder.jp

実装

/// エラトステネスの篩
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;
  }
};

Trie木の実装をライブラリ化した

前提

  • Trie木は、複数の単語(文字列)を登録でき、ある文字列(またはそのprefix)が登録済みであるかを高速に検索できるデータ構造
  • 基本的な機能はシンプルに、insert / searchの2つだけ
  • 実態は有向木
  • 基本的な動作は、単語のprefixが同じならNodeを共有し、異なる文字が現れた際に新規にNodeを追加していく
    • prefixを共有しながら枝分かれしていくイメージ
  • 各Nodeは3つの基本情報を持つ
    1. next: 次のNodeへのポインタ
      • 扱いたい文字種の数char_size(例: a-zなら26)だけ、ポインタを持つ必要がある
    2. commom: そのノードがいくつの単語で共有されているか
      • この情報を持たせることで、そのNode以降にいくつの単語が登録済みかがわかる
      • 故に、Node(0) を確認すればTrie木に登録済みの単語の総数がわかる
    3. accept: そのノードがある単語の末尾文字であるか
      • 単語検索の際にはacceptを確認し、prefix検索の際には確認しない
  • 今Trie木が管理している文字数を知りたい場合は、単にNodeの数をカウントすればよい(Node(0)は文字ではないので、正確にはsize - 1になる)

Trie木のわかりやすい解説

algo-logic.info

検証用問題

leetcode.com

※ 以下の実装を使うには、メソッドのシグネチャとtemplate部分を変更が必要

実装

/**
 * Trie-Tree
 * 
 * @param char_size Trie木で扱う文字種数
 * @param base Trieで扱う文字種のうちの先頭
 */
template<int char_size, int base>
class Trie{
  private:
  struct Node{
    int common;         // この頂点を共有する単語数
    vector<int> next;   // 次の頂点id
    vector<int> accept; // この頂点を終端とする単語idを保持する
    Node(): common(0), next(vector<int>(char_size, -1)) {};
  };
  vector<Node> nodes;
  
  public:
  Trie(): nodes(vector<Node>(1)) {};
  
  void insert(const string& s, int string_id=0){
    nodes[0].common++;
    int last=0;
    for(int i=0; i<s.size(); i++){
      int char_id=(int)(s[i]-base);
      if(nodes[last].next[char_id]==-1){
        nodes.push_back(Node());
        nodes[last].next[char_id]=nodes.size()-1;
        last=nodes.size()-1;
      }else{
        last=nodes[last].next[char_id];
      }
      nodes[last].common++;
    }
    nodes[last].accept.push_back(string_id);
  }
  
  bool search(const string& s){
    int last=0;
    for(int i=0; i<s.size(); i++){
      int char_id=(int)(s[i]-base);
      if(nodes[last].next[char_id]==-1) return false;
      else last=nodes[last].next[char_id];
    }
    return nodes[last].accept.size()>0 ? true : false;
  }
  
  bool search_prefix(const string& s){
    int last=0;
    for(int i=0; i<s.size(); i++){
      int char_id=(int)(s[i]-base);
      if(nodes[last].next[char_id]==-1) return false;
      else last=nodes[last].next[char_id];
    }
    return true;
  }
  
  int count_words(){
    return nodes[0].common;
  }
  
  int size(){
    return nodes.size();
  }
};

三分探索の実装をライブラリ化した

二分探索と比べて、区間の内分点の取り方やら、下に凸か上に凸かやら、地味に考えることが多いので汎用性を目指してライブラリ化してみた。

前提

  • 三分探索とは、ある区間[l, r]において、極値がただ一つだけ存在するとき、その極値の近似値を区間の長さをnとして、およそ O(log n) で求める手法
  • 二分探索が単調グラフにおける最大値、最小値を求めるものに対して、三分探索は凸関数における最大値、最小値を求めるためのもの
  • 三分探索の実装では、以下のことを行っている
    • 区間を等間隔に3つに内分する2点を取り、f(x) が極値とより離れている方を狭める、ということを繰り返す
    • なぜ上記によって区間内に極値が収まるのかは、以下の記事に図解がある

三分探索のわかりやすい解説

juppy.hatenablog.com

実装

  • l, rに、対象の区間を指定する
  • 第3引数fに、対象の関数を与える
  • 上に凸な関数なら、第4引数is_downward_convexに false を指定する
/// 三分探索
double ternary_search(double l, double r, const function<double(double)> f, const bool is_downward_convex=true){
  int loop=500;
  if(is_downward_convex){ // 下に凸
    while(loop--){
      auto mid_l=(l*2+r)/3, mid_r=(l+2*r)/3;
      if(f(mid_l)>=f(mid_r)) l=mid_l;
      else r=mid_r;
    }
  }else{                  // 上に凸
    while(loop--){
      auto mid_l=(l*2+r)/3, mid_r=(l+2*r)/3;
      if(f(mid_l)<=f(mid_r)) l=mid_l;
      else r=mid_r;
    }
  }
  return l;
}

検証用問題

atcoder.jp

呼び出し元

上記問題を例に、以下に実際の使い方を載せる。

double p;

/* ternary_search */

signed main()
{
  cin>>p;
  auto f=[&](double x){
    return x+p/pow(2,x/1.5);
  };
  auto ans_x=ternary_search(0,1e18,f);
  printf("%.15lf\n", f(ans_x));
}

浮動小数を文字列として受け取り、整数型に変換する実装(個人的なメモ)

これは何

最近コンテストで、入力値が小数第4位程度まで与えられ、その入力値を素直にdoubleで受け取って計算すると丸め誤差で WA となるような問題が散見されるようになった。

(直近の問題だと ABC 191 D Circle Lattice Points がそう)

atcoder.jp

この手の問題の対策はいくつかあると思う。

  1. 多倍長浮動小数型がある言語ならそれを使う
  2. long double等、longより精度の高い型 + EPSを使う
  3. 入力値を10p 倍して整数に変換して計算した後、10p で割って最終的な解を求める

以下の実装は、3における入力値を10pするためのもので、丸め誤差が生じないように文字列として入力値を受け取り、それを long long 型で返すものである。(毎回実装考える時間がもったいないので、実装を残しておく)

実装

/**
 * @param s 文字列化した小数値
 * @param p sの数値を 10^p 乗倍する
 * 
 * 10^p乗倍した結果、小数点以下が0になることと、その値がlong long型に収まる必要がある
 */
long long string_double_to_long(const string & s, const int p){
  if(s.find('.')==string::npos){
    long long p=1;
    while(p--) p*=10;
    return stoll(s)*p;
  }else{
    int n=s.size();
    for(int i=0; i<n; i++){
      if(s[i]=='.'){
        int tail_len=n-i-1;
        int loop_cnt=p-tail_len;
        auto t=s.substr(0,i)+s.substr(i+1,n);
        while(loop_cnt--) t.push_back('0');
        return stoll(t);
      }
    }
  }
}

転倒数を求める実装をライブラリ化した

要求される度に、uniqueやらeraseやら、1-indexやらの箇所で微妙に考えてしまうので、この際ライブラリ化しておく。

※ 転倒数には、i < j でかつ、a_i > a_j であるとき転倒とする場合と、a_i >= a_j であるとき転倒とする場合の2種類があり、以下の実装は前者のものに対応する。

前提

  • 転倒数は素直に実装すると、数列の長さをnとして、O(n2)
  • BITを活用することで計算量をO(n log n)に改善できる
  • 先頭から要素を見ていきながら、ループの中ではざっくり2つのことを行っている
    1. 今の要素値 a_ i 以下の個数をBITを用いてカウントする
      • 今の添え字i - ↑の個数 = a_i を末尾に追加することで生じる転倒数
    2. BITに要素 a_i を登録する

実装上での注意点

  • 扱う数列の要素値が正で最大値が小さい場合(105とか)なら、要素の取り得る最大値をBITの長さとして取ればよい。
    • 例えば数列の長さがnで、要素が1~nの連続した値になっている場合がそう。
  • 問題によっては、要素の最大値が1018などになっており、これをそのままBITの要素数とすると空間計算量オーバーとなる場合がある。(以下に載せた検証用問題がそう)
  • その対策として以下の実装では一次元座標圧縮を行っている
    • 今見ている要素が全体で何番目に大きいか(rank)で要素を各管理すればBITのサイズは数列の要素数だけで済む

検証用問題

onlinejudge.u-aizu.ac.jp

実装

※ 要 Binary Index Tree

pione.hatenablog.com

/// 転倒数
long long inversion_number(const vector<long long>& v){
  auto u=v;
  sort(u.begin(), u.end());
  u.erase(unique(u.begin(), u.end()), u.end());
  
  auto bit=BIT<long long>(u.size());
  
  long long res=0;
  for(int i=0; i<v.size(); i++){
    int rank=lower_bound(u.begin(), u.end(), v[i])-u.begin()+1;
    res+=i-bit.sum(rank);
    bit.add(rank);
  }
  return res;
}

座標圧縮の実装に関する個人的なメモ

座標圧縮の実装箇所において、地味に実装バグのトラップとなる箇所があるのでメモしておく。

1. unique 呼び出し前には必ずソート

未ソートの状態でもuniqueは当然呼び出し可能だが、この場合、隣り合った重複要素だけが取り除かれる。

つまり、↓のように完全に重複が除去されない可能性がある。

a a a a a b b b b b a => a b a

cpprefjp.github.io

2. erase の引数
  u.erase(unique(u.begin(), u.end()), u.end());

上記のeraseは↓のように第二引数を忘れてもコンパイルが通ってしまう。

  u.erase(unique(u.begin(), u.end()));

この場合のeraseの挙動は、uniqueが返す index の1つの要素のみを削除する。


1 or 2 のミスを犯すことで、重複が完全に除去されずに残ってしまい、後の二分探索が正常に動作せずバグが生じる。