ジョイジョイジョイ

ジョイジョイジョイジョイジョイ

符号付き整数を二冪で割る

符号付き整数を 2 冪で割るときに右算術シフトを使ってしまって痛い目にあったことのある皆さんこんにちは。僕はあなたの仲間です。

コンピュータで割り算を行うのは一般に重いので、できるだけシフト等の軽い操作で済ませたいですよね。

最近のコンパイラは賢いので適当に書けば最適化してくれますが、何かしらの事情で自分で最適化しなければならいないときのためにやり方を紹介しておきます。

C99 で 32 ビット符号付き整数 x1<<n で割ることを考えます。

符号なし整数なら二冪の割り算 x/(1<<n) と右シフト x>>n は一致するのですが、符号あり整数については、割り算は 0 方向への切り捨てなのに対して、右算術シフトは切り下げとなります。

例えば (-3) / 2 の結果が -1 なのに大して (-3) >> 1 の結果は -2 です。

これは困りました。

簡単な対策として、x が負のときシフトを行う前に (1<<n)-1 を足すという方法があります。

こうすると、負の数のとき切り上げになり、割り算と同じになります。

ですが、if は一般に重いので、できるだけ使わずに済ませたいです。

そこで、全体として以下の手順を踏みます

  1. x31 だけ右算術シフトしたものを y とする。
  2. y32-n だけ右論理シフトしたものを z とする。
  3. xz を足す。
  4. xn 右論理シフトする。

1, 2 が重要です。これによって zx が正のとき 0 に、負のとき (1<<n)-1 になっています。

たとえば n=6 のとき

x = 0xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y = 00000000000000000000000000000000
z = 00000000000000000000000000000000

x = 1xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
y = 11111111111111111111111111111111
z = 00000000000000000000000001111111

となっています。

正の数のときは最上位ビットが 0 なので、算術シフトによって 0 が引き出され全て 0 で埋まります。論理シフトをしても結果は変わりません。

対して負の数のときは最上位ビットが 1 なので、算術シフトによって 1 が引き出され全て 1 で埋まり、続いて論理シフトによって上位ビット 32-n ビットが全て 0 で埋まります。 1 は全部で n 個並んでいるのでこの値は (1<<n)-1 となります。

例えば以下のように実装できます。

gist.github.com

C99 では論理シフトになるか算術シフトになるかは処理系依存であることに注意してください。

上記のコードは gcc 5.4.0 で動作を確認しています。

おまけ

ちなみに最近のコンパイラでも、ここで紹介した感じになります。

以下は clang 3.8.0x / 64MIPSコンパイルした結果(の一部)です。

sra $2, $1, 31
srl $2, $2, 26
addu $1, $1, $2
sra $1, $1, 6

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] Henry S. Warren. (2012). Hacker’s Delight (2nd Edition)

シュトラッセンのアルゴリズムとその導出の仕方

 O(n^{3}) よりも小さい時間計算量  O(n^{\log_{2} 7}) \approx O(N^{2.807}) で行列積を計算するシュトラッセンアルゴリズムとその導出方法を紹介します。

シュトラッセンアルゴリズム

 n \times n の行列積  A B = C を考えます。

まず  A, B, C

$$ \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{pmatrix} $$

というように  \frac{n}{2} \times \frac{n}{2} 行列に分解します。

 n が奇数のときは余った行、列を適当にどちらかに分けてあげてください。

まず  7 つの行列積

$$ \begin{align} P_{1} &= ( A_{11} + A_{22} ) ( B_{11} + B_{22} ) \\ P_{2} &= ( A_{21} + A_{22} ) B_{11} \\ P_{3} &= A_{11} ( B_{12} - B_{22} ) \\ P_{4} &= A_{22} ( B_{21} - B_{11} ) \\ P_{5} &= ( A_{11} + A_{12} ) B_{22} \\ P_{6} &= ( A_{21} - A_{11} ) ( B_{11} + B_{12} ) \\ P_{7} &= ( A_{12} - A_{22} ) ( B_{21} + B_{22} ) \end{align} $$

を計算します。

すると、 C の各ブロックは

$$ \begin{align} C_{11} &= P_{1} + P_{4} - P_{5} + P_{7} \\ C_{12} &= P_{3} + P_{5} \\ C_{21} &= P_{2} + P_{4} \\ C_{22} &= P_{1} + P_{3} - P_{2} + P_{6} \end{align} $$

となります。

この証明は両辺を頑張って展開して一致することを確かめればできます。

ここで、各  P_{i} を計算するときに、再帰的にこのアルゴリズムを適用することを考えます。

各ステップの計算は再帰部分以外は行列の加減算なので、 O(n^{2}) でできます。

再帰は大きさが半分の問題を  7 回呼び出しています。

よって、 Master theorem より、全体の計算量は  O(n^{\log_{2} 7}) \approx O(N^{2.807}) です。

導出方法

 P_{i} C_{ij} の式が与えられれば証明は頑張って計算するだけですが、 P_{i} C_{ij} の式がどうやって出てきたのか分かりません。

ここでは簡単な(とはいっても骨の折れる)導出方法を一つ紹介します。

流れは以下のとおりです。

  1.  2 \times 2 の実行列の掛け算で、要素同士の掛け算を  7 回することを目指す
  2. 平面上の  120 度回転を表す行列  D をとる
  3. 非零ベクトル  u をとる
  4.  u と直交して  u^{\bot} D u = 1 を満たす横ベクトル  u^{\bot} をとる
  5.  X = u u^{\bot} とする
  6.  \{ D, X, D^{-1}XD, DXD^{-1} \} \{ D^{-1}, X, D^{-1}XD, DXD^{-1} \} 2 \times 2 の行列の基底になっている
  7.  A, B をそれぞれ上の基底で表現してかけると  7 回ですむ
  8.  D は実は固有多項式が同じなら同じ議論が成り立つので、簡単な行列をとってみる
  9.  u も簡単なベクトルをとってみる
  10. すると実は 7 はシュトラッセンアルゴリズムと同じだった

長いですね。

一つずつ見ていきましょう。

回転行列  D

 120 度 ( 2/3 \pi) の回転行列

$$ D = \begin{pmatrix} -\frac{1}{2} & -\frac{\sqrt{3}}{2} \\ \frac{\sqrt{3}}{2} & -\frac{1}{2} \end{pmatrix} $$

を取ります。

 tr(D) = -1,  det(D) = 1 なので、固有多項式 \lambda^{2} + \lambda + 1 です。

