わりとよくある備忘録

競プロ他雑多な私的メモ

CGOを使ったGoパッケージをRaspberryPy用にCrossCompileしたい

と、思ったのですが、Macでやろうと思うとなんだか面倒だったのでクロスコンパイル用のコンテナを作った。
かなりニッチだと思いますが、せっかく作ったので整理して下記に置いた。
github.com

使い方

事前にbuildしたいソース一式をgithubに配置して置く必要があります。

使い方は簡単で、まずコンテナをビルドし、

docker build . -t gorasp:build

そしてバイナリ出力先をマウントしつつ実行。(下記ではgorasp_builder直下のbuild)

docker run -it -v "$(pwd)"/build:/go/build gorasp:build

あとは、github repositryのURLを聞かれるので入力し、CGOを使っている場合はその旨入力すればクロスコンパイルされる。

余談

この作業をするに併せて、久々にラズパイのOSを書き込んだんですが昨今はimagerなんてものがあるんですね。
マジで楽すぎて感動した。文明の進歩に感謝...。

bit DP:TSP(ABC180E)

atcoder.jp
TSP!

Traveling Salesman Problem

巡回セールスマン問題(じゅんかいセールスマンもんだい、英: traveling salesman problem、TSP)は、都市の集合と各2都市間の移動コスト(たとえば距離)が与えられたとき、全ての都市をちょうど一度ずつ巡り出発地に戻る巡回路のうちで総移動コストが最小のものを求める(セールスマンが所定の複数の都市を1回だけ巡回する場合の最短経路を求める)組合せ最適化問題である。

by 巡回セールスマン問題 - Wikipedia

全探索を行えば勿論最適解は求まるが、あり得る訪問経路の数は(N-1)!通りになるのですぐに組合せ爆発する。(例えばN-1=16なら16!≒2e+13)
wikipediaにあるようにTSPはいわゆるNP困難な問題で、多項式時間で解けるアルゴリズムは見つかっていないので現実的な時間で最適解を求められるかどうかは都市の数(=N)に依存する。

AtCoderの時間制約的にはNが10以上程度になると既に全探索では厳しくなってくるが、N≦17程度ならばDPによる高速化で十分に解ける場合が多い。

ちなみに全探索版で解ける問題はABC183Cで出てましたね。
atcoder.jp

どう解くか

というわけでDPによる高速化をしていく。

重要なポイントは、現状態から終端までの最適な訪問順は「これまでどの都市を訪問してきたか」ということと「今いる都市」という情報のみに依存し、「これまで訪問した都市の訪問順序」には依らないということ。つまり残りの未訪問の都市集合と、スタート地点が同じなら、最適な経路は当然一意に定まる。
※このエントリの図が非常に分かりやすく参考になりました。
dalekspritner.hatenablog.com


よってあり得る状態は、「これまでどの都市を訪問してきたか」という状態Sと、「今いる都市」であるvで表現でき、dp[S][v]を状態S-vに至るまでの最小コストとするとdp[S_{all}][0]を求めれば良いということになる。
遷移は、Sから今いる都市を除いた状態(現状態の前状態としてあり得る状態)をS'と表すと、dp[S][v]dp[S'][v']+cost[v'][v],(v' \in S')の中から最小のものを選べば良いので下記のようになる。

