Think Stitch
PRINCIPIA  最近の更新


Isabelle によるヒープソート siftup の部分正当性証明

Theorem Prover Advent Calendar 2013の参加エントリです.

ヒープソートの一部である siftup という手続きの部分正当性証明をやってみたいと思います.今回も定理証明支援系 Isabelle を使います.細かいところは説明しませんが,話の流れだけお楽しみいただければ幸いです.

ヒープソートについて

ヒープソートのアルゴリズムについてはご存知あると仮定します. ヒープソートは2つの手続きから構成されています.呼び名は様々でしょうけど,ここでは siftup と siftdown と呼ぶことにします.Jon Bentley さんの本 "Programming Pearls" に合わせました.邦訳「珠玉のプログラミング」は手に入らない状態だったのですが,うれしいことに来年再出版されるそうです.

配列のインデックスは1から始まるとします.C 言語や Isabelle でもそうですが,配列のインデックスが 0 から始まる言語では 0 番目のスロットを未使用とします.インデックスを -1 して無駄なく使うこともできますけど,混乱するので,ここでは 0 を使わない方針で行きます.こういったトリックは往々にして些細なバグの原因になりますけど,必要があれば以下の証明を修正して,確認してください :P

ヒープであることの条件

配列 v[i] (1 <= i <= n) がヒープになっている条件は,おおまかにいうと,すべてのインデックス i について,2つの子供 v[2 * i] と v[2 * i + 1] がどちらも v[i] よりも大きい(逆順に並べる場合は小さい)ということです.ただし各インデックスが範囲内にあるときだけです.実際には,次のように親から子を見るよりも,子から親を見る方が条件を完結に表現できます.

\[ \forall i. 1 < i \land i \le n \Rightarrow v[i / 2] \le v[i] \]

これが,配列 v[i] がヒープを形成していることの条件です.ただし i / 2 は整数の除算(切り捨て)を表すものとします.

siftup

配列 v[i] は 1 <= i <= k (< n - 1) の範囲ですでにヒープになっているとします.範囲を1つ広くして k + 1 まで考えると,一般にはヒープになりません.手続き siftup は v[k + 1] を適切な位置に"押し上げる" ことで, 1 <= i <= k + 1 の範囲がヒープになるように調整します.Isabelle のコードを示します.

theorem "VARS (v::int list) i
  {k + 1 < length v &
   (ALL j. 1 < j & j <= k --> v!(j div 2) <= v!j)}
  i := k + 1;
  WHILE i ~= 1 & v!(i div 2) > v!i
  INV {ループ不変条件}
  DO
    v := v[i div 2 := v!i, i := v!(i div 2)];
    i := i div 2
  OD
  {ALL j. 1 < j & j <= k + 1 --> v!(j div 2) <= v!j}"

v は整数のリストです.ここでは配列として使います.配列の要素 v[i] を Isabelle では v ! i と書きます.div は整数の除算です.今回も簡単のために順列の条件は外してあります.

配列への代入 v[i] := x は v := v[i:=x] と書きます.インデックス i と j の内容を交換する場合は

v := v[i := v!j, j := v!i]

と書くことができます.ループ本体では親と子を入れ替えているわけです.

アルゴリズムは単純で,i = k + 1 からはじめて,親の方が大きい間,要素を上に向かって"押し上げて"行きます.

事前条件は 1 <= i <= k の範囲がヒープになっていることです.事後条件は1つ範囲を広げて k + 1 までがヒープになっていることです.

ループ不変条件

インデックス k + 1 にあった値をインデックス i まで押し上げた途中の状態(計算の断面)を考えると,v[i] 以外のところはヒープになっているはずです.そこでループ不変条件の候補として,次の条件を採用することができます.

\[ \forall j. 1 < j \land j \le n \land j \neq i \Rightarrow v[j / 2] \le v[j] \]

このループ不変条件は事後条件を弱めた形になっています.j = i の場合の積項が1つ除かれているからです.ループ不変条件が事後条件を弱めたものになるという現象はよく起こります.これを逆に考えると,ループ不変条件の候補を事後条件から導きだすことができるということです.

おおまかにいうと,ループを脱出した時点ではループ条件の否定 v[i/2] <= v[i] が成り立ちますから,これとループ不変条件を合わせると事後条件が成り立つことがわかります.ループ条件,ループ不変条件,そして(ループ脱出直後の)事後条件の間にはこういう関係が成り立つわけです.

最初のトライ

コマンド apply(vcg) を使って証明課題に変換し,あとは Isabelle さんにがんばってもらいます.

theorem "VARS (v::int list) i
  {k + 1 < length v & (ALL j. 1 < j & j <= k --> v!(j div 2) <= v!j)}
  i := k + 1;
  WHILE i ~= 1 & v!(i div 2) > v!i
  INV {
    k + 1 < length v & 1 <= i & i <= k + 1 &
    (ALL j. 1 < j & j <= k + 1 & j ~= i --> v!(j div 2) <= v!j)
  }
  DO
    v := v[i div 2 := v!i, i := v!(i div 2)];
    i := i div 2
  OD
  {ALL j. 1 < j & j <= k + 1 --> v!(j div 2) <= v!j}"

