2019-10-18

mod 逆元と拡張ユークリッド互除法

概要

本稿では整数の $\bmod M$ での逆元について取り扱います。
文中の合同式は $\bmod M$ です。

拡張ユークリッド互除法で逆元を求める
互除法の手順を高速化できるかもしれないという提案
$1$ から $N$ の逆元を $\Theta ( N )$ で計算する
$N$ 個の数について全ての逆元を $\Theta ( N + \log M )$ で計算する (おまけ)

拡張ユークリッド互除法で逆元を求める

$a$ の逆元を計算します。以下の形の合同式を考えます。
$s \equiv ta$
もし $s = 1$ なら $t$ が求める逆元となります。

ここで、明らかな合同式として以下の $2$ つがあります。
$M \equiv 0a$
$a \equiv 1a$
この $2$ 式を足し引きすると、新たな合同式が得られます。例えば
$M - 2a \equiv (0 - 2)a$
などです。よって、 $M$ と $a$ でユークリッドの互除法を行い $t$ も同時に管理すると、 $s = 1$ となり逆元が求まります。

例 $: M = 7 , a = 5$
$7 \equiv 0 * 5 \qquad \cdots ①$
$5 \equiv 1 * 5 \qquad \cdots ②$
$2 \equiv 6 * 5 \qquad ( ① - ② ) \qquad \cdots ③$
$1 \equiv 3 * 5 \qquad (② - 2 * ③)$
$\therefore 5 ^{-1} \equiv 3$

互除法のショートカット

$s \equiv ta$ について、 $s$ の逆元 $s ^ {-1}$ が既に分かっていたとします。すると両辺に $s ^ {-1}$ を掛けることで $1 \equiv ( s ^ {-1} t ) a$ となり、 $a$ の逆元が得られます。
普通は $s = 1$ となったときに互除法を停止しますが、例えば $1$ から $N$ までの逆元を前計算していた場合、 $s \le N$ となった時点で終了できることになります。 $s$ はおおまかに指数的に減少するため、 $N = \lceil \sqrt M \rceil$ まで計算しておけば逆元が $2$ 倍速程度になることが見込めるかもしれません。

$1$ から $N$ までの逆元を $\Theta ( N )$ で計算する

互除法を $1$ 手順だけ進めたときを考えてみます。
$M \equiv 0a$
$a \equiv 1a$
$\displaystyle M \bmod a \equiv - \lfloor \frac M a \rfloor a$
前述したように、もし $( M \bmod a ) ^ {-1}$ が計算済であれば $\Theta ( 1 )$ で逆元が得られます。ここで $M \bmod a \lt a$ であるため、 $a$ の昇順に逆元を計算することでこれは達成できます。よって全体で $\Theta ( N )$ の時間計算量となります。

$N$ 個の数について全ての逆元を $\Theta ( N + \log M )$ で計算する

数列 $a _ 0 \ldots a _ { N - 1 }$ のそれぞれについて逆元を計算します。
先頭からの累積積 $p _ i := \prod _ { k \lt i } { a _ k }$ と末尾からの累積積 $s _ i := \prod _ {i \lt k} {a _ k}$ を $\Theta ( N )$ で計算します。
次に、全要素の逆元の積 $c := \prod {a _ i ^ {-1}}$ を $p _ N ^ {-1}$ によって $\Theta ( \log M )$ で計算します。
すると、 $a _ i ^ {-1} = c p _ i s _ i$ によって $\Theta ( 1 )$ でそれぞれの逆元を計算可能です。

2019-10-10

スライド最小値を半群に一般化する

データ構造

概要

以前記事を書いたのですが消したので別の場所にもう一度書きました。今度はある程度まともな実装付きです (その分遅いですが)。

Sliding Window Aggregation - data-structures

計算量解析

Queue の場合

ポテンシャル関数を $\lbrace 末尾側の \ \rm stack \ の要素数 \rbrace$ で定義します。

Deque の場合

ポテンシャル関数を $\rm abs ( \lbrace 先頭側の \ \rm stack \ の要素数 \rbrace - \lbrace 末尾側の \ \rm stack \ の要素数 \rbrace )$ で定義します

2019-10-08

計算量について、償却/期待/平均など

アルゴリズム

本記事は皆様からのご指摘を募集しております

誤った記述があるかもしれません

概要

