最小費用流についてのちょっとしたメモ
恥ずかしながら、最小費用流をかなりブラックボックスとして利用していたのだが、今まであまり気付いてなかったことたちが少し整理できたので雑なメモを残しておきます。アルゴリズム自体は各種資料を見てください(怠惰) 双対 LP との関係も知る人ぞ知るジャンルではあるのですが詳しい資料がたくさんあるのでそちらへ
最小費用流のライブラリを整備したい人は、これが今まで見た中で最も一般的な形をしている(fが負って何?)ので頑張って通すと良さそう。 Minimum cost b-flow
フローについて詳しいことが知りたい人には Mi_Sawa さんの記事たちがかなり参考になり、今回はこれを参考にしました
ぼくの考えたさいきょうのフローライブラリ - みさわめも
みなさん大好き組合せ最適化も参考にしました。hosさんやwataさんやtokoharuさんの資料も有名です。
b-flow というのは、各頂点に需給が設定されており、各頂点周りで辺の流量のin/outのバランスがbとなっているものである。(組合せ最適化においてもこの定式化)
これをよく見るs-tフローの形にするには、b(v)>0なるvたちにはs->vに容量b(v)の辺を張り、b(v)<0なるvたちにはv->tに容量-b(v)の辺を張って、s->tにbが正のものの和だけ流せばよい。
最短路問題のようなものを考える場合でも、負の閉路があるかどうかというのが重要で、ないものは保守的(組合せ最適化より)と呼ばれ、そのような場合実行可能ポテンシャルが取れる。 (負の閉路が無いことと同値である) 保守的である場合、例えば Bellman-Ford などで適当にやると実行可能ポテンシャルの1つが取れる。
最小費用流においても残余グラフにおいて負の費用の閉路がないことが最適性を特徴付けるのだが、SSP(蟻本もこれだっけ)のアルゴリズムでは各ステップであるb'があって最適なb'-flow、となる形で流量を保持しつつ、この実行可能ポテンシャルも同時に保持しているということがある(これもあまりわかってなかった) このアルゴリズムでは、cだと負辺がある場合もあるが、前のステップのポテンシャルを用いて簡約距離 (x->yに重みwの辺があるときh(x)+w-h(y) みたいな量、ポテンシャルになってればこれは非負)についての最短距離を求めても同じであることを利用して Dijkstra を使える形にしている (天才)
ちなみに、最小費用流の双対LPを見ると実行可能ポテンシャルですね〜という形をしていて、SSPで得られるポテンシャルが双対問題の解にもなっています。
負の辺があると...怖い!というのは有名なはずだがあまり対処方法は有名でないような。結論から言えば、負のコスト(かつ容量正)の辺を逆に全部使ってしまったことにしてbの値を調整し、逆辺があったことにすると残余グラフが保守的であることにできます。あとは1回Bellman-Fordした後はDijkstraでできます。SSPの計算量は|b(v)|の値に依存してるので計算量に注意が必要ではあります。
これは、上で述べたように、1回ポテンシャルの形にできてしまえば後は非負のままできるというわけで、何もしなくても保守的なグラフであれば、すべての容量が0でもbが0ベクトルのときのb-flowの形になってるので初めのポテンシャルを求めてあげればいいわけです。
容量下限がある場合にも、下駄をはかせたりして解決しがちですが問答無用で下限分だけ使ってbを調整した後上限-下限の辺があることにしてオラクルを呼べばいいです。これも計算量に注意 (misawaさんの記事には逆辺の上限を弄る方法が紹介されていて、天才)
なお、yosupo judgeのはbがでかいので scaling 系とかをやらないと通りません。やりおるわ
拡張ユークリッドの互除法の解の範囲
Modint を実装する過程で、拡張ユークリッドで得られる解の範囲について気になったことがあったことを思い出したので (かつ、あまり言及されているのを見たことがなかったので) メモを残しておきます。
- の解
とします (そうでないときは gcd で割ったり符号を変えて解いた後戻せばいいです(gcd の方は割らずに進めてもいいです))
結論から言うと、 を満たす解が得られ、過程では overflow を気にする必要はないです。
などすると他の解は得られるので、ギリギリの解が求まっていますね。
- 雑な説明
とすると、 より
が成り立っているとき、 とすると解が得られます。
だとすると、 は言えます。
の方ですが、 を使うと言えます。 は符号が逆な感じになってます
ベースケースや境界をちゃんと言うのは許して
同様にして gcd が 1 でないときも が言えると思います。
- 宣伝
は今年の自信作なのでまだの人は解いて
ICPC 国内予選 2019 参加記
参加記を書き書き。 (激寒)
5年連続チーム catsatmat で参加して国内9位で予選通過しました。
- 練習
数回しか3人では練習してなかった (国内予選/アジア地区はかなり埋まっていたのでCF gymなど5時間セット中心) あとはopencupにたまに出た (opencupはcamypaper さんと出ることが多かったが)
- 当日
実習で遅刻しかけてヒヤヒヤしたがまあ自分がいなくても行けるのでは?とも思っていた。間に合ってホッとした。
開始したらcatupperが印刷にいき、shiatsumatに設定してもらいつつAを読み、自分が通した。Bを帰ってきたcatupperに投げて自分はCを読む。shiatsumatには後ろを読んでもらう。
Dで困っていたので概要を聞き、どうせDPでしょと言いながらCを詰めているとcatupperが合流して漸化式を立ててくれていた。自分がCの見積もりになぜか失敗して無限に実行が進まず途方にくれていたら気づいたのでDに代わる。途中で数回改善を思いついたのでパソコンを奪いつつやっていきCとDが同じような時間に終わった。C, D同時並行してる間にshiatsumatにE以降を読んでもらった。
Eは頑張る系だということでshiatsumatに実装を考えてもらいつつF以降を検討していく。F全然わからなくて焦っていたらcatupperがGの解法を思いつくが実装が大変と主張していた。
ここまで相当出遅れていたので順位表を見た感じでEをやり切る必要があるのではという結論に至り、しかし自分はEを丸投げしてF絶対解くマンをしていた。catupperと議論もしたが全然まともな解法が出ず、しかし気付いたら解法を思いついており実装も軽そうなのでEがバグったところでパソコンを奪い通す。
あとはE通さないとやばいかな〜といいつつGを完全に諦める。Eだいたいできていそうだったので手でケースを提案したり、if文コメントアウト二分探索 (サイアク) をしたりする。あと数分というところでEで判定パートではなく列挙パートがミスっていることに気付く (こんなの気付かないでしょ)。コンテスト終了。かろうじて10位に入っており安堵。
- 反省
今まで本番として3人で出た大会の中では1番相談が上手く行ったのではないか。自分がCをバグらせていたのが相当罪が重いがFが解けて±0くらいにはなった。Eも複数人でチェックしたらいい感じであと5分あれば通ったんじゃないかなあ。ABCDFをさっと通してEとGを分担して詰めるみたいなのが出来ればよかったが... コンテストの成績はかなり危なかったがまあこういうこともあるでしょくらいの範囲には収まっていそう(実力不足) まだアジアもあるので5年間の総括は今度にします。
全国統一プログラミング王決定戦本戦 G - Greatest Journey
コンテストではFで辺のコストをmin(C[i],C[j])だと勘違いして破滅しました(こちらも解けるらしいのですが) サンプルに書いてあるやん(絶望)と最後の方に気付いて悲しいね。Gの解説がア(上級者向け)なのでメモを残しておきます。
まず、どういう移動が最適かを考えます。
どの辺も複数回通ってもいいという制約のもと、通った辺のコストが順にa_1,a_2,...,a_Lだとします。
この中のmax=Mをとるindex(複数ある場合は1番左)をpとすると、a_{p+1}...,a_LをMに置き換えてよい(途中から往復し続ける)ことがわかります。
Mに最短で至るまでの経路は全て1回は使わないといけないが、M以前で複数回使うような辺があったとき重複分をMにやってもらった方がいいこともわかるので、iさんはV_iからある頂点まで最短で移動した後、同じ辺(しかもそこまでの経路でmax)を使い続ける形のみを考えればよいです。
次に、こんなの分割統治(重心分解)するしか、ないだろ。という気持ちになるところは多分慣れです。
重心分解では、重心を通る移動のみを考えて後は各部分木に分ければいいので、一旦重心に集めます。V_iから重心Gまでの辺の数をE_i, コストの和をC_iと置きます。
重心GからDFSして各頂点vまでの辺の数をD_v, コストの和をB_v, vの親をp_v, vの親とのコストをw_vとし、V_iからGを経由してvまで行き、vとp_vの間を0回以上往復します。(vまで行って親の方向に往復すると考える方が考えやすいです(子はたくさんあることもあるので))
このとき各クエリiについて、C_i+(L_i-E_i-D_v)*w_v+B_v = C_i+(L_i-E_i)*w_v+(B_v-D_v*w_v) のvについての最大値(vまで辿りつける、つまりL_i - E_i - D_v >= 0の条件の下で)を速く求めたいです。見やすさのためにL_i-E_i をL'_i, B_v-D_v*w_vをB'vとおき、iについての定数項(C_i)を無視すると max L'_i * w_v + B'v (ただし L_i - E_i - D_v >= 0) となり、直線y=w_v*x+B'vたちのうちx=L'_iで1番上にあるものがわかればよいです。
まだ考慮できる頂点にL_i - E_i - D_v >= 0の条件があるため好きな順番に直線を追加する必要があり大変に思えますが、上で見たようにGからvまでのパスにおいてmaxではない辺を複数回使うものは考えなくてよく、maxの辺のみを使うことを再び考えます。(考慮する頂点のうち)適当にvを取ります。x=d(D < D_v)でvが選ばれないことが言えると嬉しいです。素直にvまでの途中の深さdの点をu (D_u = d)とします。uまで潜ったコストはB_u、vでdを代入した値は(D_u - D_v) * w_v + B_vですが、w_vがパス上のmaxという仮定より(B_v-B_u) <= w_v*(D_v - D_u)であるのでvまで潜らない場合の方がコストが良いです。よってL_i - E_i - D_v >= 0の条件を無視して解いても答えは変わりません。(ここがなかなかわからなかった、まあライブラリが強ければ思考停止でいいので)
結局分割統治の各stepではw_vが小さい順に直線を追加していくことで蟻本のenvelopeの管理と同様にdequeででき、O((NlogN + MlogM)logN)でこの問題が解けました。
ACM-ICPC 2018 Asia Yokohama Regional 参加期
チーム catsatmat で参加して ABCDGK の6完12位でした。アジア地区予選には2016年も参加していて11位だったので、今回の結果は残念なものでした。感想文だけ送るのが面倒なので簡単な参加記も残しておくことにしました。参加記を書き書き。
まず、チームは catupper と satashun と shiatsumat で構成されていて、まあ3人とも橙のどこかみたいな実力してると思います。
- コンテスト前
shiatsumat が来なくて焦った。maroonとcatupperがスマブラをしていて草。US配列に前日必死で慣れた。
- コンテスト
catupperがemacs、satashunとshiatsumatがgeditを使う予定だったのでとりあえずAをcatupperに任せてあとの2人は後ろを読む。
読んでないけど実装が怠いみたいなことを言っていた気がするがA AC(0:14)。
satashunはBを、shiatsumatはCを読んでいたがまだどちらも解けていなかった。
shiatsumatがCはsatashunの方が得意そうというので見たら、マンハッタン距離で個数をカウントして、小さい方から繰り上がらせていけばいいのではということになった。配列サイズが半分しかなくて1WAののち、C AC(0:27) 。
Cが通ったことで落ち着いたのか、続いてBも(今の要素, 前の要素)の組は1回ずつしか見なくていいことに気づき実装。B AC (0:38)。初めmathっぽく考えてしまった(初項を固定してみたりとか)
Dが一瞬で解かれていたのでshiatsumatに読んでもらいつつ後ろを見ていく。自分はJとかKとかを読んだがDをもらう。Gはcatupperが出来たということでG AC(1:08)
この辺まではそこまで悪くなかったのに詰まり始める。
しばらくして satashun さんがDは出来たと主張して実装してバグらせてくれました。アイディアは割とすぐ出たのに遷移がとてもバグってしまった(2つともを処理して終わったかどうかが適当だった) catupperはK出来たというのでパソコンの奪い合いのち、D AC (2:41), K AC (3:04)
Eは苦手なので2人に投げて、パソコンが空いている時間にJをしばらく考えていたらLCAが書ければできることに気付いた。
Jはsatashunさんが担当してバグらせてくれました。Eはcatupper、Fはshiatsumatが書いてた。
残りの時間はJを実装したりEを実装したりFを実装したりという感じだったがどれも通せず。Hもshiatsumatにより解法っぽいものが出てたが色のflipが詰めきれてなかったっぽい。
- 反省
Fは2人とも嘘をついたせいで重いだけに見えてしまったのがミスだった(EとJが手詰まりに見えたからなんですが、冷静に幾何なんてやるべきではないだろ)
Eは自分は見てないけど2人の解法のいいところを取るといい感じだったっぽくて残念。
Jはほとんど合ってたが、lcaが変わる場合の処理が甘かった。終盤時間が足りなくて手元でサンプルを作る余裕がなかった。
Hみたいなの残り2人は得意っぽいので前半を早くして余裕が作れたらチャンスなんだけどなあ(2年前もskinny polygonが終わらなかった)
まあでも学内topが離れすぎていたためとりあえず手を付けた問題を全部解き切るための動きだったので、人を適宜マージしたら1問は増えた(上手く行けば+2)かなあというのがあるが、シンプルに実力不足だと思いました。中盤が多いコンテストだったので少し実力が伸びればだいぶ変わりそうなセットでした。
来年も出られれば5年連続同じチームということになるのでいい結果が残せるといいですね(赤くなろうね)
有理数の探し方
この記事は Competitive Programming (2) Advent Calendar 2018 - Adventar の 12/11 分です(の予定でした)
連分数展開・Stern–Brocot tree には様々な面白い構造があり、たまにプログラミングコンテストでも出題される割にあまり知らなかった(かつ日本語の記事もあまり見たことがない)ので、紹介していきます。
(Farey 数列 という概念も聞いたことがあるかもしれませんが、 Stern-Brocot tree に似たものです)
簡単のため、以下では基本的に非負の有理数を扱います。
- Stern–Brocot tree
唐突ですが、Stern–Brocot treeとは以下のような木を指します (from wikipedia)
まず作り方がイメージしやすい方 (この方法はこの図のイメージと少し違いますが)
初め端にと が書かれています。(0以上の有理数からなる区間を表しています)
ここから、以下の操作を 繰り返して数列を作ります。
- 数列で隣り合う全ての有理数の組 と の間に を書く。
のとき明らかに なので、この列は常に(通常の順序で)昇順となっています。(とが隣り合う時、 という性質も保存されます)
各頂点が区間を持っており、ある区間を分割して出来る2つの区間が子となった2分木のようなものをイメージすると良いです。
この図から、分母が小さい方から有理数が列挙されていく様子を感じてもらえると思います。
次に詳しくこの木の性質を見る前に、この木と深い関わりのある 連分数展開 を紹介します。
- 連分数展開
ここではsimpleな表現のみを考えます。(以下のように、各部分の分子に1しか使えないもののことです)
有理数 x が連分数展開として と表されるとき ] と表記します。
例えば、 なので、] です
この連分数の作り方は Euclid の互除法に似た操作となっています。各stepで、を表したいとして、 (ただし、) と表せます。
有理数 を連分数展開して]となったとします。(k : 非負整数, a_0 : 非負整数, a_1, .. : 正の整数)
のとき、] ]ですが、最後が1であることを許さなければ一意に表せます。
(雑ですが、 は正かつ1以下(1となるのは、hogeが無かつa_i=1の時だけ)であるため、毎回整数部分を入れていくしかないので)
まずこの表現の場合、根は[1]です。
そして、] ()]の子は、Stern–Brocot treeにおいて]と}です。
小さい方を左の子にするのですが、これはその数列の長さの偶奇によって変わって、奇数の時は末尾を+1したものの方が大きく、偶数の時はその逆です。
逆に、この頂点の親はakから1を引いたもの(akが1になったら撤去してa_{k-1}に1を足す)ことによって得られます。
(連分数展開をexplicitな展開を書く)
Stern-Brocot tree には綺麗な性質があって、
- 二分探索木を成す
- 全ての既約な有理数が1度ずつ現れる
ある頂点の左子孫は全て自分より小さく、右の子孫は全て自分より大きいという性質を満たします。これも深さの偶によりますが、連分数展開で今見ているところ以降を操作すると、右下の部分が変わるだけで
これまた雑ですが、全部が網羅されていそうなことの雑な説明です。
まず深さdには全て異なる2^d 個の頂点があります。
を観察すると、深さdにある頂点ではこの和がd+1となっていることがわかります。
先頭のみ0が許される & (先頭でない)末尾は1が許されない ような 和が d+1 の個数を考えます。
d=0の時は[1]のみなので良いです。それ以外の時、oがd+1個あって、末尾と、末尾の1つ手前以外に仕切りを入れると列になるのでこれは 2^d 通りです。
導入が複雑(雑)でしたがこういった性質を利用すると、有理数で近似・有理数の数当てゲーム(例えば整数であれば普通の2分探索で解ける) などができます。
(近似は結構複雑なのでそのうち追記したいですが、気になる人はwikipediaのcontinued fractionを見てください)
- 例題
XVIII Open Cup named after E.V. Pankratiev. Grand Prix of SPb
E. Decimal Form
なる整数a,bを用いてa/bを計算したものを小数点以下18桁にしたものが10^4個まで与えられるので、それぞれa,bを答えてください
Stern-Brocot treeの上で探索したいのですが、実は下る向きを変える(左の次は右、右の次は左)とすると分母がだいたい倍になるので O(log max) 回しか方向転換せず、それぞれでは doubling してまっすぐ下りることで解けます。別の解法としては、小数点以下18桁の x が与えられるは a/b が[x-10^-19, x+10^-19) に入るということで、この両端を連分数展開すると近似の問題として解くこともできます。
JOI春合宿 2008 day3 Fraction
を満たす既約分数a/bのうち小さい方からk番目を求めよ。
kが小さいので、この二分探索木に小さい方から潜るだけで良いです。(知識...) (既約なものってそんなにないのでは(トーシェント関数の累積和を評価すればよい))
- ギャグ
Pythonのfractionsにはなんと分母のmaxを指定して近似する関数 (limit_denominator)があります。なんてこった。
- 参考
Spaghetti Source - Stern-Brocot 木
Grand Prix of SPb - Codeforces
WebDiarios de Motocicleta: Rational Search
Continued fraction - Wikipedia
https://www.ipsj.or.jp/07editj/promenade/4503.pdf
説明は自分で考えたところが多いので誤りがあったり、他にも面白い応用(や問題)があったら教えてください。
ICPC2018国内予選 H問題
今年もチームcatsatmatで出場しました。4位で通過出来て嬉しかったです。
- H問題
前原先生の記事
木上のナップサック問題
でだいたいわかりますが、コードを載せておきます。
まず、素直な解法として、dp[v][c] := vの部分木についてコストcで得られるmaxとおくとO(nkk)ですが、このままでは高速化が難しいです。
ポイントは、各頂点vの子たちを、vの必要量の昇順に並べておくと、ある子に潜ることが出来る時その手前の子たちも全部使えて、このことを利用するとdpの配列を持ちつつdfsすることができるという点です。これだけでオーダーが良くなるのは不思議な感じですね。これに加えて各頂点に使える上限があるので面倒になっています(このパートはスライド最小値をやると線形で解けます) よって時間計算量はO(Case * NK)です。100*100*10^5 = 10^9を8秒なので結構実装に気を使わないといけません。
#include <bits/stdc++.h> using namespace std; #define rep(i, n) for (int i = 0; i < (int)n; ++i) typedef vector<int> vi; typedef long long ll; typedef vector<ll> vl; const ll inf = -1e18; const int MN = 110; const int MK = 100010; inline void chmax(ll &x, ll y) { if (x < y) x = y; } int n, k; vi g[MN]; int h[MN], s[MN], lo[MN]; vl dp[MN]; int deq[MK]; void dfs(int v, vl vec) { vl &res = dp[v]; int l = 0, r = 0; for (int i = 0; i <= k; ++i) { if (vec[i] != inf) { ll x = vec[i] - (ll)i * s[v]; while (l+1 <= r) { int id = deq[r-1]; if (vec[id] - id * (ll)s[v] <= x) { --r; } else { break; } } deq[r++] = i; } if (l+1 > r) { res[i] = inf; break; } else { res[i] = vec[deq[l]] - (ll)deq[l] * s[v] + (ll)s[v] * i; if (deq[l] <= i - h[v]) { ++l; } } } for (int to : g[v]) { dfs(to, vec); vl &nx = dp[to]; for (int i = 0; i <= k; ++i) { chmax(vec[i], nx[i]); } l = r = 0; for (int i = 0; i <= k - lo[to]; ++i) { if (nx[i] != inf) { ll x = nx[i] - (ll)i * s[v]; while (l+1 <= r) { int id = deq[r-1]; if (nx[id] - id * (ll)s[v] <= x) { --r; } else { break; } } deq[r++] = i; } if (l+1 > r) { nx[i] = inf; break; } else { ll val = nx[deq[l]] - (ll)deq[l] * s[v] + (ll)s[v] * i; chmax(res[i + lo[to]], val + (ll)s[v] * lo[to]); if (deq[l] <= i - (h[v] - lo[to])) { ++l; } } } } } int main() { while (cin >> n >> k, n) { rep(i, n) dp[i] = vl(k + 1, inf); rep(i, n) g[i].clear(); rep(i, n) cin >> h[i]; rep(i, n) cin >> s[i]; for (int i = 1; i < n; ++i) { int p; cin >> p; --p; g[p].push_back(i); } for (int i = 1; i < n; ++i) { cin >> lo[i]; } rep(i, n) { sort(g[i].begin(), g[i].end(), [&](int a, int b) { return lo[a] < lo[b]; }); } vector<ll> vec(k + 1, inf); vec[0] = 0; dfs(0, vec); cout << *max_element(dp[0].begin(), dp[0].end()) << endl; } return 0; }