apply(vcg)
apply(auto simp add: nth_list_update)

ここまでで,次の2つのゴールが残りました.

proof (prove): step 2

goal (2 subgoals):
 1. !!v j.
       [| Suc k < length v; Suc 0 <= j div 2;
          j div 2 <= Suc k;
          ALL ja.
             Suc 0 < ja &
             ja <= Suc k & ja ~= j div 2 -->
             v ! (ja div 2) <= v ! ja;
          j div 2 ~= Suc 0;
          v ! (j div 2) < v ! (j div 2 div 2);
          j div 2 ~= j; Suc 0 < j; j <= Suc k;
          j ~= j div 2 div 2 |]
       ==> v ! (j div 2 div 2) <= v ! j
 2. !!v i j.
       [| Suc k < length v; Suc 0 <= i; i <= Suc k;
          ALL j.
             Suc 0 < j & j <= Suc k & j ~= i -->
             v ! (j div 2) <= v ! j;
          i ~= Suc 0; v ! i < v ! (j div 2);
          i ~= j div 2; i ~= j; Suc 0 < j; j <= Suc k;
          j ~= j div 2; i div 2 = j div 2 |]
       ==> v ! i <= v ! j

先に2番目のゴールを見ると, i div 2 = j div 2 でかつ i ~= j (~ は否定)ですから,i と j は i div 2 の異なる2人の子供ということです.加えて v ! i < v ! (j div 2) という条件がありますから,これは証明できるはずです.実際以下のようにして証明できます.defer は先頭のゴールを後回しにするというコマンドです.

defer
apply(drule_tac x="j" in spec)
apply(drule mp)
apply(clarsimp)
apply(clarsimp)

結果,ゴールはあと1つになりました.

proof (prove): step 7

goal (1 subgoal):
 1. !!v j.
       [| Suc k < length v; Suc 0 <= j div 2;
          j div 2 <= Suc k;
          ALL ja.
             Suc 0 < ja & ja <= Suc k & ja ~= j div 2 -->
             v ! (ja div 2) <= v ! ja;
          j div 2 ~= Suc 0;
          v ! (j div 2) < v ! (j div 2 div 2);
          j div 2 ~= j; Suc 0 < j; j <= Suc k;
          j ~= j div 2 div 2 |]
       ==> v ! (j div 2 div 2) <= v ! j

文字が i から j に変わってしまっていますが,押し上げている要素の親と子の間の,1つ飛ばした大小関係を示せといわれています.しかしこれは無理です.j と j div 2 の間の大小関係についてわかっていない上に,v ! (j div 2) < v ! (j div 2 div 2) という条件がありますからますます無理です.

では成り立たないのかというとそんなことはありません.なぜかということ,要素を押し上げたとき,代わりに下がった要素とその親との間には大小関係があったからです.ヒープだったからです.したがって押し上げた後でも,1つ飛ばした要素の間に大小関係はあるのです.しかしループ不変条件にはそのことが入っていないので証明することができません.

Programming Pearls の中で,この点について脚注があります.本文では上のプログラムと同じループ不変条件を使っているのですが,これだけではゴールを証明できないから条件を強めるべきだと Knuth 先生からつっこみが入ったという話です.

Bentley さんは条件に気づいていなかったのか,それとも話を簡単にするためにあえて省略したのか,どちらなのかはわかりませんが,プログラムを書いているとき,暗に仮定したことを使ってしまってバグを出すということはよくありますよね(数学といっしょです).バグにならなくても,一部しか考えていなかった,考え漏れがあったのにたまたま運がよかっただけと後で気づくこともあります(気づかないままのことももちろんあります).しかし証明をしてみると条件が足りないので必ず気がつきます.紙の上で証明を考えていると間違う可能性がありますけど,ツールを使っていれば許してもらえない(?)ので必ず気がつきます.それと同時に,なぜそのプログラムが正しいのか,どういうメカニズムで望む結果が得られるのかという理由がとことんわかります.これが証明を行うもう1つの意義だと思います.

ループ不変条件の強化

さて,話を戻して,足りない条件をループ不変条件に追加します.このようにループ不変条件を強化することもよくあることです.証明を試みて,証明できないで残ったゴールを見ると,プログラムの性質に関する知見と,強化すべき条件のヒントが得られます.ここでは直接的に次の条件を追加すればいいでしょう.

1 <= i div 2 & 2 * i <= k + 1 --> v!(i div 2) <= v!(2 * i))
1 <= i div 2 & 2 * i + 1 <= k + 1 --> v!(i div 2) <= v!(2 * i + 1))

押し上げている要素をはさんだ親と子の間に大小関係があるという条件です.これらを追加して再挑戦します.

正当性証明