ケーリーハミルトンの定理より  D^{2} + D + E = 0 \cdots (1)が成り立ちます。

両辺に  D をかけると  D^{3} + D^{2} + D = 0 \cdots (2) が成り立ちます。

(1) と (2) の左辺どうしは等しいので  D^{3} = E \cdots (3) となります。これは、 D 120 度の回転行列であることを考えると当たり前のことです。

(3) より  D^{-1} = D^{2} \cdots (4) です。

よって、 (1) は  D^{-1} + D + E = 0 \cdots (5) と書き換えられます。

(5) より  tr(D^{-1}) = -tr(D) - tr(E) = 1 - 2 = -1 \cdots (6) となります。

非零ベクトル  u

非零ベクトル  u をとります。

 u と直交する非零の横ベクトル  u' をとります。

 D の固有多項式の実根は存在しないので、 u D固有ベクトルではありません。

よって、 Du u と別方向を向いているはずなので、  u'Du \neq 0 です。

 u^{\bot} = \frac{u'}{u'Du} とします。

 u^{\bot}u と直交し  u^{\bot} D u = 1 となります。

行列  X

 X = u u^{\bot} をとります。ここで、

 X^{2} = (u u^{\bot}) (u u^{\bot}) = u (u^{\bot} u) u^{\bot} = 0 \cdots(7)

 XDX = (u u^{\bot}) D (u u^{\bot}) = u (u^{\bot} D u) u^{\bot} = u u^{\bot} = X \cdots(8)

 XD^{-1}X = (u u^{\bot}) D^{-1} (u u^{\bot}) = u (u^{\bot} D^{-1} u) u^{\bot} = -u u^{\bot} = -X \cdots(9)

が成り立ちます。

一次独立性

前々節でも言ったように、 D の固有多項式の実根は存在しないので、 u D固有ベクトルではありません。

よって、 u Du は一次独立です。

同様に  u^{\bot} u^{\bot} D は一次独立です。

よって、これらの積  u \cdot u^{\bot}, u \cdot u^{\bot} D, Du \cdot u^{\bot}, Du \cdot u^{\bot} D 2 \times 2 行列の基底になります。

これらに両側から  D を掛けます。

 D^{-1} = D^{2} \cdots (4) を考えると、  DXD, DXD^{-1}, D^{-1}XD, D^{-1}XD^{-1} となります。

 D は正則なので、これらも  2 \times 2 行列の基底になります。

これらを足し合わせると  (D + D^{-1}) X (D + D^{-1}) = (-E) X (-E) = X となります。ひとつ目の等号は  D^{-1} + D + E = 0 \cdots (5) から従います。

よって、基底のうち一つを  X で置き換えても基底になります。

これより、  X, DXD^{-1}, D^{-1}XD は一次独立です。

基底表現

 X, DXD^{-1}, D^{-1}XD のトレースはそれぞれ、

 tr(X) = tr(u^{\bot}u) = tr(uu^{\bot}) = 0

 tr(DXD^{-1}) = tr(XD^{-1}D) = tr(X) = 0

 tr(D^{-1}XD) = tr(XDD^{-1}) = tr(X) = 0

となります。

よって、 \{ X, DXD^{-1}, D^{-1}XD \} から生成される  2 \times 2 行列の部分空間の要素のトレースは  0 です。

 tr(D) = -1, tr(D^{-1}) = -1 より、これらはそれぞれ  X, DXD^{-1}, D^{-1}XD と一次独立です。

よって、 \{ D, X, D^{-1}XD, DXD^{-1} \} \{ D^{-1}, X, D^{-1}XD, DXD^{-1} \} はそれぞれ  2 \times 2 行列の基底になっています。

$$ A = a_{1} D + a_{2} X + a_{3} D^{-1}XD + a_{4} DXD^{-1} $$

$$ B = b_{1} D^{-1} + b_{2} X + b_{3} D^{-1}XD + b_{4} DXD^{-1} $$

と表します。このとき基底同士の掛け算は

f:id:joisino:20170903160923p:plain

となります。ここでは (7), (8), (9) の結果をたくさん使っています。(表は [2] より引用)

同じ結果が同列、同行に現れているので、これらはまとめられます。

つまり、

$$ \begin{align} AB &=& a_{1} &\times& b_{1} &\times& id \\ &+& a_{2} &\times& ( b_{1} + b_{4}) &\times& XD^{-1} \\ &+& a_{3} &\times& ( b_{1} + b_{2}) &\times& D^{-1}X \\ &+& a_{4} &\times& ( b_{1} + b_{3}) &\times& DXD \\ &+& ( a_{1} - a_{4}) &\times& b_{2} &\times& DX \\ &+& ( a_{1} - a_{2}) &\times& b_{3} &\times& XD \\ &+& ( a_{1} - a_{3}) &\times& b_{4} &\times& D^{-1}XD^{-1} \cdots (10) \end{align} $$

となります。

もう少し簡単に

ここまで  D の性質としては  tr(D) = -1,  det(D) = 1 であることしか使っていません。よって、たとえば

$$ D = \begin{pmatrix} 0 & -1 \\ 1 & -1 \end{pmatrix} $$

ととっても構いません。このとき

$$ D^{-1} = \begin{pmatrix} -1 & 1 \\ -1 & 0 \end{pmatrix} $$

です。

また、 u にも任意性がありましたが、例えば  u = ( 1, 0 )^{{\rm T}} ととります。

このとき

$$ X = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} $$

となります。

このようにとることで上記の議論を  2 \times 2 行列の要素がブロック行列でも適用できることになります。( 1単位行列で、 0 は零行列です)

また、計算すると

$$ D^{-1}XD = \begin{pmatrix} -1 & 1 \\ -1 & 1 \end{pmatrix} $$

$$ DXD^{-1} = \begin{pmatrix} 0 & 0 \\ -1 & 0 \end{pmatrix} $$

となります。 A, B の要素ごとに考えると、

 X, DXD^{-1}, D^{-1}XD はトレースが  0 なので、  a_{1} = - tr(A) = - a_{11} - a_{22} となります。

また、 X, D^{-1}XD (2,2) 要素が  0 なので、 a_{3} = a_{22} + a_{1} = - a_{11} となります。

残りは、それぞれ異なる場所に  1, -1 が立っているだけなので、それぞれ

 a_{2} = a_{12} + a_{1} - a_{3} = a_{12} - a_{22}

 a_{4} = -( a_{21} - a_{1} + a_{3}) = - a_{21} - a_{22}