競技プログラミングをやっていると ${\rm average} \ \mathrm{O} ( N )$ などの表記を見掛けることも多いでしょう *1。 ${\rm best} , {\rm worst} , {\rm average} , {\rm expected} , {\rm amortized}$ それぞれについて、大雑把な意味をまとめました。アルゴリズムの挙動の正確な把握は競技においても重要です。
以降、全て時間計算量に付いて議論します。

注: 本稿内で用いられる $\mathrm{O}$ はほとんどが $\Theta$ に置き換えられますが、 Big O notation と同時に説明すると混乱を招くと判断し、競技プログラミングにおいて常用されている $\mathrm{O}$ を使用しています。

$\rm best :$ 最良計算量

多くのアルゴリズムは、入力によって計算量が変化します *2。例えば、ソートアルゴリズムには大まかに $N!$ 通りの入力が存在します。あり得る全ての入力のうちの計算量の最小値を最良計算量と呼び、 $\rm best$ を付けて表記します。

線形探索は ${\rm best} \ \mathrm{O} ( 1 )$ (最初に求める値が存在した場合)
マージソートは ${\rm best} \ \mathrm{O} ( N \log N )$
挿入ソートは ${\rm best} \ \mathrm{O} ( N )$ (実装次第)

例えば、ほとんどのソートアルゴリズムは最初に列がソート済かを調べることで ${\rm best} \ \mathrm{O} ( N )$ に改善できるでしょう。 (そのような改善をして嬉しいかどうかは別の話ですが。)

$\rm worst :$ 最悪計算量

最小値を取った最良計算量に対して、最大値を最悪計算量と呼び、 $\rm worst$ を付けて表記します。

線形探索は ${\rm worst} \ \mathrm{O} ( N )$
マージソートは ${\rm worst} \ \mathrm{O} ( N \log N )$
挿入ソートは ${\rm worst} \ \mathrm{O} ( N ^ 2 )$

最悪計算量は強い保証であり、最悪計算量でアルゴリズムを見積もれば想定外に時間がかかることはありません。競技プログラミングはどのようなテストケースがあるか分からないため、最悪計算量が重要となることが多いです。特に何も書かない場合、暗黙に最悪計算量を指すことがあります。(文脈に依存します)

$\rm average :$ 平均計算量

全ての入力に対しての平均を平均計算量と呼び、 $\rm average$ を付けて表記します。

先頭の値をピボットに用いるクイックソートは ${\rm average} \ \mathrm{O} ( N \log N )$

平均的に高速なアルゴリズムは、最悪計算量が悪くても実用的な場合があります。競技プログラミングで平均計算量を信用することは必ずしも得策ではありませんが、「テストケースがランダムに生成されている」場合解法を平均計算量で評価することが可能です。逆に writer に想定されていると、平均計算量のよいアルゴリズムも最悪計算量が悪いと落とされてしまいます。

$\rm expected :$ 期待計算量

アルゴリズムの中には、実行時に乱数を生成し、乱数の出目次第で計算量が変化するアルゴリズムがあります。このとき、計算量の期待値 (平均値) を期待計算量と呼び、 $\rm expected$ を付けて表記します。平均を取るという点で平均計算量と似ているように見えますが、異なる概念です。

ピボットをランダムに選ぶクイックソートは ${\rm expected} \ \mathrm{O} ( N \log N )$
ハッシュテーブルの検索は ${\rm expected} \ \mathrm{O} ( 1 )$ *3

平均計算量との違いは、平均を取る対象が入力全体か、乱数全体かという点です。入力はしばしばランダムではないですが、乱数は常にそのランダム性が保証されています *4。よって競技プログラミング的視点では「平均計算量は hack 可能」「期待計算量は hack 不可能 *5」という特徴があります。ピボットをランダムに選ぶことで、クイックソートが期待計算量に改善されている点は興味深いと言えるでしょう。

$\rm amortized :$ 償却計算量

