Think Stitch
PRINCIPIA  最近の更新


ループ不変条件を"解いて"求める

以前に,ループ不変条件を事後条件から求めるという話をしました。 今日は,ループのコードを先に思いついたときに,後からループ不変条件を求める方法について書いてみたいと思います。

ループのコードが先に得られたとしても,ちょっと複雑なコードだと正しさに不安が残ります。 ループの理解が不十分だと,テストケースの抽出にも不安が残ります。 そのようなときはやはりループ不変条件を求めて正しさを証明したくなります。 そうすれば正しいという確信だけでなく,証明の過程でコードに対する深い理解も得られるからです。

ループ不変条件そのものについては,以下の過去記事を参考にしていただければ幸いです。

ループ不変条件を未知関数とする方程式を求める

基本的なアイデアは,ループ不変条件を未知関数とする方程式を立てるということです。 百聞は一見にしかず,まずはベキ乗を計算する簡単な例を見ていただくことにします。

以下に示すコードは,与えられた自然数 x, y に対してベキ乗 x^y を求めるプログラムです。 結果を格納する変数 z を用意して,x を y 回だけ掛け合わせるという単純なコードです。 求める結果を仕様変数 C で表しています。 定数 C があること,および事前条件と事後条件の形からエネルギー保存型のループ不変条件が予想されます。 そこで,エネルギー関数を E x y z と置いてみることにします。 E は未知の関数です。

lemma "VARS (z::nat) y
  {x ^ y = C}
  z := 1;
  WHILE y ~= 0
  INV {E x y z = C}
  DO
    z := z * x;
    y := y - 1
  OD
  {z = C}"

apply(vcg)
apply(auto)

コマンド vcg を使って証明課題に変換し,auto を使って整理すると以下の3つのゴールが残ります。

goal (3 subgoals):
 1. !!y. E x y (Suc 0) = x ^ y
 2. !!z y. 0 < y ==> E x (y - Suc 0) (z * x) = E x y z
 3. !!z. z = E x 0 z

ごらんのとおり,これは未知関数 E に対する方程式です。 これを解くことができれば E を求めることができます。 この場合は以下のようにして簡単に求めることができます。 (Suc は後者関数です。つまり Suc 0 = 1 です。)

\[ \begin{split} E\ x\ y\ z &= E\ x\ (y-1)\ z\cdot x \\ &= E\ x\ (y-2)\ z\cdot x^2 \\ &= E\ x\ (y-3)\ z\cdot x^3 \\ &\cdots \\ &= E\ x\ 1\ z\cdot x^{y-1} \\ &= E\ x\ 0\ z\cdot x^y \\ &= z\cdot x^y \end{split} \]

条件2を繰り返し使って,最後に条件3を使いました。結果は条件1を満たしていることがわかります。

この結果を使うと,ループの正当性を証明することができます。

lemma "VARS (z::nat) y
  {x ^ y = C}
  z := 1;
  WHILE y ~= 0
  INV {z * x ^ y = C}
  DO
    z := z * x;
    y := y - 1
  OD
  {z = C}"

apply(vcg)
apply(auto)
apply(erule notE)
apply(case_tac y)
apply(auto)
done

ベキ乗 その2

もう少し凝ったループでやってみます。対数時間でベキ乗を計算するプログラムです。

lemma "VARS (z::nat) x y
  {x ^ y = C}
  z := 1;
  WHILE y ~= 0
  INV {E x y z = C}
  DO
    IF y mod 2 = 0 THEN
      x := x * x;
      y := y div 2
    ELSE
      z := z * x;
      y := y - 1
    FI
  OD
  {z = C}"
apply(vcg)
apply(auto)

先ほどと同じようにエネルギー保存型のループ不変条件を仮定し,エネルギー関数を未知関数 E とします。 以下の条件が得られます。

 1. !!x y. E x y (Suc 0) = x ^ y
 2. !!z x q. 0 < q ==> E (x * x) q z = E x (2 * q) z
 3. !!z x y.
       [| 0 < y; y mod 2 = Suc 0 |]
       ==> E x (y - Suc 0) (z * x) = E x y z
 4. !!z x. z = E x 0 z