となります。

同様に

 b_{1} = -b_{11} - b_{22}

 b_{2} = b_{11} + b_{12}

 b_{3} = b_{22}

 b_{4} = b_{11} - b_{21}

となります。

また、瑣末な計算なのでここでは省きますが、 id, XD^{-1}, D^{-1}X, DXD, DX, XD, D^{-1}XD^{-1} もそれぞれ  1, 0, -1 からなる簡単な  2 \times 2 行列になっています。

最後に (10) で表したようにそれぞれ係数を掛けて足し合わせると、シュトラッセンアルゴリズムが出てくるというわけです。

最後に

何か間違いや質問などがあればお願いします。

もっと簡単な導出があれば教えていただけると嬉しいです。

最後まで読んでいただきありがとうございました。

参考文献

[1] Volker Strassen. (1969) Gaussian elimination is not optimal.

[2] Christian Ikenmeyer. Vladimir Lysikov. (2017) Strassen’s 2x2 matrix multiplication algorithm: A conceptual perspective.

[3] Master theorem - Wikipedia

ケイリーの公式の証明6種類

ケイリーの公式の証明たちの紹介です。

ケイリーの公式とは

ケイリーの公式とは  n 頂点のラベル付きの木の総数  T(n) T(n) = n^{n-2} であるという公式のことです。

ここで、ラベル付きであるとは、それぞれの頂点を区別するということです。

たとえば  n = 3 のとき、頂点を区別しない場合は長さ  3 のパスのみの  1 通りですが、ラベル付きの木の場合は  1 - 2 - 3,  2 - 1 - 3,  1 - 3 - 2 3 (= 3^{3-2} = T(3)) 通りです。

証明 1 (プリューファーコード) [1]

おそらく一番有名な証明です。

 n 頂点のラベル付きの木の集合から  \{ 1, 2, \ldots, n \}^{n-2} への全単射を以下のように構成します。

最もラベルが小さい葉を木から取り除き、その葉と繋がっていた頂点のラベル  a_{1} を数列の最初の値とします。

続けて、最もラベルが小さい葉を木から取り除き、その葉と繋がっていた頂点のラベル  a_{2} を数列の  2 番目の値とします。

以下同様に頂点が  2 つになるまで操作を続けます。

こうしてできた数列  a_{1}, a_{2}, \ldots, a_{n-2} が木の値となります。

この数列をプリューファーコードといいます。

例えば、以下の木のプリューファーコードは  3, 2, 6, 6 となります。

f:id:joisino:20170819151034p:plain

これの逆写像は以下のように与えられます。

まず、辺のない  n 頂点を用意します。

次に、  a_{1}, a_{2}, \ldots, a_{n-2} に含まれない値のうち、最も小さいものを  b_{1} とし、  a_{1} b_{1} を繋ぎます。

続けて、  b_{1}, a_{2}, a_{3} \ldots, a_{n-2} に含まれない値のうち、最も小さいものを  b_{2} とし、  a_{2} b_{2} を繋ぎます。

続けて、  b_{1}, b_{2}, a_{3}, a_{4} \ldots, a_{n-2} に含まれない値のうち、最も小さいものを  b_{3} とし、  a_{3} b_{3} を繋ぎます。

以下同様に木を構成していきます。

最後に、 b_{1}, b_{2}, \ldots, b_{n-2} に含まれていない  2 頂点を繋いで完成です。

これが逆写像になっていることは以下のようにわかります。

まず、プリューファーコードに登場する頂点は元の木の内点で、かつ元の木の内点を全て尽くしています。

なぜなら

  •  a_{i} i 番目の操作の時点で葉とつながっている点なので、これが葉であると、操作後に頂点数が 1 の成分になってしまうので矛盾です。よって全ての  a_{i} は内点です。

  • 全ての内点は、操作のどこかの時点で葉になり、 2 頂点以外は取り除かれていきます。少なくともその頂点が葉になった瞬間においてコード中に登場するはずなので、全ての内点を尽くしています。

以上より、逆に  a_{1}, a_{2}, \ldots, a_{n-2} に含まれない値は初期状態において葉であって、全ての葉を尽くしています。

そのうち最もラベルの小さい頂点は順写像において最初に取り除かれた頂点で、 a_{1} と繋がっていたはずです。

以下帰納的に順写像の操作に対応する辺が見つかっていきます。

最後は、順写像の操作で取り除かれなかった頂点を結べば完成となります。

全単射が存在するということは要素の数が等しいということです。

 \{ 1, 2, \ldots, n \}^{n-2} の要素数はもちろん  n^{n-2} なので  T(n) = n^{n-2} となります。

証明 2 (もうひとつの全単射) [2]

もうひとつ全単射を使った証明を紹介します。

今度は  \{ 1, 2, \ldots, n \}^{n} から、  n 頂点からなるラベル付きの木であって  1 頂点に赤い印がついていて  1 頂点が星形である木への全単射を構成します。

ここで、赤の頂点と星の頂点は同じ頂点でもかまいません。

ラベル付きの木を固定すると赤の頂点と星の頂点の選び方は  n^{2} 通りあるので、値域の要素数 T(n) n^{2} で、定義域の要素数 n^{n} なので、これらが等しければケイリーの公式の証明が完了します。

数列  (a_{1}, a_{2}, \ldots, a_{n}) を入力とします。

 1 \le i \le n について  i から  a_{i} へ辺があるような  n 頂点  n 辺からなるグラフを考えます。

全ての出次数が  1 であるグラフは functional graph と呼ばれ、サイクルに木がついた連結成分がいくつかあるような形になります。

例えば  7, 4, 8, 1, 7, 4, 4, 3 のときグラフは以下のようになります。

f:id:joisino:20170819161216p:plain

ここで、サイクル上の頂点をソートして並べたものを  b_{1}, b_{2}, \ldots, b_{k} とします。

上の例の場合  1, 3, 4, 7, 8 となります。

出力の木としては、  a_{b_{1}}, a_{b_{2}}, \ldots, a_{b_{k}} をこの順に並べてパスとし、残りは functional graph と同様に繋げます。赤の頂点は  a_{b_{1}} を選び、星の頂点は  a_{b_{k}} を選びます。

上の例の場合は以下のようになります。

f:id:joisino:20170819162001p:plain

写像は以下のように構成します。

