2019-10-30

DAGの橋？

なんか初めて見たので。

DAG与えられます。それぞれの辺について、その辺を削除したときに、sからtにいけなくなるか答えてください(s, tは固定)

DAGを無向グラフとみなして橋を求めるのはダメです。

4頂点4辺で、辺がそれぞれ

を結ぶようなケースで辺(1,2)を除くと1から3に到達できなくなります。

これは、tから逆辺で到達できる頂点のみに注目してグラフを作り直し、無向グラフだとみなして橋を求めればいいです。

2019-10-30

サイクル基底・サイクルの個数について

チームでICPC練をしていたらサイクルの数を数える問題に遭遇した(Count Cycles: Asia Tsukuba Regional 2017)

サイクル基底を知らずに解けたけど(追記: 解けてなかった)サイクル基底を知らなかったのでついでに調べた。

サイクル基底

まず直感的な説明をします。サイクルを2つ書いて、辺のxorをとってください。辺のxorとは、2つのサイクルが同じ辺を共有している場合その辺を削除するということです。そうすると新しいサイクル(正しくはオイラリアンなグラフ。つまり全頂点の次数が偶数。非連結になっちゃうこともあるので)ができますよね。グラフのサイクル(正しくはオイラリアンな部分グラフ)をすべて考えたいときに、いくつかサイクルを用意して、それらのサイクルをxorすることによってそれを全部考えようというのがサイクル基底のアイデアです。では、どのようなサイクルの集合(基底)を考えればいいのでしょうか。

$n$ 頂点 $m$ 辺, $c$ 連結成分のグラフを考えます。全域森を考えると、 $m - n + c$ 本の辺が全域森に含まれません。これら全域森に含まれない辺のそれぞれについて、その辺と全域森の辺のみからなる閉路があります。これらの閉路の集合がサイクル基底になります。

使い道

サイクル基底のアイデアは上のとおりですが、実際に扱うときは、 $m$ 次元ベクトルで、 $i$ 番目の要素はそのサイクルが $i$ 番目の辺を含むとき1, そうでないとき0としたものとして扱うことが多いようです。

サイクルの個数カウント

一般にサイクル基底で適当にxorをとったものはオイラリアンになります。今回解いた問題では、サイクル(各頂点を1度しか通らないようなもの)を列挙するのですが、今回の問題ではサイクル基底の個数がたかだか16なので、実際にサイクルを(1 << 16)個を全列挙し、そのうち同じ頂点を一度しか通らないものをチェックするという方針で解けます(サイクルの全列挙には、縮約の前処理が必要です。計算量が指数なので一般のグラフだとできないと思う)。

stパスのカウント

stパスと(サイクル基底を適当にxorしたもの)もstパスになります。tex: XOR回廊みたいな、パスのコストが辺のコストのxorで決まる問題との相性がすごいです。XOR回廊はまんま、stコストの最大値を求めたいわけですが、これを知っているとxorのGauss Jordan Eliminationに持ち込むだけになります。

サイクルの個数(bitdpで)

上述したようにサイクル基底でもどうにかなりそうです。ただ、今回のように同じ頂点を1度しか通ってはいけないという制約の場合、bitdpで解くこともできます。ぼくたちはサイクル基底を知らなかったのでこうしました(追記: と思ったんですが、頂点数が16程度ならこの方針で解けるんですが、実際は頂点数は50程度になることもあって無理でした)。

具体的には辺を一本ずつ追加していきます。辺 $(u, v)$ を追加するとき、新たに増えるサイクルの個数は $(u, v)$ を通るサイクル、つまり追加前のグラフにおける $(u, v)$ パスの本数です。これはdp[S][v]: 頂点集合Sを回っていまvにいるときの場合の数を計算すればできるので、bitdpを辺の本数分やれば求まります。こっちのほうが思考停止でできるのでよさそうな気がします。

2019-10-28

LinkCutTreeを書きました

拾ったLinkCutTreeの使い方がわからないことがたまにあったので自分で一番使いやすいと思うものを書きました。

けんしんのLazySegmentTreeやぼくのImplictTreapとまったく同じインターフェースで使えます。

たまに便利なモノイドが追加されるのでそれが一気にセグ木にもTreapにもLCTにも乗るのはうれしすぎますね。

verify済みです。ソースコードはここからコピペしてください。

#394456 No.399 動的な領主 - yukicoder

使い方は以下のとおりです。

