SuperCon 2022 予選参加記
はじめに
3 回目のSuperConです. 今回の問題はここにあります.
問題概要
のマス目それぞれに から までの数字が書き込まれており,座標 の数字は である.
時刻 には座標 にロボットがあって,あなたは時刻が 進むことに上下左右隣接するマスにロボットを進めることができる.(時刻 でのロボットの位置を とする.)
各時刻においての目標の色のパターン が与えられるので,ロボットの進路を好きに決め,コスト関数 を最小化せよ.
書いたアルゴリズム
DP をして復元する方針とダイクストラ (または 0-1BFS っぽいもの) をする方針があると思いますが,DPの方針を取りました. 単純に時刻とその時のロボットの場所をキーにしてそこまでのコスト関数を持てばよいです.
ただこのままだとおもしろくなかったので,SIMD 化して,各行の処理を並列して実行するようにしました. 問題文の実行環境の欄にある「2.8 GHz Intel Core i5」という情報だけだとどの世代かがわからなかったので,sse4.2 相当対応の CPU で動くようにしています.
また, のように座標をいい感じに変換することで,やや計算量が下がります.(各行 の回数が 回減る.)
実行時間
手元で ケースあたり平均 秒でした.(標準偏差 秒) 愚直な DP 解が ケースあたり 秒だったのでかなり高速だと思います.
提出コード
#include <immintrin.h> #include <algorithm> #include "sc1.h" inline int to_index(int x, int y) { return (x * SC_L + y); } inline __m128i get_rotate_perm(int n, int r, int l) { alignas(32) unsigned char a[l]; for (int i = 0; i < n; i++) { a[i] = (i + r + n) % n; } for (int i = n; i < l; i++) { a[i] = i; } return _mm_load_si128(reinterpret_cast<__m128i*>(a)); } inline int mod_l(int x) { if (x >= SC_L) { return x - SC_L; } else { return x; } } __attribute__((target("sse4.2"), optimize("unroll-loops"))) int main() { SC_input(); constexpr int SIMD_L = 16; __m128i rotate_right = get_rotate_perm(SC_L, -1, SIMD_L), rotate_left = get_rotate_perm(SC_L, 1, SIMD_L); alignas(32) unsigned char dp[SC_Nt + 1][SC_L][SIMD_L]; std::fill(dp[0][0], dp[0][SC_L], 127); constexpr int center = (SC_L - 1) / 2; dp[0][0][center] = 0; __m128i C[SC_L]; for (int i = 0; i < SC_L; ++i) { alignas(32) unsigned char tmp[SIMD_L]; for (int j = 0; j < SC_L; j++) { tmp[j] = SC_C[to_index(mod_l(i + j), j)]; } C[i] = _mm_load_si128(reinterpret_cast<__m128i*>(tmp)); } for (int i = 0; i < SC_Nt; i++) { __m128i min[SC_L]; for (int j = 0; j < SC_L; j++) { min[j] = _mm_load_si128(reinterpret_cast<__m128i*>(dp[i][j])); } for (int j = 0; j < SC_L; j++) { min[j] = _mm_min_epu8(min[j], _mm_shuffle_epi8(min[j], rotate_left)); } for (int j = 0; j < SC_L; j++) { __m128i abs = _mm_abs_epi8(_mm_sub_epi8(C[j], _mm_set1_epi8(SC_Cr[i + 1]))); __m128i dist = _mm_min_epu8(abs, _mm_sub_epi8(_mm_set1_epi8(SC_Nq), abs)); __m128i res = _mm_min_epu8(min[mod_l(j - 1 + SC_L)], _mm_shuffle_epi8(min[mod_l(j + 1)], rotate_right)); res = _mm_add_epi8(res, dist); _mm_store_si128(reinterpret_cast<__m128i*>(dp[i + 1][j]), res); } } int min_x = 0, min_y = 0, min = 255; for (int i = 0; i < SC_L; i++) { for (int j = 0; j < SC_L; j++) { if (dp[SC_Nt][i][j] < min) { min_x = i; min_y = j; min = dp[SC_Nt][i][j]; } } } unsigned char next_min = 255; for (int i = SC_Nt - 1; i >= 0; i--) { SC_ans[i + 1] = to_index(mod_l(min_x + min_y), min_y); int next_x, next_y; for (int s = 0; s < 2; s++) { int x = mod_l(min_x + (s == 0 ? SC_L - 1 : 1)); for (int j = 0; j < 2; j++) { int y = mod_l(min_y + (s == 0 ? j : SC_L - j)); if (next_min >= dp[i][x][y]) { next_x = x; next_y = y; next_min = dp[i][x][y]; } } } min_x = next_x; min_y = next_y; } SC_output(); }
Miller素数判定法とMiller-Rabin素数判定法
高校の文化祭で講義をしたスライドで,日本語でMillerテストについて詳しく触れている文献があまりなさそう (Miller-Rabinなら山のようにある) なので公開しておきます.
Barrett Reduction について考えたこと
Barrett Reduction とは
えびちゃんさんの記事 がわかりやすかったです. アイデアは固定された除数 に対して 適当な 冪である定数 を使って を前計算しておくことで, を で近似すれば除算を乗算とビット演算で置き換えられるという感じです. また, の時近似の誤差が高々 であることが知られています.(ACLの実装も (実質) を要求しています)
考えたことその1 : とは限らない場合
とします.この時 です. で,求めたい値が であることに留意すると,誤差が高々 になる十分条件は より であればよさそうということがわかります. よりも扱いやすい制約になった気がします.
考えたことその2 : 誤差を にする方法
前の章と同じようなことをすると, であればよさそうなので,こうなる条件を考えることにします.これの十分条件として があるので, の最大値を として となるように をとることとします. この時オーバーフローが若干心配になりますが, なので途中で出てくる数はそこまで大きくないことも確認できました.
実際この方法のほうが ACL の実装より幾分か高速な様です.( としたとき を で切り下げ除算した結果が常に なことを利用した高速化もしています.)
https://wandbox.org/permlink/2lTdpwYkMnUw3MkG
追記 とした方がやや楽です.
https://wandbox.org/permlink/ozHZS9fb8BleVpt2
参考文献
AtCoder ARC128 C - Max Dot の双対解
問題概要
整数 と数列 が与えられる.
非負実数列 が と を満たすとき, の最大値を求めよ .
解法
とします.(ただし とします)
を使って問題を言い換えると以下のようになります.
制約
最大化
これは線形計画問題なので双対を取っても最適値は同じです.
制約
最小化
うれしいことに最適化すべき変数が 2 つになりました.
を固定すると の最小値は で出せるので, で三分探索して を最小化すれば計算量は です.
おまけ 三分探索で最適解が出ることの証明
での制約を満たす の最小値を と書くことにします.
ここで は制約を満たすので,任意の を満たす について も制約を満たします.
での の最小値は であるので であり, は凸関数です.
よって も凸関数となるので,三分探索などで最小値が求められます.
supercon2021 予選
はじめに
問題
途中経過
6月2日
問題を読む。
判定問題を考える。
個のマスを頂点としてそれぞれのマスからどの操作を使えるかを前計算して辺と見なせば、辺の個数は自明になのでDFSなどを使えばで解ける。
ところが辺を情報を生成を適当にやると になってTLE。
ではどうするかというと辺が使える条件を「操作の途中で盤面からはみ出さない」かつ「操作の途中で障害物にぶつからない」と言い換える。前者は各操作についての変位のを持っておけばok。後者は操作の途中まで決めた時ちょうど障害物に当たるような場所を逆算して求めることで調べられる。なので、全体でになる。
これ割と難しくないですか?diff1700ぐらいあります。
判定問題を解いたら、リポジトリを立てて、テストケースを100個生成した。
テストケースを眺めていると操作に不自然に同じ文字が続いていたのでgenerator.cのコードを読んだ。原因を理解したが、問題にはあまり関係なさそうだったのでそのあとは一切気にしなかった。
6月3日
判定問題のコードを書き上げた。
到達したかの判定やある場所である操作が使えるかの情報を持つための配列にvector<vector<char>>を使うかvector<vector<bool>>を使うかで迷ったが結局charにした。ちなみにあとでbool版も書いたが遅かった。
ここで頂点や辺の数を集計したファイルを作った。あとで便利だった。
このファイルを作った目的の一つにグラフの疎性の確認があって、疎なグラフでDFSをたくさん回してもいいのではという期待があったが見事に打ち砕かれた。
密なグラフということで何回もDFSを再計算するのは避け、順方向の遷移だけで済むビームサーチ方針にしようかななどと考えた。
6月4日
DFSの再計算をせずに辺を削除する方法を考えた。
最短経路木に含まれるの辺を削除するとき、到達できることの必要十分条件はは頂点からへのパスがあることであることは明らかなので、それをうまく管理したいと思った。
動的木のアルゴリズムを調べてみたが応用できそうに見えないうえ、実装もしたくないのであきらめた。
6月5日
蜜蜂にgit/githubの使い方を教えた。
環境構築大変だよね、わかる。
6月6日
簡単な貪欲を書いてみる。
操作の個数をの状態から、その操作を追加したときに到達できるようになる頂点の個数が最大のものを選ぶ、というのを繰り替えすように実装した。
そのままのDFSでは複数回走らすには遅かったので、操作を追加する前の到達済みの配列を流用して、途中からDFSを回すようにすると多少早くなった。
6月7日
DFSを毎回計算して操作の良しあしを決めるのは効率が悪いので、以下の指標を導入した。
(操作を追加する前に到達できる頂点の数その操作によって増える辺のうち、にすでに到達出来ていて、 に到達できていないものの数
前者は終了判定の時に入手可能で、後者は辺の数のオーダーで求められるので、速度が上がる。
6月11日
visualizerを書いた。
少なくともC++で書くのは苦行で、JSはよく解らないので、Pythonで書くことにした。GUIを持っているものを作るときにはOpenCVやPyGameを使うのが好みなのだけど、ライブラリを使うのは避けたかったのでtkinterを使った。
visualizerを書いたことでアルゴリズムが改善したとかはなかったが、出力形式に関するバグの発見に使えたので書いてよかったと思う。
このバグを見つけた時はmainかvisualizerのバグかわからなかったが、蜜蜂先生にcheckerを書いてもらったことでバグの場所を特定できた。感謝。
6月12日
local checkerを書いた。
pythonだとshellの実行とか並列化がやりやすくて良き。 この時の結果は貪欲だけでケースだった。
6月13日
ビームサーチの実装をした。
素直に貪欲をビームサーチにした形で、貪欲の時の指標に沿って上何個かを保持する。 ビーム幅は30ぐらいしかない割にになった、すごい。
6月14日
終点の入次数/出次数に応じて辺のスコアを変えること思いついたので、いろいろ実験を始めた。
AHCで話題に上がっていて気になっていたoptunaで入次数/出次数の重みの最適化をすることにした。実行には10並列でも半日ぐらいかかる(遅い)。
ビームサーチの重複を避けるためのzobrist hashの存在をここで思い出したので実装した。スコアが少し上がった。でのhashではあったが意外と衝突しても平気そうだった。
6月15日
探索の序盤と終盤で評価値の計算速度を縮めた。
探索の序盤、終盤はすでに到達できると到達できないの候補が少ないので、少ないほうを選べば楽に全探索できる。
個の高速化でビーム幅を70に変えた。
また、時間を食うビームサーチで時間を使いきった対策に、最初に高速な貪欲出力をすることにした。
6月16日
実装も詰めに差し掛かる。
まずTLE対策に9.5秒を過ぎた時点でビーム幅を10にするようにした。ここの境界の決定には気を使っていて、ジャッジのCPUがi3-8109Uであることを予想し手元環境の約1.4倍ぐらいとみて、それに合わせて時間を決めた。
またビームサイズを増やした時のスコアの動きを見て、ビーム幅を500ぐらいに出来ないかなと考えていた。
6月17日
実質最終日。
ここで天啓に導かれ、bitsetを使うことを思いつく。
ボトルネックである操作のスコアの計算を高速化する方法を考える。
ます操作での移動量shiftと、頂点に到達したか(visited)、頂点で操作を使えるか(edge)をbitsetとして持つ。
すると求める値はとなり、これはの計算量となる。
ビーム幅500とまではいかないものの128にすることができた。。
-algo.txt
初めに、答えがYESかNOであるかを求めます。 入力を読み取った後、各マスについてどの操作が使えるかを判定します。このために各操作について、操作の中でy,x座標が最大でどれだけ増える/減るかをそれぞれ求めることで範囲外に出るかどうか、また操作を「途中まで進めた時ある障害物とぶつかる」⇔「ある障害物から操作の途中まで逆算した位置から操作をする」であることを利用し判定しています。この判定にはO(N^2 M)かかります。 その後、DFSの要領で各マスに到達できるかを判定しています。グラフの頂点をマス目の各マスに見立て、グラフの有向辺をマスXから1回の操作によって移動できるマスYに対し、XからYへ貼ることでできるグラフでDFSを行っています。このグラフの頂点の数はN^2、辺の数はO(N^2 M)なのでDFSで十分高速に判定できます。 探索ではビームサーチで操作を一つずつ追加して、もし全ての頂点に達成できたなら終了するという方針にしました。終了判定にはYES/NO判定と同様にDFSを用いましたが、すでに見た頂点を探索で使いまわすことで高速に処理しています。 ビームサーチでは、ある状態に操作を追加したときの評価関数は次のものを用いました。 (追加する前に到達できる頂点数) + (追加する操作の辺のうち、端の一方がすでに到達でき、もう一方がまだ到達できないものの個数) 第一項は終了判定の結果が流用できるため、第二項を高速に処理することを考えます。 まず到達できる頂点、または到達できない頂点の個数が少ないときは、少ないほうの頂点から出る/入る辺を全て調べることで効率よく処理できます。 上のどちらでもない場合はbitsetにより高速化しました。到達できる頂点をvisitedで、その頂点で辺が使えるかをedgeとし、これをbitsetとして管理し、操作の移動量をshiftとして書くと求める値は popcount(((edge & visited) << shift) & (\neg visited)) となり、これはワードサイズをwとしてO(N^2 / w)で処理できます。 また、ビームサーチにおいて使用する候補の操作を重複して含むことを防ぐために、zobrist hashをもちいてビームの情報を管理しています。 ビーム幅は最初は128ですが、9.5秒を超えた場合は10に設定しなおし、すぐに終了させます。
感想
終了後他の人を見て煮た方針がなくびっくりした。 最初の貪欲から劇的な改善というほどではなかったので個人的には心残りがあったが、解法ガチャはSRぐらいは引けているようで良かった。
追記(2021/06/23)
予選通過した、やったー!