まず、木のうち、赤の頂点から星の頂点までの(一通りしかない)パスをとって、 c_{1}, c_{2}, \ldots, c_{k} と順に並べます。

パス上の頂点のラベルをソートした数列を  b_{1}, b_{2}, \ldots, b_{k}とします。

 1 \le i \le k について、元の数列の  b_{i} 番目の値が  c_{i} であることが分かります。

パス以外の頂点は、パス上の頂点を根とした根付き木を考え、頂点  v の親が  u のとき、元の数列の  v 番目の値が  u となることが分かります。

これで復元できました。

証明 3 (double counting) [2] [3]

この証明もかなり有名な証明です。

wikipedia にもこの証明が載っています。

頂点にも辺にもラベルが付いた  n 頂点からなる根付き木の総数  U(n) を二通りの方法で数えます。

一つは、 T(n) を使う方法です。ラベル付きの木を固定したとき、根の選び方が  n 通り、辺へのラベルの付け方が  (n-1)! 通りあるので、 U(n) = T(n) \cdot n \cdot (n-1)! となります。

もうひとつは辺を  1 つずつ追加していって逐次的に木を構成していく方法です。 i 番目に追加した辺のラベルを  i とします。

辺のない  n 頂点からなるグラフを用意します。

最初の辺は、親の選び方が  n 通りあって、子の選び方は自分以外の  n-1 通りです。

 2 つ目の辺は、親の選び方がこれも  n 通りあります。子の選び方は、根でない頂点を選ぶとその頂点の親が  2 つになってしまうのでだめです。自分が属している木の根も、ループができてしまうのでだめです。よって、子の候補は  n-2 個となります。

同様に、 1 \le i \le n-1 について  i 番目の辺は、親の選び方が  n 通り、子の選び方が  n-i 通りとなります。

よって、 U(n) = \Pi_{i=1}^{n-1} n (n-i) = n^{n-1} \cdot (n-1)! となります。

これが  T(n) \cdot n \cdot (n-1)! と等しいので  T(n) = n^{n-2} となります。

証明 4 (行列木定理) [2]

行列木定理 を使った証明です。

ラベル付きの木の総数と  n 頂点完全グラフ  K_{n} の全域木の総数は同じです。

 K_{n}ラプラシアン行列は

$$ \displaystyle \left[ \begin{array}{rrr} n-1 & -1 & -1 & \ldots & -1 \\ -1 & n-1 & -1 & \ldots & -1 \\ -1 & -1 & n-1 & \ldots & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1 & -1 & -1 & \ldots & n-1 \end{array} \right] $$

となります。行列木定理より、 T(n) はこの行列の任意の余因子と一致します。

 (1, 1) 余因子をとりましょう。

 (1, 1) 余因子はもとの行列から  1 行目と  1 列目を取り除いた  (n-1) \times (n-1) 行列の行列式です。

まず、 2 行目から  n-1 行目を  1 行目に足します。すると  1 行目が全て  1 になります。

$$ \displaystyle \left| \begin{array}{rrr} 1 & 1 & 1 & \ldots & 1 \\ -1 & n-1 & -1 & \ldots & -1 \\ -1 & -1 & n-1 & \ldots & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -1 & -1 & -1 & \ldots & n-1 \end{array} \right| $$

次に、 1 行目を  2 行目から  n-1 行目に足します。

$$ \displaystyle \left| \begin{array}{rrr} 1 & 1 & 1 & \ldots & 1 \\ 0 & n & 0 & \ldots & 0 \\ 0 & 0 & n & \ldots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & n \end{array} \right| $$

この式の値は  n^{n-2} なので、  T(n) = n^{n-2} となります。

証明 5 (漸化式) [2]

より一般的な問題を漸化式と帰納法によって解きます。

ラベルの集合の、要素数  k の部分集合  A を任意にとります。

 A に含まれている頂点が互いに異なる連結成分に属するようなラベル付き森の数を  S(n,k) とします。

これは  A の要素数にしかよらないので well-defined です。

便宜上  S(0,0) = 1, S(n,0) = 0 ( n > 0 ) とします。

 S(n,k) = k \cdot n^{n-k-1} (*) であることを示します。

 T(n) = S(n,1) なので、これが示せればケイリーの公式は示されます。

ラベルの集合を  \{ 1, 2, \ldots, n \} とし、 A = \{ 1, 2, \ldots, k \} としたときに  S(n, k) で数える森を  1 つ固定します。

この森で頂点  1 と繋がっている頂点を  B = \{ a_{1}, a_{2}, \ldots, a_{l} \} とします。

このときこの森から頂点  1 を取り除いてできる森は、ラベルの集合を  \{ 2, 3, \ldots, n \} とし  A = \{ 2, 3, \ldots, k, a_{1}, a_{2}, \ldots, a_{l} \} としたときに  S(n-1,k+l-1) で数える要素になります。

逆に ラベルの集合を  \{ 2, 3, \ldots, n \} とし  A = \{ 2, 3, \ldots, k, a_{1}, a_{2}, \ldots, a_{l} \} としたとき、新たに頂点  1 を加えて  a_{1}, a_{2}, \ldots, a_{l} とつなげると、先ほど考えていた  S(n, k) で数えた森に対応します。

つまり、 S(n,k) で数える森は、 A と交わらない頂点集合  B ごとに  S(n-1,k+|B|-1) 個の森と対応づいているということです。

以上の考察から

 \displaystyle S(n,k) = \sum_{l=0}^{n-k} \binom{n-k}{l} S(n-1, k+l-1)

となります。これは  B として、 A 以外の要素 n-k 個から  l 個選ぶというのを  l = 0, 1, \cdots n-k について足しあわせているということです。

