以前二分探索の正当性をCoqで証明しましたが，この証明スクリプトからOCamlのコードを抽出して実際に競プロに用いようとすると，添字の動く範囲が保証されていないので添字エラーを起こすこともしばしばでした．そこで，依存型を活用して添字の動く範囲も保証しようと思います．

以前の実装

以前Coqで書いた二分探索は以下のように，まずプログラムを書いてから

Variable P : nat -> Prop.
Hypothesis P_dec : forall n, { P n } + { ~P n }.

Function lower_bound_aux offset n { wf lt n } :=
  match n with
  | O => offset
  | _ =>
    let m := Nat.div2 n in
    if P_dec (offset + m) then lower_bound_aux offset m
    else lower_bound_aux (S m + offset) (n - S m)
  end.
Proof.
  - intros. apply Nat.lt_div2. omega.
  - intros. omega.
  - apply lt_wf.
Defined.

Definition lower_bound alpha beta :=
  lower_bound_aux alpha (beta - alpha).

満たすべき性質を証明するものでした．

Variable threshold : nat.
Hypothesis P_monotone1 : forall n, threshold <= n -> P n.
Hypothesis P_monotone2 : forall n, P n -> threshold <= n.

Lemma lower_bound_aux_spec n : forall offset error,
  threshold = offset + error ->
  error <= n ->
  lower_bound_aux offset n = offset + error.