償却計算量とはならし計算量とも呼び、データ構造のように複数の操作を繰り返し行う場合に操作一回を解析する概念です。
例として動的配列 *6 に要素の追加を繰り返し行った際の挙動を考えてみましょう。容量が限界に達したら $\mathrm{O} ( N )$ 掛けて倍の大きさのメモリを確保し、全ての要素を移動させます。拡張のタイミングの操作は当然 $\mathrm{O} ( N )$ の計算量が掛かります。しかしこのような「重い」操作は滅多に起こらず、実際には一操作辺り平均して $\mathrm{O} ( 1 )$ の計算量しか掛かっていません。
この「平均すると速い」は平均計算量や期待計算量に近い性質ですが、どのような入力でも大丈夫ですし、運の良し悪しにも左右されません。償却計算量はこのようなアルゴリズムを評価することに向いています。
(データ構造が空の状態から) $N$ 回の操作を行った際の計算量が $f( N )$ であるとき、一操作辺りの償却計算量は ${\rm amortized} \ \frac{f( N )}{N}$ となります。

動的配列への要素の追加は ${\rm amortized} \ \mathrm{O} ( 1 )$
splay tree の各種操作は ${\rm amortized} \ \mathrm{O} ( \log N )$

償却計算量は期待計算量と同様、競プロで信頼できる計算量評価の一つです。最悪計算量は悪いが償却計算量は良いようなデータ構造には、最悪計算量が同程度のものより簡潔、高速である場合があり注目に値します。

2021/12/19 追記

上記の定義は不完全で、不都合なことが多いと感じます。償却計算量について厳密に定義した文献を見つけられていませんが、私の経験上以下のような定義だと多くのケースで整合すると思われるものを紹介します。

データ構造に対して、起こり得る操作全体を $M$ とする。操作の償却計算量が $(A _ m) _ {m \in M}$ であるとは、任意の操作列 $(m _ i) _ {i = 1} ^ {n}$ に対して

\begin{align} \text{操作列} (m _ i) _ {i = 1} ^ {n} \text{を実行する実際の計算量} \leq \sum _ {i = 1} ^ {n} A _ {m _ i} \end{align}

が成り立つことと定義する。ここで、操作列 $(m _ i) _ {i = 1} ^ {n}$ は「何もない」状態から始まるもののみを指す。つまり、最初の操作 $m _ 1$ は「空の heap を作る」や「与えられた列から、segment tree を構築する」などになる。 (この条件が無いと、重い操作 1 つからなる列の計算量が壊れてしまう)

(2021/12/19 追記終わり)

余談

期待計算量や償却計算量は最良/最悪/平均計算量と同時に存在しうる概念だと思うのですが、厳密な定義や使用例が発見できなかったので詳しく書くことは控えました。ご存知の方がいらっしゃいましたら教えて頂けると嬉しいです (記事の誤りなども)。

*1: $\rm expected$ などの修飾を計算量の前と後のどちらに付けるかはどちらの例も見られよくわからなかったため、前で統一しています

*2: $\mathrm{O}$ 記法を用いるのは「計算量のオーダー」であって、「計算量」は計算回数を指すことに注意してください

*3:一部のハッシュテーブルは $\rm worst$ で保証する

*4:現実的には疑似乱数ですが、実用的には十分にランダムです

*5:乱数の seed を hack する場合は別

*6:C++ ならば std::vector

2019-10-05

区間代入/区間積 Θ(logN)/query

データ構造

概要

遅延セグメント木は様々な区間操作が $\Theta ( \log N )$ で可能です。一方で、代入と積の組み合わせは累乗を計算する必要があり、そのままでは $\Theta ( \log ^2 N )$ に悪化してしまいます。
本稿では空間計算量を $\Theta ( N )$ に抑えたまま、 $\rm amortized \, \Theta ( \log N )$ で処理する手法を提案します。また、積に限らず一般のモノイドについても同様に $\Theta ( \log N )$ 回の演算で区間への代入と和を取る処理が可能です。

アルゴリズム

基本的には遅延セグメントツリーと同様です。
$x$ を代入するクエリの時、 $N = 2 ^ {k}$ として $i \in \lbrack 0 , k \rbrack$ について $x ^ {2 ^ {i}}$ を計算し、保存しておきます。セグメントツリーはノードの要素数が $2^i$ の形しか存在しないため、これで $\Theta ( \log N )$ でクエリの処理が可能になりました。
さて、このまま愚直に処理をすると空間計算量が $\Theta ( N + Q \log N )$ になってしまいます。ここで $Q \gt N / \log N$ となったときに全ての遅延を解消し、前計算をすべて破棄します。この処理は $\Theta ( N )$ で可能なため、全体でも 1 クエリの計算量は $\rm amortized \, \Theta ( \log N )$ で保たれていることになります。