// 上のソースコードにある中ですきなタイプのクエリをとってくる
template <class T0, class T1>
struct MinUpdateQuery : public BaseLinkCutTree<T0, T1> {
    using BaseLinkCutTree<T0, T1>::BaseLinkCutTree;
    MinUpdateQuery() : MinUpdateQuery(numeric_limits<T0>::max(), numeric_limits<T1>::min()) {}
    T0 f0(T0 x, T0 y) override { return min(x, y); }
    T1 f1(T1 x, T1 y) override { return y == numeric_limits<T1>::min() ? x : y; }
    T0 g(T0 x, T1 y) override { return y == numeric_limits<T1>::min() ? x : y; }
    T1 p(T1 x, int len) override { return x; }
};

int main() {
    // ...
    MinUpdateQuery<int, int> lct;
    // 頂点追加やlink/cut, evertなどもいろいろできる。基本的にpublic関数しか触らなくてよくて、expose忘れとかevert忘れとかしにくいかも
    // uvパスにある頂点の値をxに変える
    lct.update(u, v, x);
    // uvパスにある頂点のminを求める
    lct.query(u, v);
    // ...
}

使いやすすぎないか？？？

2019-10-25

進化したImplicit Treapたち

はじめに

昨年のAdvent Calendarで書いた、Implicit Treapというとても有能なデータ構造があります。

xuzijian629.hatenablog.com

機能をもりもり盛り付けた平衡二分探索木で、セグ木でできるような操作が、配列の要素を反転したり削除したり挿入したりしながらできるデータ構造です。

去年のAdvent Calendarでは一応モノイドを乗せる抽象化をしたImplicit Treapの実装を公開したのですが、あれ以来いろいろ進化したので改めて公開します。記事後半にあるPrioritySumというのとPairQueryというやつが主に新しいです。

マクロ等を一切使っていないので特に衝突なく貼るだけで動くはずです。

プレーン

library2/implicit_treap.cpp at master · xuzijian629/library2 · GitHub

コードがかなり長いのでリンクで共有します。よくあるRMQ+RUQ, RSQ+RAQなどの組み合わせについては事前に定義済みなのでそのまま使えます。

けんしんのLazySegmentTreeとまったく同じ継承方法を採用しています。こっちも便利なので宣伝しておきます。

サンプルコード(抜粋)

MinUpdateQuery<int, int> hoge;
hoge.set_by_vector(vector<int>(n, (1LL << 31) - 1));

hoge.query(l, r);
hoge.update(l, r, x);

リファレンス

int size()

現時点でのサイズを返します。 $O(1)$

void insert(int pos, T0 x)

先頭からposの位置にxを挿入します。たとえばpos = 0なら先頭に挿入します。 $O(\log n)$

void update(int l, int r, T1 x)

[l, r)の半開区間にxを作用させます。 $O(\log n)$

T0 query(int l, int r)

[l, r)の半開区間について累積を求めます。 $O(\log n)$

int binary_search(int l, int r, T0 x, bool left = true)

[l, r)内のkでf0(tr[k], x) != xとなる最左/最右のものの位置を返し、存在しない場合は-1を返します。 $O(\log n)$

説明が非直感的かと思うので具体例で説明すると、累積用の演算がminの場合、[l, r)の範囲にあるx未満の最左/最右の要素の位置を返します。

void erase(int pos)

位置posの要素を削除します。 $O(\log n)$

void reverse(int l, int r)

区間[l, r)を反転します。 $O(\log n)$

void rotate(int l, int m, int r)

区間の[l, r)の先頭がmになるようにシフトさせます。std::rotateと同じ仕様です。 $O(\log n)$

T0 operator[](int k)

インデックスアクセスができます。 $O(\log n)$

void dump()

デバッグ用です。配列の中身をprintします。

PrioritySum

これを含め、以下のライブラリはImplicit Treap(プレーン)に依存するので、そちらも一緒に貼って使ってください。 library2/priority_sum.cpp at master · xuzijian629/library2 · GitHub

ネーミングはびーとくんのライブラリから拝借していますができることは微妙に違います。

以下の操作が $O(\log n)$ でできます。

要素の挿入
要素の削除(最大/最小要素でなくてもいい)
上位k個の累積(このkが固定だったらびーとくんのでいい。定数倍軽い)

リファレンス

void add(T a)

要素を追加。 $O(\log n)$

void erase_at(int k)

位置を指定して要素を削除。先頭の削除ならk = 0を指定。 $O(\log n)$

void erase_value(T a)

要素を指定して削除。aがすでにaddされていることが仮定されている。 $O(\log n)$

int size()

要素数を返す。 $O(1)$

T sum(int k)

先頭k要素の累積を返す。 $O(\log n)$

T operator[](int k)

インデックスアクセスができます。 $O(\log n)$