Proof.
  induction n as [[| n'] IHn] using lt_wf_ind;
    intros offset error ? ?;
    subst;
    rewrite lower_bound_aux_equation.
  - omega.
  - destruct (P_dec (offset + Nat.div2 (S n'))) as [ HP | ].
    + apply IHn; eauto.
      * apply Nat.lt_div2.
        omega.
      * specialize (P_monotone2 _ HP).
        omega.
    + destruct (le_dec (offset + error) (offset + Nat.div2 (S n'))) as [ Hle | ].
      * apply P_monotone1 in Hle.
        congruence.
      * rewrite IHn with (error0 := error - S (Nat.div2 (S n'))); try omega.
Qed.

Theorem lower_bound_spec alpha beta :
  alpha <= threshold <= beta ->
  lower_bound alpha beta = threshold.
Proof.
  unfold lower_bound.
  intros ?.
  rewrite lower_bound_aux_spec with (error := threshold - alpha); try omega.
Qed.

問題点

主定理lower_bound_specを見ると，添字の動く範囲には一切言及されていないことが分かりますね？

このようにプログラムと証明を分けて書くスタイル，入力と出力についての性質にしか言及しないなら良いんですが，添字の動く範囲なんてのはプログラムの中で何をやっているかに相当する性質なので，これを示そうと思うとちょっと困ってしまいます．

解決法

これを解決するために，述語pにforall n, l < n <= r -> boolといった感じの依存型を付けてやりましょう．このpを呼び出すためには添字nに加えて，nがl < n <= rを満たすことの証明も与える必要がありますから， lower_bound l r pの中でpが呼び出される時，添字は必ず(l, r]の半開区間に収まっていると保証されます．

新しい実装

pが証明を受け取るようになったせいで，図らずしもプログラムと証明を混ぜて書くスタイルになってしまいました．せっかくなので正当性の証明もプログラムと一緒に書くとすれば，二分探索の実装は以下のようになります．

Program Fixpoint lower_bound' l r
  (p : forall n, l < n <= r -> bool)
  (Hacc : @Acc _ (ltof _ (fun p => snd p - fst p)) (l, r))
  (_ : exists m,
    l < m <= r /\
    (forall n H, p n H = true -> m <= n) /\
    (forall n H, m <= n -> p n H = true)) :
  { m | l < m <= r /\
    (forall n H, p n H = true -> m <= n) /\
    (forall n H, m <= n -> p n H = true) } :=
  match Hacc with
  | Acc_intro _ Hacc' =>
      if le_gt_dec r (S l)
      then exist _ r _
      else
        let m := (l + r) / 2 in
        ( if p m _ as b return p m _ = b -> _
          then fun _ =>
            let (n, _) := lower_bound' l m (fun n _ => p n _) (Hacc' _ _) _ in
            exist _ n _
          else fun _ =>
            let (n, _) := lower_bound' m r (fun n _ => p n _) (Hacc' _ _) _ in
            exist _ n _ ) eq_refl
  end.

停止性を保証するためにAccなるものを使っていますが，よく読めばこれと同じ実装になっていると気付くことでしょう．

【めぐるのアルゴリズム講座】
二分探索（整数）の書き方
難しさ：４ pic.twitter.com/LGLbkS0D7l
— 因幡めぐる@競技プログラミング (@meguru_comp) 2016年2月9日

さて，Programコマンドによって証明に相当する項を虫食いのようにしてプログラムを書きましたが，それらは後で与えてやる必要があります． pが呼び出される度に証明が必要ですから，これまた証明すべき命題が多くなってしまいますが…

Next Obligation.
  repeat split; intros; try omega.
  eapply H3. omega.
Qed.
Next Obligation.
  split.
  - apply (Nat.div_le_lower_bound _ 2); omega.
  - apply (Nat.div_le_upper_bound _ 2); omega.
Qed.
Next Obligation.
  assert (l < (l + r) / 2).
  { apply Nat.div_le_lower_bound; omega. }
  unfold ltof. simpl in *. omega.
Qed.
Next Obligation.
  assert (l < (l + r) / 2).
  { apply Nat.div_le_lower_bound; omega. }
  simpl in *. omega.
Qed.
Next Obligation.
  replace n with H; eauto.
  assert (n <= H) by eapply H3, e, le_n.
  assert (H <= n) by eapply l0, H4, le_n.
  omega.
  Unshelve.
  - destruct (le_gt_dec H ((l + r) / 2)) as [ Hle | ]; simpl in *; try omega.
    generalize H1. erewrite (e _ _ Hle). inversion 1.
  - eauto.
Qed.
Next Obligation.
  exists H. repeat split; eauto.
  destruct (le_gt_dec H ((l + r) / 2)) as [ Hle | ]; simpl in *; try omega.
  generalize H1. erewrite (e _ _ Hle). inversion 1.
Qed.
Next Obligation.
  assert ((l + r) / 2 < r).
  { apply Nat.div_lt_upper_bound; omega. }
  unfold ltof. simpl in *. omega.
Qed.
Next Obligation.
  assert ((l + r) / 2 < r).
  { apply Nat.div_lt_upper_bound; omega. }
  unfold ltof. simpl in *. omega.
Qed.
Next Obligation.
  replace n with H; eauto.
  assert (H <= n) by eapply l0, H4, le_n.
  assert (n <= H) by eapply H3, e, le_n.
  omega.
  Unshelve.
  - eauto.
  - generalize (l0 _ _ H1). omega.
Qed.
Next Obligation.
  exists H. repeat split; eauto.
Qed.

停止性を保証するために付けた余計な引数を補ってやれば完成です．

Definition lower_bound l r p :=
  lower_bound' l r p (well_founded_ltof _ _ _).

型は以下の通りです．

forall l r (p : forall n, l < n <= r -> bool),
  (exists m,
    l < m <= r /\
    (forall n H, p n H = true -> m <= n) /\
    (forall n H, m <= n -> p n H = true)) ->
  { m |
    l < m <= r /\
    (forall n H, p n H = true -> m <= n) /\
    (forall n H, m <= n -> p n H = true) }

pとそれに関する性質をSectionとVariableで再帰の外に括り出せなかったのは， pの型が引数l，rに依存しているためです．

CoqにおいてProp型を持つ型の値（つまり証明）の中身を調べてSet型を持つ型の値（いわゆる，狭い意味での値）を作るのは許されていませんから，正当性のステートメントに相当する部分がややこしくなっていますね．まぁでも今は亡きanarchy proofで見た問題に似てて，ちょっといい感じです．

使用例

この実装からOCamlのコードを抽出すると，

Extract Inductive bool => "bool" ["true" "false"].
Extract Inductive sumbool => "bool" ["true" "false"].
Extract Inductive nat => int ["0" "succ"] "(fun fO fS n -> if n = 0 then fO () else fS (n-1))".

Extract Constant plus => "( + )".
Extract Constant Nat.div => "( / )".
Extract Constant le_gt_dec => "( <= )".
Extraction "bsearch.ml" upper_bound.

以下のようになります．

type 'a sig0 =
  'a
  (* singleton inductive, whose constructor was exist *)



(** val add : int -> int -> int **)

let rec add = ( + )

module Nat =
 struct
  (** val divmod : int -> int -> int -> int -> (int, int) prod **)

  let rec divmod x y q u =
    (fun fO fS n -> if n = 0 then fO () else fS (n-1))
      (fun _ -> Pair (q,
      u))
      (fun x' ->
      (fun fO fS n -> if n = 0 then fO () else fS (n-1))
        (fun _ ->
        divmod x' y (succ q) y)
        (fun u' ->
        divmod x' y q u')
        u)
      x

  (** val div : int -> int -> int **)

  let div = ( / )
 end

(** val le_gt_dec : int -> int -> bool **)

let le_gt_dec = ( <= )

(** val lower_bound' : int -> int -> (int -> __ -> bool) -> int **)

let rec lower_bound' l r p =
  if le_gt_dec r (succ l)
  then r
  else let m = Nat.div (add l r) (succ (succ 0)) in
       if p m __
       then lower_bound' l m (fun n _ -> p n __)
       else lower_bound' m r (fun n _ -> p n __)

(** val lower_bound : int -> int -> (int -> __ -> bool) -> int **)

let lower_bound l r p =
  lower_bound' l r p

Coqがあんまり賢くないのでゴミが残っていますが，例えば $\sqrt{4}$ を求めるのに使うならこんな感じです．

# lower_bound 0 10 (fun i _ -> 4 <= i * i);;
- : int = 2

提出するコードにゴミを混ぜたくないので，手動でインライン展開とかデッドコードの削除とか色々やると

let rec lower_bound l r p =
  if r <= 1 + l
  then r
  else let m = (l + r) / 2 in
       if p m
       then lower_bound l m p
       else lower_bound m r p

が得られます．こんな単純なコードを検証するために，随分と苦労させられたものです．

かわいい ↩

fetburner.core

コアダンプ

もっと二分探索に証明を付ける

以前の実装

問題点

解決法

新しい実装

使用例