わりとよくある備忘録

競プロ他雑多な私的メモ

AtCoder用のRustテンプレートを作った

Rustの勉強を最近初めた。
準備運動代わりにAtCoderの問題を解いたりすることがあるのだが、自分にとって使い勝手の良いtemplateが見当たらなかった。
というわけで作ってみた。

github.com

  • コンテスト用に1プロジェクトで5ファイルのbinを生成(a〜f.rs)
  • vscode用に整理しており各ファイル毎にF5でbuild->debug可能

使い方

$ cargo generate --git https://github.com/kashee337/rust_atcoder_template.git

参考

Cargo.tomlのdependencyAtCoderで使えるクレートについては下記のrepositryを参考にさせていただいた。
github.com

ぼちぼちやってく

座表圧縮 (ABC188D)

昨日のABC188、Dが解けず3完だった。

atcoder.jp

いもす法と座表圧縮かなと思いつつ、座標圧縮への理解が曖昧だった+本番の焦りもありACできなかった。
最近なんとなく解法が浮かぶのは成長だなと思いつつも、あと一歩詰めきれないところが甘いので書いていく。

座標圧縮

座標圧縮とはある数列に対して位置関係や大小関係などを保存しつつ、より小さい数値で表現し直す手法。
例えば、 30,100000,500000,2という数列に対して大小関係を保持しつつ数値変換を行うと2,3,4,1のように表現できる。
この時、当然スケールの情報は失われるが位置関係や大小関係が重要な問題などでは有効な場合がある。

どう適用するか

今回のABC188 Dは区間の問題である。問題文から、imos法を知っていれば明らかにそれ(いもす法で解くやつ)だなとわかる。
ただし、今回はありうる区間の座標が1e9のため、そのままではシミュレートができない。
そこで座標圧縮で区間の位置関係を保持しつつ小さい値へ変換を行う。
具体的には、あり得る座標群を取得し、それに対して昇順でindexを付ければ良い。
C++ならば例えば下記のように実現できる。

/* a,bは区間の始点終点を保持したvector */
//あり得る座標群を重複なしで取得
set<int> day_set;
rep(i, n) {
    day_set.insert(a[i]);
    day_set.insert(b[i]+1); // いもす法用に+1しておく
}
//座標群を昇順にsort
vector<int> day_list(day_set.begin(), day_set.end());
sort(day_list.begin(), day_list.end());
//この時点でday_listに保持された順序が圧縮後の座標になる

これで座標圧縮できたが、いもす法を行うために実座標->圧縮座標への変換用のhash mapを持っておくと便利なので作っておく。

//実座標->圧縮座標への変換hash mapを作る
map<int, int> real2comp;
rep(i, L) { real2comp[day_list[i]] = i; }

これで前準備ができたので、いもす法を行なっていく

vector<long long> dp(L + 1, 0);
//メモ
rep(i, n) {
    dp[real2comp[a[i]]] += c[i];
    dp[real2comp[b[i]+1]] -= c[i];
}
//シミュレート
rep(i, L) { dp[i + 1] += dp[i]; }

最終的に計算すべき料金は期間を通しての合計料金なので、実座標での期間を考慮しつつ加算していけばok

// 元座標の区間長さを考慮しながら加算していく
long long res = 0;
rep(i, L - 1) { res += (day_list[i + 1] - day_list[i]) * min(C, dp[i]); } //prime料金Cとdp[i]の低いほうを選ぶのが最適

全コード

#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++)
#define rrep(i, n) for (int i = (n - 1); i >= 0; i--)
#define rreps(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;

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    ll n, C;
    cin >> n >> C;
    vector<ll> a(n), b(n), c(n);
    rep(i, n) {
        cin >> a[i] >> b[i] >> c[i];
        b[i] += 1;  //いもす法での立ち下がり点
    }

    //あり得る座標群を重複なしで取得
    set<ll> day_set;
    rep(i, n) {
        day_set.insert(a[i]);
        day_set.insert(b[i]);
    }
    //座標群を昇順にsort
    vector<ll> day_list(day_set.begin(), day_set.end());
    sort(day_list.begin(), day_list.end());
    ll L = day_list.size();  //座標圧縮後の座標の長さ

    //実座標->圧縮座標への変換hash mapを作る
    map<ll, ll> real2comp;
    rep(i, L) { real2comp[day_list[i]] = i; }

    //圧縮後の座標でいもす法をやっていく
    vector<ll> dp(L + 1, 0);
    //メモ
    rep(i, n) {
        dp[real2comp[a[i]]] += c[i];
        dp[real2comp[b[i]]] -= c[i];
    }
    //シミュレート
    rep(i, L) { dp[i + 1] += dp[i]; }

    // 元座標の区間長さを考慮しながら加算していく
    ll res = 0;
    rep(i, L - 1) { res += (day_list[i + 1] - day_list[i]) * min(C, dp[i]); }

    cout << res << endl;

    return 0;
}

おまけ

今回の問題に関しては解説における解のように、sort+変化点だけに着目して前から見てく方式の方がコード自体はシンプル。

#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 pb push_back
#define rep(i, n) for (int i = 0; i < (n); i++)
#define reps(i, n, s) for (int i = (s); i < (n); i++)
#define rrep(i, n) for (int i = (n - 1); i >= 0; i--)
#define rreps(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;

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    ll n, C;
    cin >> n >> C;

    vector<pair<ll, ll>> event;

    rep(i, n) {
        ll a, b, c;
        cin >> a >> b >> c;
        event.pb({a, c});
        event.pb({b + 1, -c});
    }
    sort(event.begin(), event.end());

    ll res = 0, t = 0, fee = 0;
    for (auto v : event) {
        if (v.first != t) {
            res += min(C, fee) * (v.first - t);
            t = v.first;
        }
        fee += v.second;
    }

    cout << res << endl;
    return 0;
}

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