void dump()

デバッグ用です。配列の中身をprintします。

PairQuery

ネーミングは適当です。私事ですがこの名前のライブラリを何種類かもっているので(機能は違う)名前を変えたい…

library2/pair_query.cpp at master · xuzijian629/library2 · GitHub

以下の操作が $O(\log n)$ でできます。

ペアを挿入
ペアの第一要素がx以下の要素全体に対し、その第二要素の累積を求める

リファレンス

void add(const pair<T0, T1> &a)

挿入します。 $O(\log n)$

T1 query(T0 x)

第一要素がxの要素とxより左側にある要素に対し第二要素の累積を返します。 $O(\log n)$

あとはpriority_sumと似たような機能もあります。

おわり

月1ぐらいの頻度で使える。けっこうすごい。

2019-10-24

Randomized Algorithms for Computational Geometry

とりあえずいくつか面白い論文をみつけた

Randomized Incremental Construction of Delaunay and Voronoi Diagrams

Randomized Algorithms for Geometric Optimization Problems

An Introduction to Randomization in Computational Geometry

ひとつ読んだのでまとめる。

$O(n\log n)$ randomized algorithm for slope selection

問題

二次元平面上に $n$ 点ある。このうち2点を選んで直線を引く時、 $1 \le k \le {}_nC_2$ 番目に小さい傾きを求めよ。

双対問題

点 $(a, b)$ を直線 $y = -ax + b$ に写し、直線 $y = \alpha x + \beta$ を $(\alpha, \beta)$ に写す操作を考えると、問題は以下と等価になる。

2次元平面上の $L$ 本のnon-verticalな直線について、それらの交点のうち、 $x$ 座標が $k$ 番目に小さいものを求めよ。

f:id:xuzijian629:20191024143431p:plain

簡単のため、どの3直線も1点で交わらないとし、どの2直線同士も同じ $x$ 座標で交わらないとする。

考察

答えの $x$ 座標が $(\alpha, \beta)$ の区間に存在するとする。このとき、 $y$ 軸に平行な帯 $(\alpha, \beta) \times \mathbb{R}$ を考え、 $x = \alpha$ における直線の順序と $x = \beta$ における直線の順序を考える。前者の順序を $1, 2, \ldots, n$ とし、後者を $\pi(1), \pi(2), \ldots, \pi(n)$ とすると、 $i\lt j$ かつ $\pi(i) \gt \pi(j)$ のとき、またそのときに限り2直線は交わる。つまり区間内にいくつ交点があるかはこの順列の点灯数を求まることで決まる。

これを何度か繰り返して区間を狭めていくことができるが、乱択と組み合わせるともうすこし賢くできるらしい。

2019-10-21

龍が如く劇場版 (2007) 感想

YouTube Moviesでおすすめに出てきたので300円払って見た。

www.youtube.com

ストーリーは、初代龍が如く、およびそのリメイクである龍が如く極に沿った形になっているがだいぶショートカットされている。

高評価58低評価13ですでに心配だったが見終えてみると評価に納得した。

初めに感想をまとめておくと、ストーリーを知らない人はこの映画を見ないほうがいいと思うし、この映画で龍が如くを知ってほしくない。ゲームのストーリーは終盤かなり感動要素が多いが、感動要素が皆無になっている。

見たことない人に対するおすすめ度は0です。終了。

以下ストーリーを知っている人向けの感想。

真島さんの演者がすごい！！！！！！なんか意味不明な感じを完璧に表現していてすごかった。というか真島さんの尺が長すぎるｗ終盤のミレニアムタワーのところが数分しかないんだけど、序盤と中盤が真島さんなのでもはや真島さんの映画になっている。

それと龍が如くのゲーム内ではおなじみのヒートゲージだが、映画でもエフェクトでヒートゲージがついていてさすがに笑ってしまった。終盤で、ドリンクを飲んで体力全回復 && ヒートゲージ維持するシーンがあるんだけど、錦山組の組員も桐生さんがドリンクを取り出した瞬間、「ああっ」みたいなリアクション(スローモーション)をしていて、ドリンクを飲むと最強になることが暗黙のルールみたいになっていてさすがに草だった。というかあそこの前後原作だとかなりの感動シーンなのにドリンクにすべてもっていかれた。

