競技プログラミングには概念を知っておかないと解きようがない、いわゆる覚えゲーのような問題が存在します。典型的な例が 10^9+7 といった素数で割った余りを求めろといったもので、普段業務で日常的に素数で割った余りを求めている人でもなければ、割り算がしたければフェルマーの小定理や拡張ユークリッドの互除法を使えば良いと直ぐには思い付けないのではないでしょうか。

　最小共通祖先も覚えゲーで必要な概念の一種と言えます。これは読んで字のごとく、与えられた根付き木の下で2頂点に共通する祖先のうち、最も根から遠い頂点を指す概念で、例えば木の2頂点が与えられて、頂点間の経路について何かを求めろといった問題で威力を発揮することが多いです。これを用いて解ける例を挙げるとすると次の問題でしょうか。 https://atcoder.jp/contests/abc014/tasks/abc014_4

　最小共通祖先を求めるアルゴリズムとして有名なのは、後に解説するダブリングと二分探索を用いるものですが、個人的にはこのアルゴリズム、バグを仕込み易く苦手に感じています。ダブリング自体は簡単なのですが、二分探索を正しく実装することの難しさと言ったら、競技プログラミングの経験がある方なら同意して頂けるのではないでしょうか。

　本記事では、最小共通祖先を求めるアルゴリズムの正しさを、後に解説する定理証明支援系 Coq を用いて数学的に証明し、バグが無いことを確認したライブラリを用いて上に挙げた AtCoder の問題を解きます。この記事によって、定理証明支援系や形式的検証、ひいては競技プログラミングに少しでも興味を持っていただければ幸いです。

最小共通祖先を求めるアルゴリズム

　最小共通祖先を効率的に求めるアルゴリズムは、ダブリングと二分探索を用いるものやオイラーツアーと RMQ を用いるものなど複数知られていますが、特別なデータ構造を用意する必要がないのもあって、取っつきやすいのは前者のアルゴリズムの方ではないでしょうか。

ダブリングとは

　二分探索は有名ですが、ダブリングは聞いたことが無い方も多そうなので簡単に解説しておくと、有限集合 S 上の関数 f : S \mapsto S に対して、事前処理に時間計算量 O(|S| log n) を費やすことで、繰り返し2乗法の要領で f^n (x) を時間計算量 O(log n) で計算するアルゴリズムです。そもそも繰り返し2乗法はモノイドの累乗に一般化できますし、恒等関数と関数合成はモノイドをなすので、ダブリングは繰り返し2乗法の一種と言えなくもないですが。

　前処理では、log n までの整数 k について f^2^k の返り値全てを求めておきます。関数 f^2 は関数 f 同士を合成したものと見做せますから、全ての引数に対して関数 f の返り値が分かっていれば、関数 f^2 (x) は f (f (x)) で計算できますよね。以下の OCaml のコードのように、それを繰り返せば前処理は完了です。

(* f は {0, 1, ..., n - 1} 上の関数で、mat.(0).(x) には f(x) の値が入っているものとする
 * m は log n を超える整数 *)
for i = 0 to m - 1 do
  for j = 0 to n - 1 do
    (* f^2^{k+1} (x) = f^2^k (f^2^k (x)) *)
    mat.(i + 1).(j) <- mat.(i).(mat.(i).(j))
  done
done
(* for 文を実行すると、mat.(k).(x) には f^2^k (x) の値が入る *)

　与えられた自然数 n に対して f^n (x) を求める際は、n を a_{m-1} 2^{m-1} + a_{m-2} 2^{m-2} + ... + a_0 といった具合に2進数展開し、a_i が 1 の場合は予め計算しておいた f^2^i を適用することを繰り返せば良いです。例えば 10 = 8 + 2 = 2^3 + 2^1 ですから、f^10 (x) は f^8 (f^2 (x)) = f^2^3 (f^2^1 (x)) で計算できますよね。この処理を OCaml で実装するとしたらこんな感じでしょうか？