今度はそのまま解くことができませんでした。 しかし具体的な数で実験してみると予想がつきます。

\[ \begin{split} E\ x\ 10\ z &= E\ x^2\ 5\ z \\ &= E\ x^2\ 4\ z\cdot x^2 \\ &= E\ x^4\ 2\ z\cdot x^2 \\ &= E\ x^8\ 1\ z\cdot x^2 \\ &= E\ x^8\ 0\ z\cdot x^{10} \\ &= z\cdot x^{10} \end{split} \]

答えは先ほどと同じになりそうです。 実際,以下のようにして証明することができます。

lemma "VARS (z::nat) x y
  {x ^ y = C}
  z := 1;
  WHILE y ~= 0
  INV {z * x^y = C}
  DO
    IF y mod 2 = 0 THEN
      x := x * x;
      y := y div 2
    ELSE
      z := z * x;
      y := y - 1
    FI
  OD
  {z = C}"
apply(vcg)
apply(auto)
apply(metis comm_semiring_1_class.normalizing_semiring_rules(31) comm_semiring_1_class.normalizing_semiring_rules(7) power2_eq_square)
apply(metis Suc_pred power_Suc)
done

このように,直接解くことはできないにしても,ヒントを与えてくれます。

木の平坦化

もっと複雑な例として,木の平坦化してリストに変換するループを考えてみます。 ここでの目的は,やや作為的ですが再帰を使わずに自分でスタック管理をして平坦化することだとします。 まずデータ構造と,それから仕様を記述するための関数的な定義を用意します。

datatype 'a tree = Nil | Node 'a "'a tree" "'a tree"

fun flatten :: "'a tree => 'a list"
  where
  "flatten Nil = []" |
  "flatten (Node x left right)
      = (flatten left) @ (x # (flatten right))"

ここで # は cons,@ は append です。 以下のようにして計算できます。

value "flatten (Node (3::int)
                     (Node 1
                           Nil
                           (Node 2 Nil Nil))
                     (Node 5
                           (Node 4 Nil Nil)
                           Nil))"

=> "[1, 2, 3, 4, 5]" :: "int list"

ノードの構成要素を取り出すための補助的な関数を3つ定義しておきます。

fun node_val :: "'a tree => 'a"
  where
  "node_val (Node x left right) = x"

fun left :: "'a tree => 'a tree"
  where
  "left (Node x l r) = l"

fun right :: "'a tree => 'a tree"
  where
  "right (Node x l r) = r"

コードを次のように書きました。 変数 p に平坦化すべき木が入っているとします。 これを仕様変数 T で表します。 変数 r は途中までできたリスト,変数 s はおおまかにいうと処理を待っているノードのリストです。 つまりスタックです。

平坦化の戦略は次のとおりです。まず右の枝を Nil に達するまでたどっていきます。 右へ進むのはリストが後ろからできるからです。 たどる過程で,ノードについている値と左の枝をペアにしてスタックに積んでいきます。 Nil に達したら親のノードが右端ということですから,スタックから親を取り出して値を r に連結します。 続いて左の枝を平坦化します。以上の処理をスタックが空で,かつ p が Nil になるまで繰り返します。

lemma
  "VARS p r s x
  {p = T}
  r := [];
  s := [];
  WHILE p ~= Nil | s ~= []
  INV {E p r s = flatten T}
  DO
    IF p = Nil THEN
      x := fst (hd s);
      p := snd (hd s);
      s := tl s;
      r := x # r
    ELSE
      s := (node_val p, left p) # s;
      p := right p
    FI
  OD
  {r = flatten T}"
apply(vcg)
apply(auto)

事後条件の形から,flatten T を定数とするエネルギー関数を仮定しました。 証明課題を整理すると以下の5つの条件が出てきます。 このうち主要なものは 2 と 3 です。 (2 と 4 は同じ条件です。)