この漸化式を使って帰納法により (*) を示しましょう。

 
\begin{align*}
S(n,k) &= \sum_{l=0}^{n-k} \binom{n-k}{l} S(n-1, k+l-1) & \\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (k+l-1) (n-1)^{n-k-l-1} & ({\rm 帰納法の仮定より})\\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (n-1-l) (n-1)^{l-1} & (l = n - k - l {\rm と変換}) \\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (n-1) \cdot (n-1)^{l-1} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & ({\rm 分配}) \\
&= \sum_{l=0}^{n-k} \binom{n-k}{l} (n-1)^{l} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & \\
&= (n-1+1)^{n-k} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & ({\rm 二項定理})\\
&= n^{n-k} - \sum_{l=0}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & \\
&= n^{n-k} - \sum_{l=1}^{n-k} \binom{n-k}{l} l \cdot (n-1)^{l-1} & (l=0 {\rm のとき} 0) \\
&= n^{n-k} - (n-k) \sum_{l=1}^{n-k} \binom{n-k-1}{l-1} \cdot (n-1)^{l-1} & ({\rm 二項係数をほぐす})\\
&= n^{n-k} - (n-k) \sum_{l=0}^{n-k-1} \binom{n-k-1}{l} \cdot (n-1)^{l} & (l=l+1 {\rm と変換})\\
&= n^{n-k} - (n-k) (n-1+1)^{n-k-1} & ({\rm 二項定理})\\
&= n^{n-k} - (n-k) n^{n-k-1} & \\
&= k \cdot n^{n-k-1} &
\end{align*}

証明 6 (母関数とラグランジュの反転公式) [1]

 n 頂点からなるラベル付きの根付き木の総数を  R(n) とします。便宜上  R(0) = 0 とします。

根なし木を固定したとき根の選び方は  n 通りなので、 R(n) = n T(n) となります。 R(n) = n^{n-1} を示すのが目標です。

根の選び方は  n 通りあります。

根の子の個数を  m に固定します。

子に順番があると考え、それぞれの部分木の大きさが左から順に  k_{1}, k_{2}, \ldots, k_{m} (k_{1} + k_{2} + \ldots + k_{m} = n-1) であるとします。

根の子の部分木について、どのラベルを割り当てるかが  \frac{(n-1)!}{k_{1}! k_{2}! \ldots k_{m}!} 通りあります。

各部分木について根付き木の作り方はそれぞれ  R(k_{1}), R(k_{2}), \ldots, R(k_{m}) 通りあります。

子の順番を取り除くために、これらの積を  m! で割る必要があります。

以上より

 
\begin{align*}
R(n)
&= n \sum_{m=0} \frac{1}{m!} \sum_{k_{1} + k_{2} + \ldots + k_{m} = n-1} \frac{(n-1)!}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) \\
&= \sum_{m=0} \frac{1}{m!} \sum_{k_{1} + k_{2} + \ldots + k_{m} = n-1} \frac{n!}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) \\
\end{align*}

という漸化式が成り立つことが分かります。

次に  R(n) の指数型母関数

 \displaystyle f(x) = \sum_{n=0}^{\infty} \frac{ R(n) x^n }{ n! }

を考えます。

 
\begin{align*}
x e^{f(x)}
&= x \sum_{m=0}^{\infty} \frac{ f(x)^m }{ m! } \\
&= x \sum_{m=0}^{\infty} \frac{1}{m!} (\sum_{n=0}^{\infty} \frac{ R(n) x^n }{ n! })^m \\
&= x \sum_{m=0}^{\infty} \frac{1}{m!} \sum_{k_{1}=0}^{\infty} \sum_{k_{2}=0}^{\infty} \cdots \sum_{k_{m}=0}^{\infty} \frac{1}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) x^{k_{1} + k_{2} + \ldots + k_{m}} \\
&= \sum_{m=0}^{\infty} \frac{1}{m!} \sum_{k_{1}=0}^{\infty} \sum_{k_{2}=0}^{\infty} \cdots \sum_{k_{m}=0}^{\infty} \frac{1}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) x^{k_{1} + k_{2} + \ldots + k_{m}+1} \\
&= \sum_{m=0}^{\infty} \frac{1}{m!} \sum_{k_{1}=0}^{\infty} \sum_{k_{2}=0}^{\infty} \cdots \sum_{k_{m}=0}^{\infty} \frac{(k_{1} + k_{2} + \ldots + k_{m} + 1)!}{k_{1}! k_{2}! \ldots k_{m}!} R(k_{1}) \cdot R(k_{2}) \cdot \ldots \cdot R(k_{m}) \frac{x^{k_{1} + k_{2} + \ldots + k_{m}+1}}{(k_{1} + k_{2} + \ldots + k_{m} + 1)!} \\
&= f(x)
\end{align*}

最後の等式は先ほどの漸化式から従います。

よって、 R(n) の指数型母関数  f(x) f(x) = x e^{f(x)} という単純な方程式に従うことがわかりました。

次にラグランジュの反転公式を使います。

ラグランジュの反転公式とは、べき級数

 \displaystyle f(x) = a + \frac{a_{1}}{1!} x + \frac{a_{2}}{2!} x^{2} + \ldots

 \displaystyle f(x) = a + x \phi( f(x) )

という方程式を満たすとき、

 \displaystyle a_{n} = \left| \frac{{\rm d}^{n-1}}{{\rm dt}^{n-1}} (\phi(t)^n) \right|_{t=a}

を満たすというものです。

これに先ほどの方程式を当てはめます。

(本当は収束性などを調べないといけないはずですが、参考文献では特に言及されていませんでした。実際  R(n) = n^{n-1} なので収束半径は  \frac{1}{e} になるはずなのですが、答えを使わずに証明する良い方法が思いつかなかったのでここでは省略します。僕が勘違いしていたり、良い方法があったりしたら教えてください)

ここで、 \phi(t) = e^{t}, a = 0 なので、

 \displaystyle \frac{{\rm d}^{n-1}}{{\rm dt}^{n-1}} e^{nt} = n^{n-1} e^{nt}

よって、 t = a = 0 を代入すると  R(n) = n^{n-1} となり示されました。

最後に

何か間違いや質問などがあればお願いします。

他の証明もあれば教えていただけると嬉しいです。

最後まで読んでいただきありがとうございました。

参考文献

[1] パラモノヴァ(著). ヴィンベルク(著). コーハシ(著). 武部尚志(訳). (2007). モスクワの数学ひろば 3 代数篇/対称性・数え上げ.

[2] Aigner Martin. Ziegler Günter M. (2013). Proofs from THE BOOK 4th edition. pp.201-206

[3] ケイリーの公式 - Wikipedia

前処理O(n)クエリO(1)のLCAと静的RMQ

時間計算量 <O(n), O(1)> の LCA(Lowest Common Ancestor) と RMQ(Range Minimum Query) を C++ で実装しました。

アルゴリズムの解説はDさんのスライド [1] LCA and RMQ ~簡潔もあるよ!~ がとても分かりやすいのでそちらを参照してください。

概要だけ説明します。

LCA の概要