let rec power_aux k n x =
  if k < 0
  then x
  (* n が 2^k より小さい場合、f^2^k は合成しなくて良い *)
  else if n < 1 lsl k
  then power_aux (k - 1) n x
  (* n が 2^k 以上の場合、f^{n - 2^k} (f^2^k (x)) で f^n (x) を計算する *)
  else power_aux (k - 1) (n - (1 lsl k)) mat.(k).(x)

(* f^n (x) を計算する関数 *)
let power n x = power_aux (m - 1) n x

最小共通祖先の計算

　ダブリングと二分探索を用いて最小共通祖先を求める際の基本的なアイデアは、高さが同じ頂点同士であれば、両方の頂点から同時に親に移動するのを繰り返して、両者が最初に同じになった頂点が最小共通祖先になるということです。高さが違う2頂点が与えられた際は、高い方の頂点から親を辿って行って高さを合わせてやれば良いです。

　ここで述べたアイデアをナイーブに実装しては、頂点数 |V| の根付き木に対して最小共通祖先を求めると毎回 O(|V|) の時間計算量が必要になってしまいそうですが、二分探索と先程紹介したダブリングを用いた場合、根付き木が与えられた段階で時間計算量 O(|V| log |V|) の前処理を行っておけば、与えられた2頂点ごとに O(log |V|) の時間計算量で最小共通祖先を求められます。高い方の頂点から親を辿って行って高さを合わせる操作は、与えられた頂点の親を返す関数を特定の回数適用していると考えられますから、いかにもダブリングで高速化できそうですね。また、同時に親に移動するのを繰り返して得られた2頂点が共通しているかは、同じ頂点から親を辿っても同じ頂点にしか行かないので単調性を満たす性質になっていて、最初に共通になった頂点は二分探索で高速に求められます。

　ここまでで紹介した最小共通祖先を求めるアルゴリズムを OCaml で実装すると、以下の通りです。

(* 根付き木の頂点は 0 から n - 1 までの整数
 * m は log n を超える整数
 * parent は与えられた頂点の親を返す関数
 * height は与えられた頂点の高さを返す関数
 * 根にあたる頂点は r *)
(* 2^i 回親を辿った際の頂点を覚えておくための2次元配列 *)
let mat = Array.make_matrix (m + 1) n 0 in
(* v から 2^0 回親を辿った頂点は、当然 v の親 *)
for v = 0 to n - 1 do
  mat.(0).(v) <- parent v
done;
(* ダブリングの前処理 *)
for i = 0 to m - 1 do
  for v = 0 to n - 1 do
    mat.(i + 1).(v) <- mat.(i).(mat.(i).(v))
  done
done;
let rec ancestor_aux i n v =
  if i < 0
  then v
  else ancestor_aux (i - 1) n @@ if n land (1 lsl i) = 0 then v else mat.(i).(v) in
(* 与えられた頂点から n 回親を辿った頂点をダブリングで求める関数 *)
let ancestor = ancestor_aux m in
let rec bisect_aux i u v =
  if i <= 0
  then parent u
  else
    let u' = mat.(i).(u) in
    let v' = mat.(i).(v) in
    if u' = v'
    then bisect (i - 1) u v
    else bisect (i - 1) u' v' in
(* 高さこそ同じなものの異なる2頂点が与えられた時、最小共通祖先を二分探索で求める関数 *)
let bisect = bisect_aux m in
(* 与えられた2頂点の最小共通祖先を返す関数 *)
let lca u v =
  let (u, v) =
    if height u <= height v
    then (u, ancestor (height v - height u) v)
    else (ancestor (height u - height v) u, v) in
  if u = v
  then u
  else bisect u v in

定理証明支援系 Coq とは？

　フランス国立情報学自動制御研究所で開発が進められている、数学の定理証明を支援するソフトウェアです。

　証明なんて紙とペンさえあれば書けるのに、何故ソフトウェアの助けを借りる必要があるのかと思うかもしれませんが、自然言語で書かれた証明は往々にして間違うもので、正しさを担保するために査読の形で人的資源が費やされてきました。これは言わば動的型付き言語でプログラムを書いて、バグが無いか目で確認するようなものですが、最初から静的型付き言語で書いておけば手間が省けるのにと思うことでしょう。定理証明支援系は、数学の証明の世界にも型チェッカを導入したようなものと言えます。