goal (5 subgoals):
 1. E T [] [] = flatten T
 2. !!p r s.
       [| E p r s = flatten T; p ~= Nil |]
       ==> E (right p) r ((node_val p, left p) # s) = flatten T
 3. !!r s.
       [| E Nil r s = flatten T; s ~= [] |]
       ==> E (snd (hd s)) (fst (hd s) # r) (tl s) = flatten T
 4. !!p r s.
       [| E p r s = flatten T; s ~= []; p ~= Nil |]
       ==> E (right p) r ((node_val p, left p) # s) = flatten T
 5. !!r. E Nil r [] = flatten T ==> r = flatten T

条件 2 では p が Nil ではないので,ノードに置き換えて整理します。 条件 3 では s が空リスト [] ではないので,ペアに置き換えて整理します。 これは以下のコマンドで実行できます。

prefer 2
apply(case_tac p)
apply(clarsimp)
apply(clarsimp)
prefer 3
apply(case_tac s)
apply(clarsimp)
apply(clarsimp)

その結果,2つの条件は以下のようになりました。 ループの動きと対応していることがわかると思います。

 1. !!r a b list.
       E Nil r ((a, b) # list) = flatten T
       ==> E b (a # r) list = flatten T
 2. !!r s a tree1 tree2.
       E (Node a tree1 tree2) r s = flatten T
       ==> E tree2 r ((a, tree1) # s) = flatten T

これを解くことができるかといわれれば,正直できません。 しかし条件が整理されたおかげで,議論を進めやすくなったとは思います。 ベキ乗の例のように,具体的な木で計算してみると,E 自身 flatten と同じ働きを持つことがわかります。 処理を保留したノードをスタック s に積んでいるので,それらすべてを平坦化したものが最終的なリストの一部になります。 実際,エネルギー関数を構成する要素は以下の3つです。

  1. s に含まれているノード(逆順)
  2. p が指しているノード
  3. 途中までできたリスト r

これらをこの順序で連結すれば flatten T になります。 それを上の2つの条件に当てはめてみれば,たしかに成り立つことが確かめられます。

最終的な結果は次のようになります。 スタックに格納されているノードをすべて平坦化するために関数 map を使いました。 % は λ です。 さらにその結果をすべて連結するために,補助関数 lf を定義して使いました。 ちょっと複雑ですけど,3つの要素が連結されていることは見て取れると思います。

fun lf :: "'a list list => 'a list"
 where
  "lf [] = []" |
  "lf (x#xs) = x @ (lf xs)"

value "lf [[1::nat,2,3],[4,5],[],[6],[7,8]]"
=> "[1, 2, 3, 4, 5, 6, 7, 8]" :: "nat list"

lemma x0: "lf (xs @ ys) = (lf xs) @ (lf ys)"
apply(induct_tac xs)
apply(auto)
done

lemma
  "VARS p r s x
  {p = T}
  r := [];
  s := [];
  WHILE p ~= Nil | s ~= []
  INV {lf (map (%(x,l). flatten l @ [x]) (rev s)) @ (flatten p) @ r = flatten T}
  DO
    IF p = Nil THEN
      x := fst (hd s);
      p := snd (hd s);
      s := tl s;
      r := x # r
    ELSE
      s := (node_val p, left p) # s;
      p := right p
    FI
  OD
  {r = flatten T}"
apply(vcg)
apply(auto)
apply(case_tac p)
apply(auto)
apply(unfold x0)
apply(auto)
apply(case_tac s)
apply(auto)
apply(unfold x0)
apply(auto)
apply(case_tac p)
apply(auto)
done

手品はなしです。機械的に求まるわけではありません。しかし,ループ不変条件を思いつくためのヒントになりますし,確認すべき条件を整理するための手続きとして使うことができます。 それになにより,コードの性質を深く理解するために使うことができると思います。

2014/02/01

このエントリーをはてなブックマークに追加

© 2013,2014,2015 PRINCIPIA Limited