LCA は頂点を dfs 順で訪れた順に並べると深さの列の RMQ に帰着されます。このことは [2] 蟻本 などに載っています。

この列は隣り合う数の差がちょうど  1 になっています。

この列を  \frac{\log n}{2} 個ずつのブロックに分け、それぞれのブロック内の最小値を求めます。

ブロックの数は  \frac{2n}{\log n} 個になるので、ブロックの区間の最小値を求めるクエリは sparse table を使うと前処理  \frac{2n}{\log n} \log \frac{2n}{\log n} = O(n)、クエリ  O(1) で処理できます。

ブロックの中についてですが、各ブロックは 1 つ目の数がわかればあとは +1 されるか -1 されるかの 2 通りです。この遷移の種類は高々  2^{\frac{\log n}{2} - 1} = O(\sqrt{n}) 個しかありません。

この全ての種類について、最初の数が  0 だと思って愚直に前計算すると全て合わせて  O(\log^{2} n \cdot \sqrt{n}) で計算できます。

あとはクエリが来た時に、完全に含まれるブロックを sparse table で処理して、はみ出た部分を前計算した表を参照して答えればよいです。

RMQ の概要

列を Cartesian Tree という木に変換します。

Cartesian Tree は、列のインデックスをノードの値にもつ二分探索木で、列の値についてヒープ条件を満たしています。ヒープ条件を満たしているとは、列の値について、親の値が子の値以下であるということです。

Treap のような感じです。

すると、RMQ はこの Cartesian Tree 上の LCA クエリに帰着されることがわかると思います。

Cartesian Tree の構成の仕方は [1] には載っていないので少し詳しく説明します。

Cartesian Tree を O(n) で構成する

列を左から右に見ていって、順に Cartesian Tree を構成していきます。

今まで構成した Cartesian Tree の頂点の親を配列で管理し、根から右の子を辿るパスをスタックで管理します。

新しく列の最後に値を追加するとします。

このとき、二分探索木の条件から、新しい頂点はスタックで管理されているパス上の頂点の右の子になるか、根になるかのどちらかです。

どの頂点の子になるかは、スタックの頂点の値よりも新しい頂点の値の方が小さければ、スタックから頂点を取り除くことを繰り返し、最終的にスタックの一番上に残っている頂点が新しい頂点の親になります。スタックが空ならば新しい頂点は根になります。

新しい頂点のせいで親から引き離された頂点は、新しい頂点の左の子にすればよいです。

例えば 4, 2, 5, 1, 3, 8, 7, 9 の後ろに 6 を付け加える場合は以下のようになります。

ノードは (index, val) のように表しています。赤いノードがスタックで管理されているパスです。

f:id:joisino:20170813142841p:plain

f:id:joisino:20170813142849p:plain

実装

C++ で実装しました。

LCA

gist.github.com

RMQ

gist.github.com

LCA は GRL_5_C に投げたところ、 10^{5} 頂点  10^{5} クエリで AOJ 上で 0.06s でした。入出力にもある程度時間がかかっていると思います。

RMQ の方は、DSL_3_D (スライド最小値) に投げたところ、  10^{5} 要素  10^{5} クエリで AOJ 上で 0.34s でした。 RMQ は何段階も帰着させている分線形とはいえ遅いです。

メルセンヌツイスタで生成した乱数に対して RMQ を実行すると、手元の環境だと

  •  10^{6} 要素  10^{6} クエリで 0.425s
  •  10^{7} 要素  10^{7} クエリで 5.206s

でした。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] LCA and RMQ ~簡潔もあるよ!~

[2] 秋葉拓哉. 岩田陽一. 北川宜稔. (2012). プログラミングコンテストチャレンジブック [第2版]. pp.294 - 295.

テヘラン旅行

IOI2017(イラン大会)に副団長として参加してきました。

ioi2017.org

IOIの参加記は選手・役員の感想(まだ公開されていないようなのであとでリンクを張ります)に書いたのでそちらを参照してください。

第29回国際情報オリンピック イラン大会 速報 も面白いのでまだ見てないひとはぜひ。

この記事では、これからテヘランに行く人の参考のために、大会以外の部分についてメモしておきます。

事前準備

買い物をして事前準備をしました。

地球の歩き方

E06 地球の歩き方 イラン ペルシアの旅2017~2018

念の為買っておきましたが、ツアー旅行のようなものなのでほとんど使いませんでした。

イランの一般的な情報と、ペルシャ語(特に数字)については少し役に立ったと思います。

圧縮袋

【日本製】旅行 出張 に便利な 衣類 圧縮袋 リムーブエアー ●特許製法 逆止弁 使用● M・L 各5枚 10枚組

服を圧縮する袋です。

かなり容量が小さくなるので助かりました。これはテヘラン旅行以外でも役に立つと思います。

盗難防止

首下げバッグ 貴重品ケース iPhone 6/6S/7 Plus収納可 パスポート入れ スキミング防止 バッグインバッグ 斜め掛け アウトドア 海外旅行 トラベル ネックポーチ

これのウエストポーチを買いました。

そんなに大きくはなく、パスポートと現金と鍵とカードなどを入れていました。

テヘランはそこまで治安が悪くなかったので、そこまで厳重に警戒する必要は無かったかもしれません。(スマホなどはポケットに入れていましたが取られませんでした)

バザールに行くなどすると危なかったのかもしれません。

貴重品を一箇所に集められるという点ではかなり重宝しました。

マルチビタミン

DHC マルチビタミン (60日分) 60粒

イランの食事はインディカ米、肉、インディカ米、肉、インディカ米という感じなのでビタミンが不足しがちです。

これのおかげ(?)で風邪を引きませんでした。

サイリウム

キングブレードX10III Neo シャイニング

帰国したあとすぐに映画を見るために持って行きました。

テヘランで使うつもりは無かったのですが、懐中電灯の代わりに使えて便利でした。

二人部屋で夜中起きたときにトイレの電気ガチャ(外れると部屋全体の電気が付く)をせずに済みました。

現地の気候

7 月末から 8 月頭と夏真っ盛りに行ってきました。

気温は最高気温が 35 度ちょっとの日が多く、最高気温 40 度の日もありましたが、湿度がかなり低いので、そこそこ過ごしやすかったです。個人的には夏の京都よりも過ごしやすいと思いました。

テヘランは砂漠なので木陰があまり無かったですが、建物の影や数少ない木陰に隠れると日中でも割と快適でした。