　定理証明支援系は Coq 以外にも Isabelle/HOL、Agda、HOL、Lean など色々ありますが、Coq の特色としてカリー・ハワード同型対応に基づいて設計されていることが挙げられます。

　カリー・ハワード同型対応とは、数学の世界における命題及び証明と、プログラミング言語の世界における型及びプログラムの間に存在する、一対一の対応関係です。要するに数学の命題はプログラミング言語の型に、数学の証明はプログラムに対応していると言っていて、例えば A ならば B のような命題は A を受け取って B を返す関数の型に対応します。

　Coq はこのカリー・ハワード同型対応に基づいて設計されていて、単なる定理証明支援系ではなく、述語論理も表現できるくらい強力な型システムを持った関数型言語と見ることもできます。プログラムは証明に対応するので Coq で記述できるし、そのプログラムの性質を Coq 上で証明することもできるのです。

最小共通祖先を求めるアルゴリズムの形式検証

　定理証明支援系 Coq 上ではプログラムの記述やその性質の証明ができると言われれば、実際にアルゴリズムの正しさを証明したくなるのが人情でしょう。ちょうど最小共通祖先を求めるアルゴリズムの実装が正しいのか悩んでいた所でした。

Coq 上での最小共通祖先の定式化

　最小共通祖先を求めるアルゴリズムの正しさを証明する前に、そのアルゴリズムを Coq 上で実装しなければなりませんし、その前に最小共通祖先とはどう言うものか、そもそも根付き木とは何かから定義しなくてはなりません。

　一般的な木の定義と言えば連結な閉路を持たないグラフで、根付き木は根 (root) と呼ばれる頂点が存在する木ですが、本記事では最小共通祖先を求めるアルゴリズムの実装や検証を簡単にするために、以下のような定義を採用します。

(* 頂点を表す型 *)
Variable V : eqType.
(* 根にあたる頂点 *)
Variable root : V.
(* 与えられた頂点の親を返す関数 *)
Variable parent : V -> V.
(* 与えられた頂点の高さを返す関数 *)
Hypothesis height : V -> nat.
(* 根にあたる頂点、親を返す関数と、高さを返す関数が満たすべき公理 *)
Hypothesis height_spec : forall n v, (iter n parent v == root) = (height v <= n).

頂点を表す型があって、根があって、親を返す関数や高さを返す関数が定義されている…と一見当たり前の事実しか書かれていないように見えますが、一つ公理が要求されていますね。この公理のステートメントを自然言語で書くと「ある頂点から n 回親を辿ると根に行き着くことと、n がその頂点の高さ以上であることは同値である」と言ったものですが、実はこの公理から根付き木の基本的な性質が全て導けます。

　例えば、根にあたる頂点の高さは 0 ですし、逆に高さが 0 である頂点は根しか存在しません。

Lemma height_root_eq v : (height v == 0) = (v == root).
Proof. by rewrite -leqn0 -height_spec. Qed.

この性質は、公理のステートメント「ある頂点から n 回親を辿ると根に行き着くことと、n がその頂点の高さ以上であることは同値である」の n を 0 で置き換えると「ある頂点が根であることと、頂点の高さが 0 であることは同値である」となるので直ちに得られます。（頂点の高さは自然数なので、高さが 0 以下であることと高さが 0 であることは同値なことに注意してください）

　また、子にあたる頂点の高さは、親の高さに 1 足したものになります。

Lemma height_parent v : height (parent v) = (height v).-1.
Proof.
  case (leqP (height v) 0) => [ | /prednK Hpredn ].
  - rewrite leqn0 height_root_eq => /eqP ->.
    by rewrite parent_root height_root.
  - apply /anti_leq /andP.
    split.
    + by rewrite -height_spec -iterSr Hpredn height_spec.
    + by rewrite -ltnS Hpredn -height_spec iterSr height_spec.
Qed.