最後があまりにあっけなくて、牛沢さんの極実況の最後の3話ぐらいを見るほうがストーリーはいいレベルなんだけど(でも実況見るときは最初から全部見てね！)、ギャグとしてはかなり面白かった。とくに最後の10分ぐらいはかなり笑っていた(笑う話じゃないんだけど…) なんか、映画製作陣がどういう気持ちで作っているかはわからないんだけど、ぼくなりの解釈では、演者さんたちがストーリーをすでに知っている龍が如くファンを楽しませるために精一杯やってくれたように感じました。ストーリーを伝えるというよりも、ファンサービス的な。まあじゃないと最後のシーンあんな数分で終わらせないよね。

ということで、見たことある人向けのおすすめ度は5段階評価で3か4ぐらいです。牛沢さんがこれ見ながらツッコミを入れてる姿を見たい(無理だろうけど)

2019-10-15

Python(特にNumPy)の配列操作

紹介

NumPyの配列操作の仕組みについて紹介します。生のPythonよりもできる操作が多いのでこちらで統一します。さて、NumPyで

a = [1,...,n]
b = a[::-1]

としたとき、2行目の操作は $O(1)$ か $O(n)$ かご存知ですか？まあこれは $O(1)$ です。他にもNumPyは配列に関する多くの操作が $O(1)$ でできるよう配列のデータの管理が工夫されています。この記事ではいくつか具体例をあげて説明していきます。

知っている人向けに書くと、shapeとstridesの関係を書くだけです。ここらへんのattributesは開発者なら毎日のように使っていそうだけど一般ユーザはあんまり知らなそう(特にstrides)

shapeとstrides

多くのプログラミング言語で配列のデータは連続領域に確保されます。すなわち、a[k]とa[k+1]は隣り合ったメモリ領域に確保されます。しかし、NumPyの配列は必ずしも個々のデータを連続領域にもちません。NumPyの配列は、data, shape, stridesという3要素によって管理されます。

data

データの先頭のポインタです。

array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a.data
<memory at 0x110417948>

shape

データの多次元サイズ情報です

>>> a.shape
(10,)

strides

データアクセスのindexが1増えるとメモリ領域が何byteずれるかを表します。

>>> a.strides
(8,)

これは、配列aがnumpy.int64の型なので、a[k]とa[k+1]が8byteずれていることを表します。numpy.int32型に変換してみると、stridesが変わることが確認できます。

>>> a.astype(np.int32).strides
(4,)

多次元配列の場合の例も見てみましょう。

>>> a
array([[0, 1, 2],
       [3, 4, 5]])
>>> a.shape
(2, 3)
>>> a.strides
(24, 8)

これは、この場合データがメモリ領域に連続して0,1,2,3,4,5と並んでおり、二次元配列の[i,j]番目の要素が連続領域の3*i + j番目の要素に対応するからです。このときjを一つずらすと、int64 1つ分の8byteがずれ、iを一つずらすとint64 3つ分の24byteがずれます。

要素のアクセス方法

例として3次元配列で説明します。a[i][j][k]にアクセスするためには、メモリのどこを見ればいいでしょうか？

先程の例と同じことですが、C言語風に答えを書くと*(a + strides[0] * i + strides[1] * j + strides[2] * k)です。

a[::-1]の表現方法

さて、ここまで説明をして、stridesいらなくないかと感じた人が多いと思います。実際、データが連続している場合、 $\mathrm{strides}[k] = データサイズ * \prod_{i = k + 1}^{n} \mathrm{shape}[i]$ が成り立つからです。

自分でa = [1,2,3]にように配列を宣言する場合や、np.arangeなどによって宣言する場合についてこれは正しいですが、適宜shapeやstridesを変えることで同じデータをもとに幅広い配列を表現できます。いくつか例を紹介していきます。

たとえばa[::-1]のstridesを見てみましょう。

>>> a = np.arange(10)
>>> a.strides
(8,)
>>> a[::-1].strides
(-8,)

a[::-1]はメモリ上のデータ配置を変えずにstridesのみを変えることによってメモリ上のデータを後ろ向きに走査しています。データは変わらないので配列の宣言時に新しいメモリの確保は行われず $O(1)$ です。

broadcastの表現方法

NumPyの嬉しさの一つとしてブロードキャストがあります。例えば以下のような演算ができます。

>>> a = np.arange(6).reshape(2,3)
>>> a
array([[0, 1, 2],
       [3, 4, 5]])
>>> b = np.arange(3)
>>> b
array([0, 1, 2])
>>> a + b
array([[0, 2, 4],
       [3, 5, 7]])

このとき、内部的に一度bをaと同じサイズに拡張して計算しています(np.broadcast_to)。

>>> a
array([[0, 1, 2],
       [3, 4, 5]])
>>> b
array([0, 1, 2])
>>> c = np.broadcast_to(b, (2,3))
>>> c
array([[0, 1, 2],
       [0, 1, 2]])
