昨日のABC188、Dが解けず３完だった。
いもす法と座表圧縮かなと思いつつ、座標圧縮への理解が曖昧だった＋本番の焦りもありACできなかった。
最近なんとなく解法が浮かぶのは成長だなと思いつつも、あと一歩詰めきれないところが甘いので書いていく。
座標圧縮とはある数列に対して位置関係や大小関係などを保存しつつ、より小さい数値で表現し直す手法。
例えば、という数列に対して大小関係を保持しつつ数値変換を行うと
のように表現できる。
この時、当然スケールの情報は失われるが位置関係や大小関係が重要な問題などでは有効な場合がある。
今回の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; }
と、思ったのですが、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なんてものがあるんですね。
マジで楽すぎて感動した。文明の進歩に感謝...。
atcoder.jp
TSP！
巡回セールスマン問題（じゅんかいセールスマンもんだい、英: traveling salesman problem、TSP）は、都市の集合と各2都市間の移動コスト（たとえば距離）が与えられたとき、全ての都市をちょうど一度ずつ巡り出発地に戻る巡回路のうちで総移動コストが最小のものを求める（セールスマンが所定の複数の都市を1回だけ巡回する場合の最短経路を求める）組合せ最適化問題である。
全探索を行えば勿論最適解は求まるが、あり得る訪問経路の数は通りになるのですぐに組合せ爆発する。(例えば
なら
)
wikipediaにあるようにTSPはいわゆるNP困難な問題で、多項式時間で解けるアルゴリズムは見つかっていないので現実的な時間で最適解を求められるかどうかは都市の数()に依存する。
AtCoderの時間制約的にはが10以上程度になると既に全探索では厳しくなってくるが、
程度ならばDPによる高速化で十分に解ける場合が多い。
ちなみに全探索版で解ける問題はABC183Cで出てましたね。
atcoder.jp
というわけでDPによる高速化をしていく。
重要なポイントは、現状態から終端までの最適な訪問順は「これまでどの都市を訪問してきたか」ということと「今いる都市」という情報のみに依存し、「これまで訪問した都市の訪問順序」には依らないということ。つまり残りの未訪問の都市集合と、スタート地点が同じなら、最適な経路は当然一意に定まる。
※このエントリの図が非常に分かりやすく参考になりました。
dalekspritner.hatenablog.com
よってあり得る状態は、「これまでどの都市を訪問してきたか」という状態と、「今いる都市」である
で表現でき、
を状態
-
に至るまでの最小コストとすると
を求めれば良いということになる。
遷移は、から今いる都市を除いた状態(現状態の前状態としてあり得る状態)を
と表すと、
は
,(
)の中から最小のものを選べば良いので下記のようになる。
上記までで遷移は分かったので、あとは,
を配列のindexで取り扱えるような値で表せればok。ここで
がなんだったかを振り返ると、ある都市を訪問したか否かを表現できれば良いのでbitで取り扱えば良さそう。
例えば、都市1,3,4を訪問した場合だと下記のように表せば良い。
あり得る状態数はであり、各状態に対して遷移
を集約するので計算量は
になる。
ここまでのロジックを忠実に実装すると下記のようになるが、実際最内ループの条件分岐は無くても動く。（初期値を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楽しい。
atcoder.jp
有名なDPの問題なので知ってはいたが、ややこしそうなので後回しにしていた。
先日買ったけんちょんさん(@drken1215)の本を進めていたら章末問題に載っていたので、もりもり鍛えるべく重い腰を上げて解いてみた。
アルゴリズムとデータ構造をげっつしたのでもりもり鍛える pic.twitter.com/lmv1eW06JN
— KLEE (@kashee337) October 1, 2020
二つの異なる系列が与えられ、その共通部分列を答える問題。
ここでいう系列としては文字列などの記号列を取り扱う場合が多い。
また部分列は、順序は保持する必要はあるが隣接関係は崩れても良いというあたりがややこしい。
詳しくはwikiを参照。
最長共通部分列問題 - Wikipedia
dpの基本である部分問題に分割してそれを順に解いていくということをする。(当たり前..)
具体例として、EducationalDPコンテストの入力である二つの文字列,
を考える。
まず、,
の
,
文字目をそれぞれ
,
で表すとする。
そしてを
,
の最長共通部分列長さと定義する。
そうすると、,
文字目に注目し①
ならば、その文字を除いた前の文字列
,
に今の一致した分の長さを加えれば良いので、
遷移がここまで書ければあとは実装するだけだが、上記までの処理では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を実装するとになりTLEとなる(なってしまった....)。
例えば下記の入力だと遷移数がマス目の長さと同じオーダになるので通らない。
200000 1 1 200000
こういう時は累積和でdpを高速化するのが定番らしい。
解説動画が非常に分かりやすいが、貰うdpで解説されていたので自分が本番の時に書こうとしていた配るdpで実装する場合どうすれば良いのか調べてACしたのでメモする。
■解説動画
www.youtube.com
結論から書くとimos法で累積和をシミュレートしながら更新していけば良い。
imos法自体については考案者の方のページが非常に分かりやすい。
imoz.jp
自分の理解だと、状態が遷移する境界値だけメモしておきシミュレートすることで区間系累積和の高速化などに応用できる手法という理解。
区間数を,マス目の長さを
とすると、境界値のメモは境界を生成する区間数に比例する
となり、シミュレートは頭から足してくだけなので
となる。
よって全体としてはで累積和が計算できる。
今回のD問題ではマス目の頭から見ていき、各位置毎に自分の位置から右側へ区間の境界値を配りシミュレートしながら更新していけば良い。
よってマス目の長さ分だけ遷移区間K個分の境界値のメモを行うので、全体のオーダは
となり間に合う。
全体のコードを載せるとロジックが見にくいので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用の配列を共通にして実装したが、バグを仕込みかねないので分けた方がいいのかもしれない。
包除原理を応用する問題が出題されて本番中は手も足も出なかった。
いろんな解説を見てACしたので理解をまとめる。
atcoder.jp
ものすごくざっくり書くと、複数の集合から構成される全体集合を重複なく求めるための考え方。基本的には全体を足して、重複する部分に関しては別途足し引きして辻褄を合わせる。
詳細は下記を参照。
mathtrain.jp
今回の例では、ユニークな数字で構成される2つの数列の内、任意の
で
となるパターンを数え上げる。
このような数列を直接求めるのは難しいので、全体の数列の組み合わせを全通り求めてそこから条件を満たさないものを除いていくことを考える。
全体の数列を表す集合は、m個の中からn個を２数列分並べればよいので次の式で求まる。
今回の条件を満たさない集合とは、任意ので
となる
が１つ以上ある集合のことなのでこれを数え上げていく。
ここで、直接「ちょうど箇所の
で重複する集合」を求められれば良いが、これも難しい。「少なくとも
箇所の
で重複する集合」なら数えやすそうだが、これを全て足してしまうと重複する部分が発生する。例えば「少なくとも
箇所の
で重複する集合」には「少なくとも
箇所の
で重複する集合」も含まれている。
ここで、包除原理を使って重複を削除し「箇所の
で値が重複する集合
,
」を数え上げればよい。
上記を計算するために、「少なくとも箇所の
で重複する集合
」について考える。
これを数え上げるのは簡単で、次のように考える。
よって重複がのときの集合
は下記のようになる
ここまでくれば、を求めるには
を包除原理に従って足し合わせればよい。
よって今求めたい集合は全体からこれを引けばよいので、
#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