この性質は反対称律を使って高さの大小関係二つに落とし込んで、さらに公理によって、親を辿って行った時に根にたどり着くかどうかに帰着する感じで証明できます。

　加えて、最小共通祖先を求めるアルゴリズムの検証に使わなかったので Coq で証明はしていませんが、ここで定義した根付き木には（根から根への自己辺を除いて）閉路が存在しません。これは何故かと言うと、公理から全ての頂点は有限回親を辿ると根に行き着くことが保証されていますが、根を通らない閉路が存在すると、その閉路上の頂点は何回親を辿っても根に辿り着けないことになってしまうからです。（parentの型から、ある頂点の親はただ一つしか存在しないことが導かれることに注意してください）

　ちなみに根の親は根になっていて厳密には閉路が存在しているのですが、これは parentの定義上どんな頂点にも親が必ず一つ存在することにしたかったので、番兵のようなものです。

Corollary parent_root : parent root = root.
Proof.
  move: (height_spec 1 root).
  by rewrite /= height_root => /eqP ->.
Qed.

　Coq 上での根付き木の定義がちゃんとグラフ理論の根付き木と対応していることを簡単に確かめたところで、祖先の定義をしていきましょう。頂点 u が頂点 v の祖先であることを表す二項関係は、v から u と同じ高さになるまで親を辿って行くと u と一致するかどうかで定義できます。

Definition ancestor u v := iter (height v - height u) parent v == u.

　明らかに頂点 v の親は v の祖先ですし、v から何回か親を辿って行った頂点も v の祖先ですね。

