今となっては少々昔のことですが，AtCoder Beginner Contest 167のD問題で数列の循環検出が必要になった事がありました（問題の制約的にはダブリングでも解けますが）．

数列の循環検出はナイーブにに実装すると，循環している部分に至るまでの長さを $\mu$ ，循環の長さを $\lambda$ とすると $O(\mu+\lambda)$ のメモリが必要になってしまいますが，これを $O(1)$ で計算するアルゴリズムとして，フロイドの循環検出法が知られています．中々巧妙なアルゴリズムなもので，これを実装して本当に循環検出ができているのか少々不安に思うところもありますし，この際Coqで正当性を保証してやりましょう．

フロイドの循環検出法

フロイドの循環検出法とは， $a_{i+1}=f(a_i)$ で与えられる数列の循環を空間計算量 $O(1)$ で検出し，循環している部分に至るまでの長さ $\mu$ と循環の長さ $\lambda$ を求めるアルゴリズムです．

正直Wikipediaより分かりやすく解説できる気がしませんが*1，具体的な手順を挙げると

$a_m$ を覚えておくための変数と $a_{2m}$ を覚えておくための変数を用意し， $a_{m+1}=f(a_m)$ と $a_{2(m+1)}=f^2(a_{2m})$ を利用して $a_m=a_{2m}$ となる ${0\lt m}$ を線形探索
$a_\mu$ を覚えておくための変数と $a_{\mu+m}$ を覚えておくための変数を用意し， $a_{\mu+1}=f(a_\mu)$ と $a_{(\mu+1)+m}=f(a_{\mu+m})$ を利用して $a_{\mu}=a_{\mu+m}$ となる最小の ${0\leq \mu}$ を線形探索
$a_{\mu+\lambda}$ を覚えておくための変数を用意し， $a_{\mu+(\lambda+1)}=f(a_{\mu+\lambda})$ を利用して $a_\mu=a_{\mu+\lambda}$ となる最小の ${0\lt \lambda}$ を線形探索

となります．1. で速く動く変数 $a_{2m}$ と遅く動く変数 $a_m$ を使うことから，「ウサギとカメのアルゴリズム」とも呼ばれているらしいです．

このアルゴリズムを循環検出アルゴリズムたらしめる重要な事実は，あらゆる循環の長さらしきもの—— $a_{i+m}=a_{i}$ となる $m$ ——は実際の循環の長さ $\lambda$ の整数倍となっていることです．それ故に $a_{\mu+m}=a_{\mu+\lambda\lfloor\frac{m}{\lambda}\rfloor}=a_{\mu}$ より， $a_{\mu}=a_{\mu+m}$ となる $\mu$ を線形探索することで循環している部分に至るまでの長さを求められ，それを足がかりとして循環の長さも求められるのです． $a_m=a_{2m}$ となる $m$ を線形探索しているのは，単に $a_{i+m}=a_{i}$ となる $m$ の一例を得たいからに過ぎませんが，循環している数列 $a_i$ に対しては必ず $a_m=a_{2m}$ となる $m$ が存在することが保証されます．*2

Coqによる検証

このアルゴリズムはCoq上で検証する以前に，実装するだけでも大変です．例えば $a_m=a_{2m}$ となる $m$ を探索する処理を単純に実装すると

Fixpoint tortoise_and_hare m a_m a_2m :=
  if a_m == a_2m
  then (m, a_m)
  else tortoise_and_hare (S m) (f a_m) (f (f a_2m)).

のような再帰関数を書くことになると思いますが，停止性が明らかではないとCoqに突っぱねられてしまいます． Coqの言い分ももっともで， $a_m=a_{2m}$ となる $m\gt 0$ が存在しなければこの関数は停止しませんから，そのような $m\gt 0$ が存在することを加味した型の関数にしなければなりませんね．

こういう停止性の怪しい再帰関数を3回も書かなければならないので，もう少し一般化して使い回しが効くような関数を定義してしまいましょう．