dp [S][v]=min(dp[S][v],dp[S'][v']+cost[v'][v])

bitによる状態Sの表現

上記までで遷移は分かったので、あとはS,S'を配列のindexで取り扱えるような値で表せればok。ここでSがなんだったかを振り返ると、ある都市を訪問したか否かを表現できれば良いのでbitで取り扱えば良さそう。

例えば、都市1,3,4を訪問した場合だと下記のように表せば良い。

S=\{1,3,4\} => 0b11010

これで各状態を整数値で取り扱えるようになったので、前述の遷移を実装すれば終了。

計算量

あり得る状態数は2^nであり、各状態に対して遷移n^2を集約するので計算量はO(2^n*n^2)になる。

実装例

ここまでのロジックを忠実に実装すると下記のようになるが、実際最内ループの条件分岐は無くても動く。(初期値をINFにしていることでどちらにしても更新されないので。)

#include <algorithm>
#include <cmath>
#include <cstring>
#include <deque>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <string>
#include <vector>

#define rep(i, n) for (int i = 0; i < (n); i++)
#define reps(i, n, s) for (int i = (s); i < (n); i++)

using ll = long long;
using namespace std;
constexpr long long INF = 1LL << 60;
constexpr int MAX_N = 17;

ll dp[(1 << MAX_N)][MAX_N];

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

    int n;
    cin >> n;

    vector<ll> x(n), y(n), z(n);
    rep(i, n) cin >> x[i] >> y[i] >> z[i];
    // calc cost
    vector<vector<ll>> cost(n, vector<ll>(n, 0));
    rep(i, n) {
        reps(j, n, i + 1) {
            ll c = abs(x[i] - x[j]) + abs(y[i] - y[j]);
            cost[i][j] = c + max(0LL, z[j] - z[i]);
            cost[j][i] = c + max(0LL, z[i] - z[j]);
        }
    }

    // TSP
    fill(dp[0], dp[(1 << MAX_N)], INF);
    dp[0][0] = 0;
    reps(S, 1 << n, 1) {
        rep(v, n) {
            if ((S >> v) & 1) {
                rep(j, n) {
                    ll S_hat = S - (1 << v);
                    if (S_hat >> j & 1 || S_hat == 0) { dp[S][v] = min(dp[S][v], dp[S_hat][j] + cost[j][v]); }
                }
            }
        }
    }

    cout << dp[(1 << n) - 1][0] << endl;
    return 0;
}

bit DP楽しい。

LongestCommonSubsequence(EDP-F)

atcoder.jp
有名なDPの問題なので知ってはいたが、ややこしそうなので後回しにしていた。
先日買ったけんちょんさん(@drken1215)の本を進めていたら章末問題に載っていたので、もりもり鍛えるべく重い腰を上げて解いてみた。


自分的には結構ややこしくて遷移の理解に苦しんだ。
文字列系アルゴリズムへの苦手意識が拭えない。

最長共通部分列問題(LCS Problem)

二つの異なる系列が与えられ、その共通部分列を答える問題。
ここでいう系列としては文字列などの記号列を取り扱う場合が多い。
また部分列は、順序は保持する必要はあるが隣接関係は崩れても良いというあたりがややこしい。
詳しくはwikiを参照。
最長共通部分列問題 - Wikipedia

どう解くか

dpの基本である部分問題に分割してそれを順に解いていくということをする。(当たり前..)
具体例として、EducationalDPコンテストの入力である二つの文字列S=axyb,T=abyxbを考える。

まず、S,Ti,j文字目をそれぞれs_{i},t_{j}で表すとする。
そしてdp [i][j]S_{sub_i}=s_{1}s_{2}...s{i},T_{sub_j}=t_{1}t_{2}...t_{j}最長共通部分列長さと定義する。
そうすると、i,j文字目に注目し①s_{i}==t_{j}ならば、その文字を除いた前の文字列S_{i-1},T_{j-1}に今の一致した分の長さを加えれば良いので、

dp [i][j]=dp[i-1][j-1]+1

となる。
次に②s_{i}!=t_{j}の場合を考えると、dp配列の値はS_{i-1},T_{j}のLCSとS_{i},T_{j-1}のLCSのうちより長いほうの値がそのまま現状態の長さになり、遷移は下記のようになる。
dp [i][j]=max(dp[i-1][j],dp[i][j-1])

この遷移を適用し、dp配列を埋めていくと下記表のようになる。
ここで矢印は更新に適用された遷移であり、""は空文字を表している。
f:id:KLEE:20201010180407j:plain
この表の上では斜めの遷移が①(注目文字の追加)の遷移、縦横の遷移が②(注目文字の削除)となっている。
これはS_{1}=a,T_{4}=axybを一致させるにはt_{2}t_{3}t_{4}を削除する必要があることを考えると分かりやすい。

実装例

遷移がここまで書ければあとは実装するだけだが、上記までの処理ではLCSの長さしか分からない。
復元が求められる場合は、表の右下から見ていって丸が着いた部分に遭遇したら斜め左上に遷移し、そうで無い部分のときは縦横のうち自信と同じ値を持つ方向に遷移するというポリシーで表の左上まで辿っていけばLCSそのものが得られる。
この時、縦横のうち両方とも同じ値の場合にどちらを優先するかで復元されるLCSは異なる物になる場合がある(LCSは一般に複数パターンあり得る)。
復元も含めて実装した例は下記のようになる。

#include <iostream>
#include <algorithm>
#include <string>
#include <vector>

#define rep(i, n) for (int i = 0; i < (n); i++)
using namespace std;

int main() {
  string s,t;
  cin>>s>>t;
  int n=s.size();
  int m=t.size();
  vector<vector<int>>dp(n+1,vector<int>(m+1,0));

  //各状態でのLCS長さをメモ
  rep(i,n){
      rep(j,m){
        if(s[i]==t[j]){//同じ文字なら
            dp[i+1][j+1]=dp[i][j]+1;
        }
        else{
            dp[i+1][j+1]=max(dp[i][j+1],dp[i+1][j]);
        }
      }
  }
  //後ろから辿っていきLCSを復元
  int i=n,j=m;
  string ans="";
  while(i>0&&j>0){
      if(s[i-1]==t[j-1]){
          ans+=s[i-1];
          i--;j--;
      }
      else{
          if(dp[i][j]==dp[i-1][j]){//横方向優先
              i--;
          }
          else{
              j--;
          }
      }
  }
  //後ろ向きにメモしたので逆順に直す
  reverse(ans.begin(),ans.end());

  cout<<ans<<endl;
  return 0;
}

ややこしい...

DP~imos法による高速化(ABC179D)

atcoder.jp

素朴に配るDPを実装するとO(N^2)になりTLEとなる(なってしまった....)。
例えば下記の入力だと遷移数がマス目の長さと同じNオーダになるので通らない。

200000 1
1  200000

こういう時は累積和でdpを高速化するのが定番らしい。

解説動画が非常に分かりやすいが、貰うdpで解説されていたので自分が本番の時に書こうとしていた配るdpで実装する場合どうすれば良いのか調べてACしたのでメモする。
■解説動画
www.youtube.com

結論から書くとimos法で累積和をシミュレートしながら更新していけば良い。

imos法による配るdpの高速化

imos法自体については考案者の方のページが非常に分かりやすい。
imoz.jp

自分の理解だと、状態が遷移する境界値だけメモしておきシミュレートすることで区間系累積和の高速化などに応用できる手法という理解。
区間数をK,マス目の長さをNとすると、境界値のメモは境界を生成する区間数に比例するO(K)となり、シミュレートは頭から足してくだけなのでO(N)となる。
よって全体としてはO(N+K)で累積和が計算できる。

今回のD問題ではマス目の頭から見ていき、各位置i毎に自分の位置から右側へ区間の境界値を配りシミュレートしながら更新していけば良い。
よってマス目の長さN分だけ遷移区間K個分の境界値のメモを行うので、全体のオーダはO(NK)となり間に合う。

実装例

全体のコードを載せるとロジックが見にくいのでmain部分のみ。全コードはここ

int main()
{
  int n, k;
  cin >> n >> k;
  vector<pair<int, int>> v(k);
  rep(i, k)
  {
    cin >> v[i].first >> v[i].second;
  }
  vector<mint> dp(n + 1, 0);
  dp[0] = 1;
  reps(i, n, 1)
  {
    //現在マスの値を更新
    dp[i] += dp[i - 1];

    //現在マスの値を配る
    for (auto p : v)
    {
      //境界インデックスを計算
      int l = i + p.first;
      int r = i + p.second + 1;
      if (l > n) continue;//左境界の領域外判定
      dp[l] += dp[i];
      if(r > n+1) continue;//右境界の領域外判定
      dp[r] -= dp[i];
    }
  }
  cout << dp[n] << endl;
  return 0;
}

累積和用の配列とdp用の配列を共通にして実装したが、バグを仕込みかねないので分けた方がいいのかもしれない。

包除原理(ABC172E)

包除原理を応用する問題が出題されて本番中は手も足も出なかった。
いろんな解説を見てACしたので理解をまとめる。
atcoder.jp

包除原理とは

ものすごくざっくり書くと、複数の集合から構成される全体集合を重複なく求めるための考え方。基本的には全体を足して、重複する部分に関しては別途足し引きして辻褄を合わせる。
詳細は下記を参照。
mathtrain.jp

どう応用したか

今回の例では、ユニークな数字で構成される2つの数列A,Bの内、任意のiA_i≠B_iとなるパターンを数え上げる。
このような数列を直接求めるのは難しいので、全体の数列の組み合わせを全通り求めてそこから条件を満たさないものを除いていくことを考える。

全体集合

全体の数列を表す集合Uは、m個の中からn個を2数列分並べればよいので次の式で求まる。

U=( _m P_n)^2

条件を満たさない集合

今回の条件を満たさない集合とは、任意のiA_i=B_iとなるiが1つ以上ある集合のことなのでこれを数え上げていく。
ここで、直接「ちょうどk箇所のiで重複する集合」を求められれば良いが、これも難しい。「少なくともk箇所のiで重複する集合」なら数えやすそうだが、これを全て足してしまうと重複する部分が発生する。例えば「少なくとも1箇所のiで重複する集合」には「少なくとも2箇所のiで重複する集合」も含まれている。

ここで、包除原理を使って重複を削除し「k箇所のiで値が重複する集合S,k=\{1,2,...,n\}」を数え上げればよい。
上記を計算するために、「少なくともk箇所のiで重複する集合S_k」について考える。
これを数え上げるのは簡単で、次のように考える。

  1. まず、長さnの数列の中で重複させる位置をk箇所選ぶ ->  _n C_k
  2. 次に、そのk箇所に重複させる数字をm個の中から選んで並べる ->  _m P_k
  3. 最後に、残ったn-kの位置に当てはめる数字をm-k個の中から選んで並べる(数列A,B2つ分を選ぶため二乗となることに注意) -> ( _{m-k} P_{n-k})^2

よって重複がkのときの集合S_kは下記のようになる

 S_k=_n C_k* _m P_k*( _{m-k} P_{n-k})^2

ここまでくれば、Sを求めるにはS_kを包除原理に従って足し合わせればよい。

S=\sum_{k=1}^{n}(-1)^{(k-1)}* _n C_k* _m P_k*( _{m-k} P_{n-k})^2

よって今求めたい集合は全体からこれを引けばよいので、

U-S=( _m P_n)^2 - \sum_{k=1}^{n}(-1)^{(k-1)}* _n C_k* _m P_k*( _{m-k} P_{n-k})^2

となる。ここで、k=0のときは、
_n C_k* _m P_k*( _{m-k} P_{n-k})^2=( _m P_n)^2

となるので、上式を整理すると結局、
\sum_{k=0}^{n}(-1)^k* _n C_k* _m P_k*( _{m-k} P_{n-k})^2

を求めればよいことになる。

実装例

#include <bits/stdc++.h>

#define rep(i, n) for (int i = 0; i < (n); i++)
#define reps(i, n, s) for (int i = (s); i < (n); i++)

using ll = long long;
using namespace std;
constexpr long long MAX = 5100000;
constexpr long long INF = 1LL << 60;
constexpr int MOD = 1000000007;

ll fac[MAX], finv[MAX], inv[MAX];

void COMinit() {
  fac[0] = fac[1] = 1;
  finv[0] = finv[1] = 1;
  inv[1] = 1;
  for (int i = 2; i < MAX; i++) {
    fac[i] = fac[i - 1] * i % MOD;
    inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
    finv[i] = finv[i - 1] * inv[i] % MOD;
  }
}

ll nCk(int n, int k) {
  if (n < k) return 0;
  if (n < 0 || k < 0) return 0;
  return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD;
}

ll nPk(ll n, ll k) {
  if (n == 0 || k == 0) return 1;
  return fac[n] * finv[n - k] % MOD;
}

class mint {
  long long x;

 public:
  mint(long long x = 0) : x((x % MOD + MOD) % MOD) {}
  mint operator-() const { return mint(-x); }
  mint& operator+=(const mint& a) {
    if ((x += a.x) >= MOD) x -= MOD;
    return *this;
  }
  mint& operator-=(const mint& a) {
    if ((x += MOD - a.x) >= MOD) x -= MOD;
    return *this;
  }
  mint& operator*=(const mint& a) {
    (x *= a.x) %= MOD;
    return *this;
  }
  mint operator+(const mint& a) const {
    mint res(*this);
    return res += a;
  }
  mint operator-(const mint& a) const {
    mint res(*this);
    return res -= a;
  }
  mint operator*(const mint& a) const {
    mint res(*this);
    return res *= a;
  }
  mint pow(ll t) const {
    if (!t) return 1;
    mint a = pow(t >> 1);
    a *= a;
    if (t & 1) a *= *this;
    return a;
  }
  // for prime MOD
  mint inv() const { return pow(MOD - 2); }
  mint& operator/=(const mint& a) { return (*this) *= a.inv(); }
  mint operator/(const mint& a) const {
    mint res(*this);
    return res /= a;
  }

  friend ostream& operator<<(ostream& os, const mint& m) {
    os << m.x;
    return os;
  }
};
int main() {
  cin.tie(0);
  ios::sync_with_stdio(false);
  ll n, m;
  cin >> n >> m;

  COMinit();
  // nCk*mPk*(m-kPn-k)^2
  mint ans = 0;
  rep(k, n + 1) {
    mint now = nPk(m - k, n - k);
    now *= now;
    now *= nPk(m, k);
    now *= nCk(n, k);
    if (k % 2 != 0) now = -now;
    ans += now;
  }
  cout << ans << endl;
  return 0;
}

参考

以下、参考にさせて頂いたサイトです。
www.youtube.com
qiita.com

回文検出(Manacher)

何を解きたいか

ご存じのように回文とはaba,cladalcのような逆順でも変化しない文字列のことです。このとき、ある文字列の部分文字列で作られる回文の最大長さや数を問う問題を解きたい場合を考えます。ナイーブに解こうとすると二重のfor文で全ての部分文字列S_{sub}を生成し、S_{sub}=reverse(S_{sub})をチェックすればよいですが、O(|S|^2)となるので短い文字列にしか適用できません。

この問題を回文の特性を上手に利用してO(|S|)で解く方法があり、それが今回説明するManacherのアルゴリズムです。

どう解くか

上述したように回文の特性を上手いこと利用します。 説明を簡単にするために奇数長の回文のみを検出対象とします。偶数への拡張は別途最後にコード例を示します。

今回は、具体例としてabacababaという文字列で考えます。

この文字列を頭から見ていって、各文字を中心に両端に何文字同じ文字を持つかをメモしていきます。これは、メモされた値をL_{i}とすると、その文字を中心に2L_{i}-1の長さの回文があるという情報をメモしたことになります。

具体的に図で見ていくと次のような形で処理が進んでいきます。

◆STEP1
1文字目はaで、右側にしか文字が存在せず回文の長さは自分自身だけのため、1をメモします。
f:id:KLEE:20200617231150p:plain
◆STEP2
2文字目はbで、両側にaが存在するため、回文の長さは自分自身+1として2をメモします。
f:id:KLEE:20200617231154p:plain
◆STEP3
3文字目はaで、両側には同じ文字は存在しないため、1をメモします。
f:id:KLEE:20200617231156p:plain
◆STEP4
4文字目はcで、両側は矢印のaの部分まで同じ文字が続くので、4をメモします。
f:id:KLEE:20200617231138p:plain

 この先STEP5,6と進めていきたいところですが、文字列をよく見るとSTEP4でcを中心にabacabaという回文が構成されていることが分かっているため、STEP5,6の結果はSTEP2,3と同じになり探索不要なことが分かります。

◆STEP5,6
f:id:KLEE:20200617231141p:plain
f:id:KLEE:20200617231144p:plain


このように、ある位置iで長さL_{i}の回文が検出されたときi-1までの計算結果を用いることで探索を効率化しようというのが基本アイデアです。
ただここで、利用できる過去の計算結果の範囲には気を付ける必要があります。
例えばSTEP7ではわかりやすくそういったことが起きます。

◆STEP7
このSTEPのメモ値は3となり、着目文字の対称にある左端のaのメモ値とは一致しません。
f:id:KLEE:20200618204048p:plain

ある文字を中心に一文字ずつ両側に同じ文字があるかどうかを探索していくとき、現在着目してる文字を含む最長の回文(STEP7ではabacaba)より内側までしか対称性が保証されません。そのため、最長の回文より外の範囲へ探索範囲が伸びた時点でそこから先は未知の世界のため、真面目に探索する必要があります。

実装例

これらを実装したコード例を下記に示します。

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

int main() {
  string s = "abacababa";
  int n = s.size();

  vector<int> rad(n);
  int c = 0, r = 0;
  
  while (c < n) {
    // cを中心に同じ文字がどこまで連続するか
    while (0 <= c - r && c + r < n && s[c - r] == s[c + r]) r++;
    rad[c] = r;

    //回文の長さに応じて利用可能な範囲を確認しつつメモ
    int k = 1;
    while (0 <= c - k && k + rad[c - k] < r) {
      rad[c + k] = rad[c - k];
      k++;
    }
    //すでに計算が終わった分だけ中心と探索半径をずらす
    c += k;
    r -= k;
  }
  
  int Lmax = *max_element(rad.begin(), rad.end()) * 2 - 1;
  cout << Lmax << endl;
  return 0;
}

偶数長への拡張

偶数長回文への拡張は簡単で、元の文字列の各文字の間に#を挟んで同じアルゴリズムで検出すればよいです。そうすると#のインデックスには偶数長の回文の長さが、元の文字のインデックスには奇数長の回文の情報がメモされます。
ただし、最後に#分の長さを差し引きしなければならないことに注意が必要です。

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

int main() {
  string _s = "abacabaaa";
  int n = _s.size();
  //'#'の追加
  string s(2 * n + 1, '#');
  for (int i = 0; i < n; i++) s[2 * i + 1] = _s[i];
  n = 2 * n + 1;

  vector<int> rad(n);
  int c = 0, r = 0;

  while (c < n) {
    // cを中心に同じ文字がどこまで連続するか
    while (0 <= c - r && c + r < n && s[c - r] == s[c + r]) r++;
    rad[c] = r;

    //回文の長さに応じて利用可能な範囲を確認しつつメモ
    int k = 1;
    while (0 <= c - k && k + rad[c - k] < r) {
      rad[c + k] = rad[c - k];
      k++;
    }
    //すでに計算が終わった分だけ中心と探索半径をずらす
    c += k;
    r -= k;
  }
  //'#'分の補正
  for (int i = 0; i < n; i++) {
    if (i % 2 == 0)
      rad[i] = (rad[i] - 1) / 2;
    else
      rad[i] /= 2;
  }

  int Lmax = *max_element(rad.begin(), rad.end()) * 2 - 1;
  cout << Lmax << endl;
  return 0;
}

参考

以下、参考にさせて頂いたサイトです。
snuke.hatenablog.com

よくある理由で、

はじめました。
競プロやったり、自作キーボード組んだりするのが好きです。

一応今年中にAtCoderで水色到達を当面の目標にしてるんですが、アルゴリズムの勉強とかしたそばから忘れていくので、記憶定着目的でいろいろ書いてく。