実装

https://noshi91.github.io/Library/library/gomi/assign_segment_tree.cpp.html

2019-10-03

2-SAT のアルゴリズムの証明

アルゴリズム

概要

蟻本で 2-SAT のアルゴリズムを最初に見たとき、私はその正当性について何ら疑問を抱きませんでした。これは私の頭がよく回ったからではなく、単に行間に気づいていなかったからということに最近気づきましたので、反省を兼ねて記録します。

反例のようなもの (反例ではない)

まず $x \lor y \Leftrightarrow (\lnot x \Rightarrow y) \land (\lnot y \Rightarrow x)$ に辺を張る部分を考えます。そもそも $x \lor y \Leftrightarrow (\lnot x \Rightarrow y)$ であり、 $\lnot x \rightarrow y$ へ辺を張るだけで十分なように思われます。 (これがまずい理由がわかる方はこの記事を読む意味がありません、申し訳ありません)

$x \lor y$ のみからなる 2-SAT のグラフを書いてみます。 $\lnot x \rightarrow y$ だけに辺を張ると、例えば以下のような順序のトポロジカルソートが考えられます。

$x,y$ の真偽を決めてみます。 $x \lt \lnot x$ なので $x = 0$ であり、同様に $y = 0$ となります。ところがこの割り当ては $x \lor y$ を満たしていません。これは「 $x$ から $\lnot x$ に到達可能なら $x = 0$ 」が必要条件であって、十分条件ではないことに起因しています。トポロジカルソートを一つ選んだ時対応する割り当てが定まることから、その割り当てが全体で整合することは (直ちには) 導かれません。
一方で $\lnot y \Rightarrow x$ にも辺を張るとこのようなことは発生しません。裏を返せば 2-SAT の正統性の証明は、明示的あるいは暗黙的に「 $x \rightarrow y$ の辺があるなら $\lnot y \rightarrow \lnot x$ の辺がある」事実を使用していることになります。

証明

トポロジカルソートして大小で真偽を割り当てたとき、それが全体で整合することを示します。

不整合を仮定すると、 $x = 1 \land y = 0$ なる $x,y$ であって $x \rightarrow y$ の辺を持つものが存在します。このとき $x \le y$ です。また、 $x,y$ の真偽の割り当てから $\lnot x \lt x$ 及び $y \lt \lnot y$ が分かります。ここで $x \rightarrow y$ より $\lnot y \rightarrow \lnot x$ の辺が存在しますが、 $\lnot x \lt x \le y \lt \lnot y$ と矛盾します。

よって割り当ては整合します。

まとめ

理解したつもりは怖い。思考の短絡には気を付けましょう

2019-09-29

JSC2019 参加記

A を読む　解けない
B を読む　解けない
C を読む　セグ木で丁寧だけど L/2 付近の挙動が良く分からなかったので放置
D を読む　解けない
E を読む　パターンマッチング系なのでどうせ S か T で Trie 作って、Aho Corasick すると終わり　ひいひい言いながらライブラリ漁ったら WA　丁度 1 時間が経過していて死にました
B を読むと解けるので解く
A を読むと NTT で、ほんまか？と言いながら書くと 1case WA　mod が悪かったかなと思って mod 変えると WA　N と M を間違えていた　通す
D を式変形したら CHT になったのでライブラリ漁ったらサンプルが合わない　うしさんのを取ってきたら通った
E を睨むと failure link 周りで処理が足りてなかったので書くと通る

2019-09-23

添え字 gcd での畳み込みで AGC038-C を解く

アルゴリズム

概要

数列 $a, b$ に対し、その gcd 畳み込み*1 $c$ を以下のように定めます

$\displaystyle c_{k} := \sum_{\gcd(i, j) = k} a_{i} * b_{j}$

この畳み込みの計算方法と、応用として AGC038-C を解く手順を示します。

数列の変換で要素毎の乗算に変換

数列から数列への関数 $f$ を導入します。 $f$ は倍数の和を取る変換で、以下のように定義します。

$\displaystyle f(x)_{k} := \sum_{k \mid i} x_{i}$

高速ゼータ変換を知っている人なら、上位集合に対する高速ゼータ変換の類似物と言えば分かりやすいかもしれません。実は、 $a$ と $b$ の gcd 畳み込み $c$ に対して以下の式が成立します。