Lemma ancestor_iter_parent n v : ancestor (iter n parent v) v.
Proof.
  rewrite /ancestor height_iter_parent.
  case (leqP n (height v)) => [ /subKn -> // | /ltnW Hleq ].
  have -> : height v - n = 0 by apply /eqP; rewrite subn_eq0.
  move: height_spec (Hleq) => <- /eqP ->.
  by rewrite subn0 height_spec.
Qed.

Corollary ancestor_parent v : ancestor (parent v) v.
Proof. exact /(ancestor_iter_parent 1). Qed.

（明らかに成り立ってそうな性質を何故証明するのかと首を傾げるかもしれませんが、Coq 上での定義が我々の定義しようとしている概念と対応しているか、仕様のデバッグを行うためにも、直感的に成り立ってほしい性質を証明するのは形式的検証では大事なことです）

　頂点 u が頂点 v の祖先であることを表す二項関係さえ定義できてしまえば、頂点 w が頂点 u, v の共通の祖先であることを表す三項関係は直ちに定義できますし、

Definition common_ancestor u v w := ancestor w u && ancestor w v.

「頂点 u と 頂点 v の共通祖先である」という性質を満たす最小の頂点になっていることを以って、頂点 w が u と v の最小共通祖先であることを定義できます。

Definition lowest_common_ancestor u v w :=
  forall w', common_ancestor u v w' = ancestor w' w.

最小共通祖先を求めるアルゴリズムの実装

　根付き木、祖先、最小共通祖先と言った概念の定義が終わって「地固め」が済んだところで、Coq 上で最小共通祖先を求めるアルゴリズムの実装を行いましょう。

　まず、高さが同じになるまで高い方の頂点の親を辿って行く操作で必要になる、与えられた回数だけ親を辿って行く関数を定義します。これはダブリングで実装できる訳ですが、Coq は副作用を扱うのが苦手なので、前処理は既に済んでいて 2^i 回親を辿った頂点は全て求まっているものとします。

(* 2^i 回親を辿った頂点を返す関数 *)
Variable climb2exp : nat -> V -> V.
Hypothesis climb2exp_spec : forall n v, climb2exp n v = iter (2 ^ n) parent v.

この前処理の結果を用いると、log n より十分大きな自然数 e と 2^e に相当する自然数 p が与えられた時に、頂点 v から n 回親を辿って行き着く頂点を返す関数 climb e p n v は以下のように定義できます。

Fixpoint climb e p n v :=
  if e is e'.+1
  then
    let p' := p./2 in
    if n < p'
    then climb e' p' n v
    else climb e' p' (n - p') (climb2exp e' v)
  else v.

左シフトに相当する Coq の関数を見付けられなかったので引数が増えてしまいましたが、最小共通祖先の解説をした時に書いた OCaml のコードと殆ど同じですね。

　また、高さは同じであるものの異なる2頂点 u, v と、それらの高さに log を取ったものより十分大きな自然数 e が与えられた時に、二分探索で u と v の最小共通祖先を求める関数 bisect e u v は以下のように実装できます。

Fixpoint bisect e u v :=
  if e is e'.+1
  then
    let u' := climb2exp e' u in
    let v' := climb2exp e' v in
    if u' == v'
    then bisect e' u v
    else bisect e' u' v'
  else climb2exp 0 u.

　これらを用いると、2頂点 u, v と、それらの高さに log を取ったものより十分大きな自然数 e、及び 2^e に相当する自然数 p が与えられた時に、u と v の最小共通祖先を返す関数 lca e p u v は次のように実装できます。

(* 高さが同じ2頂点の最小共通祖先を求める関数 *)
Definition lca_aux e u v :=
  if eq_vertex u v
  then u
  else bisect e u v.

Definition lca e p u v :=
  if height u <= height v
  then lca_aux e u (climb e p (height v - height u) v)
  else lca_aux e v (climb e p (height u - height v) u).

アルゴリズムの形式検証

　与えられた2頂点の最小共通祖先を返す関数 lca を Coq 上で実装できたので、その正しさ、lca の返り値が2頂点の最小共通祖先であることを証明していきましょう。lca が最小共通祖先を返すことを Coq のステートメントで書くと以下の通りですが、これを証明するには、lca の実装に用いた climb、bisect といった種々の関数の性質の証明が必要になることは言うまでもありません。

Theorem lca_spec : forall e u v,
  height u < 2 ^ e ->
  height v < 2 ^ e ->
  lowest_common_ancestor u v (lca e (2 ^ e) u v).

　まず、log n より十分大きな自然数 e と 2^e に相当する自然数 p が与えられた時に、頂点 v から n 回親を遡って辿り着く頂点を返す関数 climb e p u v の満たすべき性質は、タブリングを使わずにナイーブに n 回遡る実装の返す値と一致することです。これに関しては割と簡単に示せます。

Lemma climb_spec : forall e n v,
  n < 2 ^ e ->
  climb e (2 ^ e) n v = iter n parent v.
Proof.
  elim => /= [ ? ? | e IH n ? ].
  - by rewrite expn0 -subn_eq0 subn1 succnK => /eqP ->.
  - rewrite expnS mul2n half_double -addnn => Hltn.
    case (leqP (2 ^ e) n) => [ Hleq | /IH // ].
    rewrite climb2exp_spec IH.
    * by rewrite -iterD subnK.
    * by rewrite ltn_subLR.
  Qed.

　次に、高さは同じであるものの異なる2頂点 u, v と、それらの高さに log を取ったものより十分大きな自然数 e が与えられた時に、二分探索で u と v の最小共通祖先を求める関数 bisect e u v の満たすべき性質は、bisect の返り値が u と v の最小共通祖先であることです。最小共通祖先を求めるアルゴリズムの本質的な処理なのもあって、割と証明は重めですね。

Lemma bisect_spec : forall e u v,
  u != v ->
  height u = height v ->
  iter (2 ^ e) parent u = iter (2 ^ e) parent v ->
  lowest_common_ancestor u v (bisect e u v).
Proof.
  elim => /= [ u v Huv Hheight Hparent w' | e IH u v Huv Hheight ].
  - rewrite climb2exp_spec /= /common_ancestor /ancestor -Hheight height_parent -predn_sub.
    case (posnP (height u - height w')) => [ Hzero | /prednK <- ].
    { rewrite Hzero /=.
      apply /eq_bool.
      split => [ /andP [ /eqP <- /eqP Heq ] | /eqP /(f_equal height) ].
      - by move: Heq eq_refl Huv => -> ->.
      - rewrite height_parent => Hw'.
        move: subn1 Hw' Hzero => <- <-.
        case (posnP (height u)) => [ /eqP Hroot | /subKn -> // ].
        move: height_root_eq (Hroot) Huv => -> /eqP ->.
        by move: Hheight height_root_eq Hroot eq_refl => -> -> /eqP -> ->. }
    by rewrite succnK !iterSr Hparent andbb.
  - rewrite expnS mul2n -addnn !iterD !climb2exp_spec.
    case (vertex_eqP (iter (2 ^ e) parent u) (iter (2 ^ e) parent v)) => [ ? ? | Hneq Heq w' ].
    + exact /IH.
    + have /(_ w') <- : lowest_common_ancestor
        (iter (2 ^ e) parent u) (iter (2 ^ e) parent v)
        (bisect e (iter (2 ^ e) parent u) (iter (2 ^ e) parent v)).
      { apply /IH => //.
        - exact /eqP.
        - by rewrite !height_iter_parent Hheight. }
      rewrite /common_ancestor /ancestor !height_iter_parent -!iterD !(subnAC _ (2 ^ e)) -Hheight.
      case (leqP (2 ^ e) (height u - height w')) => [ /subnK -> // | /ltnW Hleq ].
      move: subn_eq0 (Hleq) add0n Hneq => <- /eqP -> -> Hneq.
      apply /eq_bool.
      split => [ | /andP [ /eqP <- /eqP Hcontra ] ].
      { move: (subnK Hleq) Hneq => <-.
        rewrite !iterD => Hneq /andP [ ].
        rewrite eq_sym => /eqP /(eq_trans _) Htrans /eqP /Htrans Hcontra.
        by move: Hcontra (Hneq) => -> /(_ erefl). }
      by move: Hcontra Hneq => -> /(_ erefl).
Qed.

　これらの関数の性質さえ証明してしまえば、主定理である lca が最小共通祖先を返すことは、補題を適用するのに必要な仮定の部分を補う程度で証明できてしまいます。そもそもダブリングと二分探索で最小共通祖先を求めるアルゴリズム自体、二分探索の所が本質で、それ以外の部分は二分探索を適用できるようにするためのお膳立てみたいなところがあるので当たり前と言えば当たり前なのですが。

Lemma lca_aux_spec : forall e u v,
  height u = height v ->
  height u <= 2 ^ e ->
  lowest_common_ancestor u v (lca_aux e u v). 
Proof.
  move => e u v Hheight Hu. 
  rewrite /lca_aux.
  case (vertex_eqP u v) => [ <- ? | /eqP ? ].
  - by rewrite /common_ancestor andbb.
  - apply /bisect_spec => //. 
    move: height_spec (Hu) => <- /eqP ->. 
    by move: Hheight height_spec Hu => -> <- /eqP ->. 
Qed.

Definition lca e p u v :=
  if height u <= height v
  then lca_aux e u (climb e p (height v - height u) v)
  else lca_aux e v (climb e p (height u - height v) u)

Theorem lca_spec : forall e u v,
  height u < 2 ^ e ->
  height v < 2 ^ e ->
  lowest_common_ancestor u v (lca e (2 ^ e) u v). 
Proof.
  move => e u v Hu Hv. 
  rewrite /lca !climb_spec.
  - case (leqP (height u) (height v)) => [ /minn_idPr ? | /ltnW /minn_idPl ? ].
    + apply /lowest_common_ancestor_climb /lca_aux_spec.
      * by rewrite height_iter_parent -minnE.
      * exact /ltnW.
    + apply /lowest_common_ancestor_comm /lowest_common_ancestor_climb /lca_aux_spec.
      * by rewrite height_iter_parent -minnE minnC.
      * exact /ltnW.
  - exact /(leq_ltn_trans (leq_subr _ _)).
  - exact /(leq_ltn_trans (leq_subr _ _)).
Qed.

正当性を証明したライブラリの使用例

　定理証明支援系 Coq の良いところは、Coq 上でプログラムの実装も、その性質の証明も行えることもさることながら、Coq で書いたプログラムを OCaml、Haskell、Scheme といった実用的な言語のプログラムに変換できる所にあります。折角最小共通祖先を求めるプログラムを実装して、実装の正しさを数学的に検証までしたのに Coq の中だけで遊ばせておくのも勿体ないですから、OCaml のコードに変換して実際に競技プログラミングの問題を解くのに使ってみましょう。

　Coq で書いたプログラムを他の言語に変換するのはそう難しい手順が必要な訳ではなく、単に Extraction というコマンドを実行するだけで良いです。例えば、今回実装した最小共通祖先を求める関数 lca を OCaml のコードに変換したければ

Extraction lca.

と書くだけで

(** val climb : (int -> 'a1 -> 'a1) -> int -> int -> int -> 'a1 -> 'a1 **)
 
let rec climb climb2exp e p n v =
  (fun fO fS n -> if n = 0 then fO () else fS (n-1))
    (fun _ -> v)
    (fun e' ->
    let p' = half p in
    if ( <= ) (succ n) p'
    then climb climb2exp e' p' n v
    else climb climb2exp e' p' (( - ) n p') (climb2exp e' v))
    e
 
(** val bisect :
    ('a1 -> 'a1 -> bool) -> (int -> 'a1 -> 'a1) -> int -> 'a1 -> 'a1 -> 'a1 **)
 
let rec bisect eq_vertex climb2exp e u v =
  (fun fO fS n -> if n = 0 then fO () else fS (n-1))
    (fun _ -> climb2exp 0 u)
    (fun e' ->
    let u' = climb2exp e' u in
    let v' = climb2exp e' v in
    if eq_vertex u' v'
    then bisect eq_vertex climb2exp e' u v
    else bisect eq_vertex climb2exp e' u' v')
    e
 
(** val lca_aux :
    ('a1 -> 'a1 -> bool) -> (int -> 'a1 -> 'a1) -> int -> 'a1 -> 'a1 -> 'a1 **)
 
let lca_aux eq_vertex climb2exp e u v =
  if eq_vertex u v then u else bisect eq_vertex climb2exp e u v
 
(** val lca :
    ('a1 -> 'a1 -> bool) -> ('a1 -> int) -> (int -> 'a1 -> 'a1) -> int -> int
    -> 'a1 -> 'a1 -> 'a1 **)
 
let lca eq_vertex height climb2exp e p u v =
  if ( <= ) (height u) (height v)
  then lca_aux eq_vertex climb2exp e u
         (climb climb2exp e p (( - ) (height v) (height u)) v)
  else lca_aux eq_vertex climb2exp e v
         (climb climb2exp e p (( - ) (height u) (height v)) u)

のようなコードを生成してくれます。

　OCaml のように競技プログラミングのサイトでも使える言語での実装が手に入ったとなれば試したくなるもので、これを使って記事冒頭で紹介した最小共通祖先を用いる問題を実際に解いてみました。

https://atcoder.jp/contests/abc014/submissions/24778992

ダブリングの前準備であったり頂点数の log を取る部分であったりは、Coq 上で扱うのが面倒でサボっていたので OCaml 上でお膳立てをする必要がありましたが、それ以外の部分については Coq でバグがないことを確認し、全く不安に思う所もなく安心して AC できました。もちろんコンテスト中に Coq で証明を書くのは時間的に厳しいとは思いますが、自分でライブラリを整備するならこれほど心強い味方もないでしょう。

まとめ

　本記事では、最小共通祖先を求めるアルゴリズムの正しさを、定理証明支援系 Coq を用いて数学的に証明し、バグが無いことを確認したライブラリを用いて上に挙げた AtCoder の問題を解きました。今回 Coq で記述したプログラムや証明の全文は Gist にアップロードしてあるので、手元で動かしてみたい方があればそちらをご利用下さい。

　最後に少し宣伝ですが、最近僕は JavaScript で書かれた既存のフロントエンドに、TypeScript で型を付けることで生産性を改善する仕事をしています。流石に Coq ほど強力な型システムを採用した言語を業務で書く機会はありませんが、強い静的型付きの言語を愛好する人が来てくれると、型の付いたコードベースが増えて僕が喜びます。