強化したループ不変条件を使うと,以下のように証明に成功します.2つ補助定理を用意して使いました. 証明の過程を追うと,さらに認識を深めることができるのですが,長いですし Isabelle の知識が必要なのでやめておきます.重要なことはツールで証明がチェックされているので,アルゴリズムは仕様に対して正しいということです.

lemma x0: "
  [| k + 1 < length v; 1 <= i; i <= k + 1;
     ALL j. 1 < j & j <= k + 1 & j ~= i --> v!(j div 2) <= v!j;
     1 < m; m <= k + 1; m ~= i
  |]
  ==> v ! (m div 2) <= v ! m"
apply(auto)
done

lemma x1: "ALL k::nat. (EX m. k = 2 * m | k = 2 * m + 1)"
apply(auto)
apply(arith)
done

theorem "VARS (v::int list) i
  {k + 1 < length v & (ALL j. 1 < j & j <= k --> v!(j div 2) <= v!j)}
  i := k + 1;
  WHILE i ~= 1 & v!(i div 2) > v!i
  INV {
    k + 1 < length v & 1 <= i & i <= k + 1 &
    (ALL j. 1 < j & j <= k + 1 & j ~= i --> v!(j div 2) <= v!j) &
    (1 <= i div 2 & 2 * i <= k + 1 --> v!(i div 2) <= v!(2 * i)) &
    (1 <= i div 2 & 2 * i + 1 <= k + 1 --> v!(i div 2) <= v!(2 * i + 1))}
  DO
    v := v[i div 2 := v!i, i := v!(i div 2)];
    i := i div 2
  OD
  {ALL j. 1 < j & j <= k + 1 --> v!(j div 2) <= v!j}"
apply(vcg)
apply(clarify)
apply(simp)
defer
apply(clarsimp)
apply(case_tac "i=j")
apply(clarsimp)
apply(clarsimp)
(* loop inv*)
apply(clarsimp)
(* range k *)
apply(rule conjI)
apply(arith)
apply(rule conjI)
apply(arith)
(* loop inv main *)
apply(rule conjI)
apply(clarsimp)
(*
 j ~= i div 2 ==>
    v[i div 2 := v ! i, i := v ! (i div 2)] ! (j div 2)
 <= v[i div 2 := v ! i, i := v ! (i div 2)] ! j
*)
apply(simp add: nth_list_update)
apply(rule conjI)
apply(clarsimp)
apply(rule conjI)
apply(clarsimp)
apply(clarsimp)
apply(rule conjI)
apply(clarsimp)
apply(clarsimp)
(*
  v ! i < v ! (j div 2)
  i div 2 = j div 2;
  i ~= j
  ==> v ! i <= v ! j
*)
apply(blast intro: x0 less_imp_le less_le_trans)
apply(clarsimp)
apply(rule conjI)
apply(clarsimp)
(*
  v ! (j div 2 div 2) <= v ! j
*)
apply(insert x1)
apply(rotate_tac -1)
apply(drule_tac x="j" in spec)
apply(clarsimp)
apply(erule disjE)
apply(drule mp)
apply(clarsimp)
apply(arith)
apply(simp)
apply(clarsimp)
apply(clarsimp)
(* loop inv sub 1 *)
apply(rule conjI)
apply(clarsimp)
apply(simp add: nth_list_update)
apply(clarsimp)
apply(frule_tac x="2 * (i div 2)" in spec)
apply(drule_tac x="i div 2" in spec)
apply(clarsimp)
apply(rotate_tac -1)
apply(drule mp)
apply(arith)
apply(rotate_tac -2)
apply(drule mp)
apply(arith)
apply(erule order_trans)
apply(clarsimp)
(* loop inv sub 2 *)
apply(clarsimp)
apply(simp add: nth_list_update)
apply(clarsimp)
apply(frule_tac x="Suc (2 * (i div 2))" in spec)
apply(drule_tac x="i div 2" in spec)
apply(clarsimp)
apply(rotate_tac -1)
apply(drule mp)
apply(arith)
apply(erule order_trans)
apply(clarsimp)
done

補足

こんなかんたんなコードでも,ずいぶんと奥が深いものだと思います.あらためてプログラミングは難しいなと思います.

ここまで厳密に考えなくても,動いているから大丈夫,ということにはならないと思います.しかし証明をするのはたいへんで,時間もかかるので,トレードオフとしてテストだけにするということはあるでしょう.場合分けを考慮したテストならば,かなり完全に近い,満足のいく検証ができるでしょう.ただ,そのような網羅性の高いテストを作るためには,プログラムの性質をよく知らなければならないという,にわとりとたまごみたいな問題がありますけど...仕様(事前・事後条件)からテストケースを生成してくれるようなプログラムがあればいいんですけどね.

ヒープソートを構成するもう1つの手続き siftdown についても同じように証明することができます.難しさは同じ程度なんですが,場合分けがあるのでちょっと面倒になります.

2013/12/11
© 2013,2014,2015 PRINCIPIA Limited