$f(c)_{i} = f(a)_{i} * f(b)_{i}$

よって $f$ の逆変換 $f^{-1}$ が存在すれば

$c = f^{-1}(f(a) * f(b))$ ( $*$ は要素毎の乗算)

によって $c$ を計算することが出来ます。このような変換を考えることは、 $\gcd$ に限らず $|$ や $+$ による他種の畳み込みにも共通しています。

変換の概説

$\gcd(i, j) = k$ が成立しているとき、 $k \mid i$ かつ $k \mid j$ です。条件を緩めて、 $k \mid i$ かつ $k \mid j$ なる $i,j$ について $\sum a_{i} * b_{j}$ ならば容易に計算することが出来ます。これは $i,j$ が独立なので $f(a)_{k} * f(b)_{k}$ となるためです。
このような $i,j$ について $\gcd(i,j) = k$ は必ずしも成立しませんが、 $k \mid \gcd(i,j)$ が成立します。これをよく見ると、 $f(c)_{k}$ に等しいことが分かります。

$f,f^{-1}$ を $\Theta(N \log N)$ で計算する

・ $f$

各 $k$ について $k \mid i$ なる $i$ を全て足しても、調和級数により計算量は $\Theta(N \log N)$ となります。

void f(std::vector<int> &a) {
  int n = a.size();
  for (int k = 1; k < n; ++k) {
    for (int i = k * 2; i < n; i += k) {
      a[k] += a[i];
    }
  }
}

$f^{-1}$ を考える都合により、in-place な実装例を示しました。
値を重複して足し合わせないよう、 $k$ のループは昇順に回す必要があります。

・ $f^{-1}$

$f$ のアルゴリズムをそのまま逆順になぞればよいです。
同様に、 $k$ の昇降に気を付けてください。

void inv_f(std::vector<int> &a) {
  int n = a.size();
  for (int k = n - 1; k >= 1; --k) {
    for (int i = k * 2; i < n; i += k) {
      a[k] -= a[i];
    }
  }
}

$i$ のループは並列なのでどちらに回しても大丈夫です。

AGC038-C を解く

数列 $x$ を $x_{i} := i * (i が A に含まれる個数)$ として、 $x$ と $x$ を gcd 畳み込みしたものを $y$ とします。

$\displaystyle y_{k} = \sum_{\gcd(A_{i},A_{j}) = k} A_{i} * A_{j}$

であるため、基本的には

$\displaystyle \sum \frac{y_{k}}{k}$

を計算すればよいです。
ただし $i \lt j$ に反する $i=j,i \gt j$ が足されていることに注意して、補正します。

以上を実装した提出がこちらです。

Submission #7659496 - AtCoder Grand Contest 038

先頭に貼ってある modint 構造体については以下の記事を参照してください。

modint 構造体を使ってみませんか？ (C++) - noshi91のメモ

更なる高速化 $\Theta(N (\log(\log(N)))$

高速ゼータ変換の約数版 - noshi91のメモ

こちらのコードなどを参照してください。

*1:正式名称が分かっていません、ご存知の方がいらっしゃいましたら教えてください

概要

拡張ユークリッド互除法で逆元を求める

互除法のショートカット

から までの逆元を で計算する

個の数について全ての逆元を で計算する

概要

計算量解析

Queue の場合

Deque の場合

本記事は皆様からのご指摘を募集しております

誤った記述があるかもしれません

概要

最良計算量

最悪計算量

平均計算量

期待計算量

償却計算量

2021/12/19 追記

余談

概要

実装

概要

反例のようなもの (反例ではない)

証明

まとめ

概要

数列の変換で要素毎の乗算に変換

変換の概説

を で計算する

・

・

AGC038-C を解く

関連する文献

更なる高速化

$1$ から $N$ までの逆元を $\Theta ( N )$ で計算する

$N$ 個の数について全ての逆元を $\Theta ( N + \log M )$ で計算する

$\rm best :$ 最良計算量

$\rm worst :$ 最悪計算量

$\rm average :$ 平均計算量

$\rm expected :$ 期待計算量

$\rm amortized :$ 償却計算量

$f,f^{-1}$ を $\Theta(N \log N)$ で計算する

・ $f$

・ $f^{-1}$

更なる高速化 $\Theta(N (\log(\log(N)))$