>>> a + c
array([[0, 2, 4],
       [3, 5, 7]])

この拡張もstridesを変更するだけで済むので $O(1)$ です。

>>> b.strides
(8,)
>>> c.strides
(0, 8)

言われれば当たり前ですが、対応するstridesを0にすると、そのindexを変えてもアクセス位置が変わらないということであり、実質次元を拡張するような操作になっています。

transposeの表現方法

配列のtranspose(軸の入れ替え)も実は非常に簡単に表現できます。

>>> a = np.arange(6).reshape(2,3)
>>> a
array([[0, 1, 2],
       [3, 4, 5]])
>>> a.strides
(24, 8)
>>> a.T
array([[0, 3],
       [1, 4],
       [2, 5]])
>>> a.T.strides
(8, 24)

このように、strides(とshape)の要素の順序を入れ替えるだけです。メモリ領域には触らないので $O(1)$ です。参照渡しは最高ですね〜〜〜

もうちょっと詳しく

普通に読み飛ばしてもらって結構です。

transposeというと2次元配列の行と列を入れ替えるイメージが強いと思いますが、numpyでのtransposeはもっと幅広く、一般次元の配列(テンソル)について、次元の順序を任意に入れ替える操作のことを指します。

配列の中身は省略しますが

>>> a = np.arange(12).reshape(2,2,3)
>>> b = np.transpose(a, (1,0,2))

のように扱えます。np.transposeの第2引数には(0, 1, ..., len(a) - 1)のpermutationが入ります。

これはなかなか便利で、たとえばnp.sumのような演算(特定の軸について潰すような演算)の実装は、配列の初めのk軸を潰すような演算として実装されていますが、np.sum(a, axis=1)のような演算が来たとしても、初めに軸を回転させて1番目のaxisを0番目に持ってくることにより、そうした実装を適用することが出来ます。transposeはこうしたリダクション演算の実装には欠かせない機能になっています。

indexing

先程のセクションに比べてちょっとレベルが下がります。

たとえば要素を2つおきにとってきた配列を作るとき、a[::2]としますよね。これもstridesを変えるだけで行えるので $O(1)$ で生成できます。

>>> a = np.arange(10)
>>> a[::2].strides
(16,)

応用編

とりあえずここまででも、計算に便利な配列を $O(1)$ で生成できるうれしさはまあまあ伝わっていると思いますが、一つ本質的に役に立つ例を紹介します。

多次元配列のindexアクセスではa[i,j,k]要素にアクセスするのにsum(a.strides * numpy.array([i, j, k]))のoffsetだけデータのポインタからずれた値を見る、と言いましたが、配列のアクセスのたびにこのような演算をしていてはオーバーヘッドが非常に大きくなってしまいます。そこで、配列の次元圧縮というものを考えます。

そのために、データがcontiguousかどうか、というフラグを考えます。これはnp.ndarray.data.contiguousに対応します。これは、データがメモリ上で連続しているか、つまりa[i][j][k]を(i,j,k)の辞書順に見た時、すべての要素が等間隔で並んでいるか、のフラグを指します。

>>> a = np.arange(10).reshape(2,5)
>>> a
array([[0, 1, 2, 3, 4],
       [5, 6, 7, 8, 9]])
>>> a.data.contiguous
True
>>> b = a[:,::2]
>>> b
array([[0, 2, 4],
       [5, 7, 9]])
>>> b.strides
(40, 16)
>>> b.data.contiguous
False

例えば上のように生成した配列bは(0,2)番目と(1,0)番目の要素が8byteしか離れていないため(ほかは16byteだけど)contiguousではありません。

さて、配列の次元圧縮とは、配列がcontiguousの場合に1次元の配列に潰すことです(より広義に、配列の一部の次元がcontiguousである場合にその次元を潰す、ということもあります)。

1次元の配列のindexアクセスは高速なので、実質1次元配列の多次元配列に関しては1次元配列のようにアクセスしてしまおうということです。

とくに要素ごとの足し算(np.ndarrayに対するa + b)をイメージしてもらえると、a, bがcontiguousな場合に計算が楽になることがわかると思います。

応用編(hoge)

実はnp.ndarray.stridesはreadonlyではないので自分で勝手に書き換えて遊ぶことができます。

>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a.strides
(8,)
>>> a.strides = (4,)
>>> a
array([          0,  4294967296,           1,  8589934592,           2,
       12884901888,           3, 17179869184,           4, 21474836480])

とくに応用することはなさそうですけど遊んでみると楽しいです。メモリを生で触っている感覚を楽しめるかもしれないです。