一方、砂が舞っているのと、自動車の排気ガスと、タバコのせいで空気が悪く、乾燥との二重苦で鼻が痛くなってきます。

これを防ぐ手立てはよくわからないのですが、大気汚染が特に酷そうな場所を歩くときは布で鼻を押さえるなどの対策をした方がよいと思います。

また、乾燥のせいで静電気もすごかったので心構えしておきましょう。

現地の買い物

まず一番気をつけないといけないのが、イランではリアルとトマーンという二種類の通貨が使われていることです。

日常的にはどちらかというとトマーンが多く使われていて、外国人向けの店ではリアルで書いてあるものが多い気がしましたが、どちらのパターンもあって、しかも通貨が書いてない場合が多かったので、よく気をつける必要があります。

僕が行った当時は 300リアル = 30 トマーン  \approx 1 円くらいでした。

また、もう一つ気をつけないといけないのが、多くの店では、値札がペルシャ数字で書かれていることです。

イランに行く前にペルシャ数字を覚えるか、少なくともペルシャ数字表を携帯していくことをオススメします。

物価は日本よりかなり安く、1 / 3 から 1 / 2 くらいな印象でした。

500mL ミネラルウォーターが 10000 リアル  \approx 30円くらいでした。

ホテル

Azadi Hotel というホテルに泊まりました。

一泊二万円以上するホテルだけあって、かなり綺麗で快適でした。

食事とトイレ以外は他の国の一級ホテルとほとんど変わらないと思います。

食事

前述しましたが、イランの食事はインディカ米、肉、インディカ米、肉、インディカ米という感じでした。

f:id:joisino:20170806163526j:plain

こういう感じです。

ホテルの朝食ではパンが出たので、朝食はパンをよく食べました。

f:id:joisino:20170806163637j:plain

こういう感じです。

風習など

まず、イスラム教の色々な風習がありますが、ここに書くには重いので他の文献に任せます。

一番驚いたのは、イラン人はかなり気さくで、街中を歩いているとよく話しかけてきたことです。

特に IOI 期間中は日本国旗が描かれた名札を下げていたこともあって、"Japon?“ と話しかけてくるケースが多かったです。その後の会話は様々ですが、適当に日本っぽい単語を言ってくることが多かったです。ポテトチップス的なお菓子を勧められたりもしました(怖かったので断りました)

あと気をつけたほうが良いのが、トイレ事情です。

イランのトイレはアナログウォシュレット(シャワーがトイレの横についている)を使って洗って、紙ナプキンのようなトイレットペーパーで拭いてゴミ箱に捨てる形式です。

外国人が多く泊まりそうな一級ホテルの Azadi Hotel でもこの形式だったのでイランのトイレはだいたいどこもこうなんだと思います。

IOI 期間中は外国人が集まっていたこともあり、ウォシュレットのせいで公共トイレの床はどこも大洪水でした。

気になる人はトイレに流せるトイレットペーパーを持って行くとよいと思います。

インターネット

イランはインターネット規制が厳しい国で、twitter などは規制されていてアクセスできません。

なので、規制されているサービスを利用したい場合は VPN を用意する必要があります。

僕は大学の VPN (KUINS-PPTP) を利用しました。

また、irancell という会社の SIM カードが IOI から配られました。

irancell の店はテヘランの至るところにあり、テヘラン旅行の際はこれを使うことになることが多そうです。

SIM カードをセットすると SMS でいろいろ情報が送られてきますが、全てペルシャ語なので USSDコード一覧 の Switching Language into English を使うことをオススメします。

回線は容量に応じて課金しました。

RechargeOption の Online Recharge からオンラインで課金できます。

課金したあとは Packages からパッケージを買います。

3GB の通信で 13000 トマーン  \approx 400 円 と激安です。

回線もそこそこ安定していました。

ホテルの wifi もあったのですが、こちらは通信が絶望的に不安定で使い物になりませんでした。具体的には 8.8.8.8 への ping が半分くらいパケットロスしました。

終わりに

テヘラン旅行はとても楽しく、また、IOI がたくさんお金をかけていたこともあって、想像以上に快適でした。

また行きたいかと言われると返答に困りますが。

もしテヘラン旅行に行くなら、楽しんできてください。

最後まで読んでいただきありがとうございました。

ミラーラビン素数判定法とロー法

ミラーラビン素数判定法について

ミラーラビン素数判定法はフェルマーテストを拡張した感じです。

フェルマーテストでは  p素数のとき、フェルマーの小定理より  a p が互いに素なとき  a^{p-1} \equiv 1 \pmod p となることを利用して、  n と互いに素な整数  a を取って  a^{n-1} \not \equiv 1 \pmod n であれば合成数と判断するのでした。

ここから更に拡張します。

 x^{2} \equiv 1 \pmod p のとき、 (x-1)(x+1) \equiv 0 \pmod p で、 p素数より  x \equiv 1 または  x \equiv -1 となることを利用します。

 (n-1) = 2^{q} d というように  d に素因数  2 が出ないように分解します。

このとき、( q \ge 1 であれば) a^{2^{q-1}d} = 1 または  a^{2^{q-1}d} = -1 となります。

さらに、 a^{2^{q-1}d} = 1 のときは( q \ge 2 であれば)  a^{2^{q-2}d} = 1 または  a^{2^{q-2}d} = -1 と続きます。

つまり、 p素数ならば、 a^{2^{q-i}d} i = 0, 1, 2, ... q のとき  1, 1, ..., 1, -1, ?, ?, ... となるか、全部  1 になるます。

これの対偶を使って、 a^{d} = 1 となるか、 a^{2^{q-i}d} (1 \le i \le q-1) のどこかで  -1 となるかのどちらかにならなければ、合成数です。

この  a をランダムに決めて、テストしてくいくのがミラーラビン素数判定法です。

ここで、テストする数  n が小さければいくつかの  a を使うだけで、確実に素数判定を実行できることが保証できます。

例えば、 n \le 4759123140 (\ge 2^{32}) なら、 a = 2, 7, 61 についてのみ調べれば良いです。

実装

C++ で実装しました。

gist.github.com

メルセンヌツイスタでランダムに生成した  10^{6} 個の  10^{9} 以下の整数について素数判定するのに手元の環境だと 1.854s かかりました。

ロー法について

乱択で素因数を発見するアルゴリズムです。

[2] 素因数分解の高速なアルゴリズム(ロー法) | 高校数学の美しい物語 がとても分かりやすいので、アルゴリズムの内容はそちらを見てください。