(* P(n)を満たす最小の自然数nを線形探索する *)
Lemma nat_eps (P : nat -> Prop)
  (Hdec : forall n, { P n } + { ~ P n })
  (Hex : exists n, P n) :
  { n | P n & forall n', P n' -> n <= n' }.
Proof.
  refine (@Fix _
    (fun n m => S m = n /\ (forall m, m < n -> ~ P m)) _
    (fun n => (forall m, m < n -> ~ P m) -> { n | P n & forall n', P n' -> n <= n' })
    (fun n eps HnP =>
       match Hdec n with
       | left _ => exist2 _ _ n _ (fun n' => _)
       | right _ => let Hlb := fun m => _ in eps (S n) _ Hlb
       end) 0 _); eauto.
  Unshelve.
  - case Hex => n Hp m.
    remember (n - m) as d.
    generalize dependent m.
    induction d as [ | ? ] => m ?; constructor => ? [ ? Hnp ]; subst.
    + (have : (n < S m) by lia) => /Hnp. congruence.
    + apply /IHd. lia.
  - by case (le_lt_dec n n') => // Hlt /(HnP _ Hlt).
  - lia.
  - move => /le_S_n ?.
    case (Nat.eq_dec m n) => [ -> // | ? ].
    apply /HnP. lia.
Defined.

このnat_epsなる関数は，ある自然数nに対しては真になることが分かっていて，適当な自然数に対して真か偽か判別できるような述語Pを受け取って，Pが真になるような最小のnを線形探索して返す関数です．肝心の処理は

Fixpoint nat_eps P n :=
  if P n
  then n
  else nat_eps P (S n).

ぐらいの単純なことしかしていませんが，再帰の停止性をCoqにも分かりやすく伝えるために面倒な記述が増えています．

ではこのnat_epsを使って，肝心のフロイドの循環検出アルゴリズムを実装していきましょう．

Section CycleDetection.
  Variable A : Set.
  Variable f : A -> A.
  Variable eq_dec : forall x y : A, { x = y } + { x <> y }.

  Definition eventually_periodic m0 l n : Prop :=
    l > 0 /\ forall m, m0 <= m -> iter (l + m) f n = iter m f n.

  Theorem tortoise_and_hare v
    (Hperiod : exists m l, eventually_periodic m l v /\ forall m' l', eventually_periodic m' l' v -> m <= m' /\ l <= l') :
    {'(m, l) | eventually_periodic m l v /\ forall m' l', eventually_periodic m' l' v -> m <= m' /\ l <= l' }.
  Proof.
    refine
      (let (m, Hm, _) :=
         (* ウサギとカメ *)
         nat_eps (fun m => @iter (S m + S m) _ f v = @iter (S m) _ f v)
                 (fun m => eq_dec (@iter (S m + S m) _ f v) (@iter (S m) _ f v)) _ in
       let (mu, Hmu, Hmu') :=
         (* ループの開始地点を見つける *)
         nat_eps (fun mu => @iter (S m + mu) _ f v = @iter mu _ f v)
                 (fun mu => eq_dec (@iter (S m + mu) _ f v) (@iter mu _ f v)) _ in
       let (l, Hl, Hl') :=
         (* ループの周期を測る *)
         nat_eps (fun l => @iter (S l + mu) _ f v = @iter mu _ f v)
                 (fun l => eq_dec (@iter (S l + mu) _ f v) (@iter mu _ f v)) _ in
       exist _ (mu, S l) _);

eventually_periodic m l vは， $a_0=v, a_{i+1}=f(a_i)$ で与えられる数列 $a_i$ が循環しており，循環に至るまでの長さがm，循環の長さがlであることを意味する3項関係です．ただし，我々が図で見てああ明らかに循環に至るまでの長さmだな，循環の長さlだなと思うものは，実際にはeventually_periodic m l vとなるような最小のm，lであることに注意が必要です．

フロイドの循環検出法は，何かしらの循環に至るまでの長さmと循環の長さlで循環している数列 $a_i$ を受け取って，そのmとlを具体的に計算する関数tortoise_and_hareとして実装されています．型だけ見るとそのままmとlを返せば良いように見えますが，Coqの型システムではProp型の型を持つ値をパターンマッチしてSet型の型を持つ値は作れないので，何かしらのmとlで循環していることの情報は，単にアルゴリズムの停止性を保証することだけに使われます．

証明とプログラムを一緒に書くスタイルにしてしまったので，このまま正当性の証明も書く必要があります．

    move: Hperiod => [ mu' [ l' [ [ Hpos Hperiod ] Hleast ] ] ]; eauto.
    { (* a_{m+1} = a_{2(m+1)}となるようなm>=0が存在することの証明 *)
      exists ((mu' + 2) * l' - 1).
      set j := (mu' + 2) * l' - 1.
      have->: S j + S j = (mu' + 2) * l' + S j by lia.
      elim: (mu' + 2) => [ // | ? ? /= ].
      rewrite -plus_assoc Hperiod //. nia. }
    have : eventually_periodic mu (S l) v => [ | /Hleast [ ? Hlel ] ].
    { split => [ | m0 ? ]; auto with arith.
      have->: m0 = (m0 - mu) + mu by lia.
      have->: S l + (m0 - mu + mu) = (m0 - mu) + (S l + mu) by lia.
      by rewrite !(@iter_plus _ _ _ _ (m0 - mu)) Hl. }
    have : eventually_periodic mu (S m) v => [ | /Hleast [ Hlemu Hlem ] ].
    { split => [ | m0 ? ]; auto with arith.
      have->: m0 = (m0 - mu) + mu by lia.
      have->: S m + (m0 - mu + mu) = (m0 - mu) + (S m + mu) by lia.
      by rewrite !(@iter_plus _ _ _ _ (m0 - mu)) Hmu. }
    (* m+1が循環の長さl'で割り切れるかどうかで場合分け *)
    case (zerop (S m mod l')) => Hmod.
    - have : mu <= mu' => [ | /(le_antisym _ _ Hlemu) ? ]; subst.
      { apply /Hmu'.
        rewrite (Nat.div_mod (S m) l'). { lia. }
        rewrite Hmod -plus_n_O mult_comm.
        elim: (S m / l') => //= ? IH.
        by rewrite -plus_assoc iter_plus IH -iter_plus Hperiod. }
      have : S l <= l' => [ | /(le_antisym _ _ Hlel) <- // ].
      { rewrite (S_pred_pos l') //.
        apply /le_n_S /Hl'.
        by rewrite -S_pred_pos // Hperiod. }
    - (* m+1が循環の長さl'で割り切れなかった場合，長さ(m+1)mod l'の循環が存在することになる
         これはl'が最も短い循環の長さであることに矛盾する *) 
      have : forall i, i * l' <= S m -> iter (S m - i * l' + mu) f v = iter mu f v => [ | /(_ (S m / l')) ].
      { elim => // i IH ?.
        rewrite -IH. { lia. }
        have->: S m - i * l' + mu = l' + (S m - S i * l' + mu) by lia.
        rewrite Hperiod //. lia. }
      rewrite (mult_comm _ l') -Nat.mod_eq => [ | Hperiod' ]. { lia. }
      have : eventually_periodic mu (S m mod l') v => [ | /Hleast [ ? /le_not_lt [ ] ] ].
      { split => [ // | m0 ? ].
        have->: m0 = (m0 - mu) + mu by lia.
        have->: S m mod l' + (m0 - mu + mu) = (m0 - mu) + (S m mod l' + mu) by lia.
        rewrite !(@iter_plus _ _ _ _ (m0 - mu)) Hperiod' //.
        apply /Nat.mul_div_le. lia. }
      apply /Nat.mod_upper_bound. lia.
  Defined.
End CycleDetection.

線形探索が停止することを保証するために $a_m=a_{2m}$ となる $m\gt 0$ が少なくとも1つは存在することを証明している部分と，そのような $m$ が循環の長さ $\mu$ で割り切れることを証明する部分が少しややこしいですが，あとは前の節で書いたようなノリですね．

まとめ

フロイドの循環検出アルゴリズムをCoqで実装するとともに，その実装で循環に至るまでの長さと循環の長さを求められることも証明しました．今回Coqで書いたコードの全文はGistで見られます．

循環があればこのアルゴリズムで検出できることはTopProverの過去問で実は既に証明していたのですが，循環に至るまでの長さ $\mu$ と循環の長さ $\lambda$ を具体的に求めることはできていませんでした．今回の記事では $\mu$ と $\lambda$ を計算する処理も含めて検証しています．

実は今回のフロイドの循環検出アルゴリズムの実装は $a_m$ や $a_{2m}$ などの値をメモ化しておらず，これをextractしようものなら甚だ非効率的なコードが出力されるのですが，実用に堪える実装とするのは今後の課題とします*3．

*1:Wikipediaを参考にするのであれば，英語版の記事の方が僕には参考になりました

*2:後ほどCoqで証明します

*3:つらいので

fetburner.core

コアダンプ

Coqでフロイドの循環検出法を検証する

フロイドの循環検出法

Coqによる検証

まとめ