実装

二分探索を実装しようとした際にまず悩むのが、どのような使い方を想定するかでしょう。

二分探索はソート済の配列から要素を探すアルゴリズムとして解説される事が多いと思いますが、競技プログラミングではむしろ、ある条件を満たす値の最小値や最大値を求めるために使われる事が多いと思います。

後者を求める関数が得られれば前者も求められるので、ここでは後者を求める関数を実装します。

Require Import Arith Div2 Omega Recdef.

Function bsearch (p : nat -> bool) n { wf lt n } :=
  match n with
  | O => 0
  | _ =>
      let m := div2 n in
      if p m then bsearch p m
      else S m + bsearch (fun x => p (S m + x)) (n - S m)
  end.
Proof.
  + intros. apply lt_div2. omega.
  + intros. omega.
  + apply lt_wf.
Defined.

一般的な二分探索の実装では[a, b)に属し、述語p(x)を満たす整数xのうち最小のものを求めるものが多いと思いますが、停止性を証明する都合上[0, n]から探すようにしています。

まぁ述語pを工夫すれば同じ事はできるので良いかと。

正当性の証明

二分探索で配列から目的の要素を探したい場合、正しい結果を得るためには配列がソートされている必要がありましたよね？

これを今の実装に当てはめるとすれば、与えられた述語p(x)はあるしきい値n0まで偽となり、それ以上のxについて真になる事でしょうか。

それに加えてしきい値が探索範囲に含まれていれば、先ほどの二分探索を行う関数はそのしきい値を返すはずです。

これをCoqで証明すると以下のようになります。

Lemma bsearch_correct : forall n p n0,
  (forall n, n0 <= n -> p n = true) ->
  (forall n, p n = true -> n0 <= n) ->
  n0 <= n ->
  bsearch p n = n0.
Proof.
  intros n.
  induction n as [[| n'] IHn] using lt_wf_ind;
    intros ? ? H H' ?;
    rewrite bsearch_equation.
  + omega.
  + remember (p (div2 (S n'))) as b.
    destruct b.
    - apply IHn; eauto.
      apply lt_div2.
      omega.
    - destruct (le_dec n0 (div2 (S n'))) as [ Hle |].
      * apply H in Hle.
        congruence.
      * { rewrite IHn with (n0 := n0 - S (div2 (S n'))); try omega.
          + intros.
            apply H.
            omega.
          + intros ? Hp.
            apply H' in Hp.
            omega.
        }
Qed.

実装がバグっていたら停止性の証明や正当性の証明で行き詰まることでしょう。

一般に、Coqで証明ができない場合は証明しようとしている命題が本当に成り立つのか疑うのが良いとされています。 Isabelleなんかは証明を始める際にQuickCheckを走らせて、反例があれば注意してくれるのでちょっと羨ましい。

使用例

Coqでプログラムを書いて証明を付けるだけでも十分面白いのですが、実際に使ってみるともっと面白そうです。

CoqにはExtractionなる機能があり、Coqで書いたプログラムと等価なOCaml等のコードを出力できるので実際にやってみましょう。

以下のように書けばOCamlのコードが得られます。

Extract Inductive bool => "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 minus => "( - )".
Extract Constant div2 => "(fun x -> x / 2)".
Extraction "bsearch.ml" bsearch.

得られたOCamlのコードは次の通りです。

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

let rec plus = ( + )

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

let rec minus = ( - )

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

let rec div2 = (fun x -> x / 2)

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

let rec bsearch p n =
  (fun fO fS n -> if n = 0 then fO () else fS (n-1))
    (fun _ ->
    0)
    (fun n0 ->
    if p (div2 (succ n0))
    then bsearch p (div2 (succ n0))
    else plus (succ (div2 (succ n0)))
           (bsearch (fun x -> p (plus (succ (div2 (succ n0))) x))
             (minus (succ n0) (succ (div2 (succ n0))))))
    n

試しに使ってみるとこんな感じ。

# (* ソート済の配列から810を探す *)
  let a = [| 114; 514; 810 |];;
val a : int array = [|114; 514; 810|]
# bsearch (fun i -> 810 <= a.(i)) (Array.length a - 1);;
- : int = 2

ARC50のB問題もこれで解いてみました。 arc050.contest.atcoder.jp

今回書いたCoqのソースコード全文が見たければこちらからどうぞ。 github.com

fetburner.core

コアダンプ

二分探索に証明を付ける

実装

正当性の証明

使用例