実装

これの素数判定にミラーラビンを使って、C++ で実装しました。

gist.github.com

メルセンヌツイスタでランダムに生成した  10^{5} 個の  10^{9} 以下の整数について素数判定するのに手元の環境だと 1.670s かかりました。

最後に

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] ミラー–ラビン素数判定法 - Wikipedia

[2] 素因数分解の高速なアルゴリズム(ロー法) | 高校数学の美しい物語

クリーネの不動点定理(とベルマンフォード法)

言語処理系で講義でクリーネの不動点定理をやった時にベルマンフォードの証明にも使えるなぁと思ったので紹介します。

クリーネの不動点定理

 (D,\le) を完備半順序集合とし、 f をその上のスコット連続写像とする。このとき、 f は最小不動点を持つ。

ここで、完備半順序集合とは最小元  \bot を持ち、任意の有向部分集合  X について  X の上限  \sqcup X が存在する半順序集合とします。

 X が有向集合であるとは任意の  a,b \in X についてある元  c \in X が存在して  a \le c \land b \le c が成り立つこととします。

 f D 上のスコット連続写像であるとは、 D の任意の有向部分集合  X について  f(\sqcup X) = \sqcup \{ f(x) \mid x \in X \} が成り立つこととします。

何言ってんだコイツと思った人もブラウザバックを早まらないでください。

あとで有限半順序の場合を紹介しますが、有限の場合は簡単で、ベルマンフォードの証明に応用するのはそこだけで十分なので、意味がわからないひとは少しスクロールしてみてください。

とりあえず一般の完備半順序集合についてクリーネの不動点定理を証明しましょう。

証明

まず、 f がスコット連続なら  f は単調です。

なぜなら  X = \{ a, b \}, a \le b とすると、 X は有向集合で、  \sqcup X = b となり、スコット連続性を考えると  f(\sqcup X) = f(b) = \sqcup \{ f(a), f(b) \} となるので、 f(a) \le f(b) が導かれるからです。

 A = \{ f^{n}(\bot) \mid n \in \mathbb{N} \} とします。( \mathbb{N} 0 を含むものとします)

このとき、最小性より  \bot \le f(\bot) が成り立ち、単調性より  \bot \le f(\bot) \le f^{2}(\bot) \le ... となります。

よって、 A は有向集合で、完備性より上限  \sqcup A が存在します。これが目的の最小不動点となることを示します。

まず不動点性ですが、スコット連続性より

 f(\sqcup A) = \sqcup \{ f(f^{n}(\bot)) \mid n \in \mathbb{N} \} = \sqcup \{ f^{n}(\bot) \mid 1 \le n \in \mathbb{N} \} = \sqcup A より不動点です。

また、任意の不動点  x について、 \bot \le x と単調性より任意の  n \in \mathbb{N} について  f^{n}(\bot) \le f^{n}(x) = x となり、 \sqcup A \le x が成立します。

以上より示せました。

有限半順序の場合

おまたせしました。

有限の場合、クリーネの不動点定理は

 (D,\le) を最小元  \bot を持つ有限半順序集合とし、 f をその上の広義単調増加写像とする。このとき、 f は最小不動点を持つ。

となります。

 f が広義単調増加写像であるとは、任意の  a, b \in D に対して  a \le b \Rightarrow f(a) \le f(b) が成り立つことです。

それでは証明しましょう(一般の場合を読んだ人は繰り返しになるので読む必要は無いかもしれません)

有限半順序の場合の証明

 D の要素数 |D| と表します。

任意に  x \in D を一つ取ります。

 \bot, f(\bot), ..., f^{|D|}(\bot) を考えると、最小性より  \bot \le f(\bot) が成り立ち、単調性より  \bot \le f(\bot) \le ..\ \le f^{|D|}(\bot) が成り立ちます。

そして、鳩の巣原理から  0 \le i <  j \le |D| f^{i}(\bot) = f^{j}(\bot) となるものが存在し、  f^{i}(\bot) \le f^{i+1}(\bot) \le ...\le f^{j}(\bot) = f^{i}(\bot) の間が全て等号で成り立つので、  f^{i}(\bot) = \bot'不動点になります。

次に、任意に不動点  y を取ります。

 \bot の最小性より  \bot \le y です。

また、単調性より、 \bot' = f^{i}(\bot) \le f^{i}(y) = y となるので、  \bot'不動点の中で最小になります。

ベルマンフォード法について

負の閉路が無い場合のベルマンフォード法の正当性をクリーネの不動点定理で証明します。

双対(ポテンシャル)を考えます。

 D = \{ 0 \} \times \{ -{\rm INF}, ..., -1, 0, 1, ..\, {\rm INF} \}^{|V|-1} とします。

ここで、 {\rm INF} は十分大きい数(例えば辺の重みの絶対値の合計)とします。

順序は  (a_{1}, ..., a_{n}) \le (b_{1}, ..., b_{n}) \Leftrightarrow a_{1} \le b_{1} \land ... \land a_{n} \le b_{n} と定めます。

 D の元はポテンシャルを表しています。

最初の  0 は開始頂点のポテンシャルを  0 に固定しているということです。

 f をベルマンフォードの一回のループによるポテンシャルの変化とします。

 f を適用するとポテンシャルは広義単調減少します。

よって、クリーネの不動点定理より、最大元  \top = (0,{\rm INF}, ..., {\rm INF}) から  f を有限回適用すると、最大不動点  \top' に到達します。

不動点に到達すると、変化が起こらなくなるので、ベルマンフォードではループを抜けます。

このとき、ベルマンフォードの更新式から  \top' は LP の条件を満たします。

ここで、LP の条件を満たす任意の  D の元を  x とすると、これはベルマンフォードの処理  f不動点になっているはずです。

よって、最大性から  \top' \ge x となり、 \top' が最適です(最短路問題の双対なので、求めるものは最大値となります)

最後に

なんか回りくどいだけになってしまった気もします。(双対を考えればベルマンフォードの正当性はすぐに分かるので)

クリーネの不動点定理は変化しなくなるまで変化させ続けるような他のアルゴリズムについても停止性や正当性の証明に使えると思うので、雰囲気を感じ取ってもらえれば幸いです。

何か間違いや質問などがあればお願いします。

最後まで読んでいただきありがとうございました。

参考文献

[1] 横内 寛文. (1994). プログラム意味論 (情報数学講座).