Just $ A sandbox

プログラミングと計算機科学とかわいさ

無限ループによる証明

Software FoundationsのIMPあたりをやってて疑問が湧いた.

w = while True do Skip endが停止しないことを示したい.
つまり, $\langle w, s \rangle \to s' \Rightarrow False$ を示したい. ($s, s' : \mathbb{Z} \to \mathbb{Z}$ は変数の値を保持する環境)

普通に証明するなら以下のようになると思う:
証明図の高さに関する帰納法による. まず高さが0のとき, 明らかに<w,s>→s'は公理から導くことはできないので高さ0の証明図を持たない. よって$\langle w, s \rangle \to s' \Rightarrow False$である.
次に, $\forall s, s'. \langle w, s \rangle \to s'$が高さnの証明図をもたないと仮定し, $\langle w, s \rangle \to s'$が高さ(n+1)の証明図をもたないことを証明する. $\langle w, s \rangle \to s'$が高さ(n+1)の証明図をもつことを仮定すると, wの定義とwhile-loopに関する推論規則より$\langle w, s \rangle \to s'$が高さnの証明図で示せる. これは帰納法の仮定に矛盾するので証明が完了する.

Agdaなどの依存型を持つ言語では, prop :: w -> Falseに対応する函数を作ってやればよい.
<w,s>→s'にあたるデータ型を再帰的に定義している場合, 例えばこの証明は無限ループを使ってprop (While-Loop True _ k) = kみたいにすると, 証明図の仮定に当たる, 「一つ前の」<w,s>→s'を使うことができる.
けれど一体これは何を証明していることになるんだろう?
明らかに上の証明とは違う気がするし, そもそも無限ループになっているのに証明が完了しているというのがよく分からない. 直観的には「<w,s>→s'を仮定すると全く同じ形が仮定の部分に出現する」ということが重要なんだろうという気がするけれど、そこからどうやって矛盾を導いているんだろう.

ちなみにIsabelleでは以下のようにやった.
上の証明の気持ちがある程度反映されているつもり. (と言っても証明だけでは伝わらない気もするけど)

lemma while_induction:
  assumes red: "(WHILE b DO c END) %: st ⇓ st'"
  and base: "⋀t t'. ⟦ beval3 t b = False; t' = t ⟧ ⟹ P t t'"
  and step: "⋀t t' t''. ⟦ beval3 t b = True; c %: t ⇓ t'; (WHILE b DO c END) %: t' ⇓ t''; P t' t'' ⟧ ⟹ P t t''"
  shows "P st st'"
proof-
  have "(WHILE b DO c END) %: st ⇓ st' ⟹ P st st'"
    using ceval.induct [of _ _ _ "λc0 t t'. WHILE b DO c END = c0 ⟶ P t t'"] apply rule
    apply (simp, auto, simp add: base)
    apply (fastforce simp add: step)
    apply (fastforce simp add: step)
    done
  thus ?thesis by (simp add: red)
qed
 
theorem loop_never_stops: "¬ (loop %: st ⇓ st')"
proof (unfold loop_def, auto)
  assume p: "WHILE BTrue DO SKIP END %: st ⇓ st'"
  show "False"
    using while_induction [of BTrue SKIP st st' "λ_. λ_. False"]
    by (rule, simp add: p, auto)
qed

de Bruijn Indexとβ変換

$ KMN \equiv (\lambda xy. x)MN \longrightarrow_{\beta} M $を示したいとする.

α変換

$y \notin \text{FV}(M)$ならば, $ (\lambda xy. x)MN \to (\lambda y. M) N \to M[y:=N] \equiv M $となる.
問題は$y \notin \text{FV}(M)$ならば$ M[y:=N] \equiv M $の部分で, これを示すにはMについての帰納法を使う.

$M \equiv x$(変数)の時は$x=y, x \neq y$で場合分けして明らか. $M \equiv (\lambda x. M_0)$の時は, α変換によって$x \neq y$としてよく, このとき$M[y:=N] \equiv (\lambda x. M_0)[y:=N] \equiv (\lambda x. M_0[y:=N])$となり, あとは帰納法の仮定よりOK. $M \equiv (M M')$のときは帰納法の仮定と代入の性質から明らか.

さてこれをIsabelle(proof assistant)で証明しようとすると, α変換によって$x \neq y$としてとれるというところを示すのが難しい.
α変換の煩わしさを軽減するためにλ式の定義を少し変えて, de Bruijn Indexによる定義をする.

de Bruijn Index

<λ-expr> ::= <var :: nat> | λ. <λ-expr> | (<λ-expr> <λ-expr>)

によって定義する. 変数の代わりに自然数を使い, さらに函数抽象をするときに束縛変数を指定しない.
<var :: nat>は, その変数のいくつ前にあるλによって束縛されているかを表す(この変数の添字の割り振り方をde Bruijn Indexと呼ぶらしい).

例として, $\lambda. \lambda. 1 \equiv \lambda xy. x$, $\lambda. \lambda. 0 (\lambda. 0 1) \equiv \lambda xy. y (\lambda z. z y)$である.

de Bruijn Indexのλ式のβ変換をどのように定義するかを考える.
λ式Mの中で, Mの中の変数nがあったとき, nに対応するλがMに含まれていることとnがMの束縛変数であることは等しい.

このとき, $(\lambda. M) N$のβ変換は, Mの中の自由変数の値を一つ減らし(λがはずれるため), Nの中の自由変数の値はλに対応するMの中の変数の深さだけ増やす.
例: $(\lambda. \lambda. \lambda. 5 2) (\lambda. 0 1) \to (\lambda. \lambda. (5-1) (\lambda. 0 (1+2))) \equiv (\lambda. \lambda. 4 (\lambda. 0 3)) $

de Bruijn Indexによるλ式の表現は, 束縛変数の添字の振り方に自由度がないのでα変換を考える必要がない.

元の問題へ

結局元の問題に戻ると, $ KMN \longrightarrow M $を示すのが目的だった.
de Bruijn Indexでは, $ (\lambda. \lambda. 1) M N \to (\lambda M) N \to M $ となることを示せばよい.
けれど, これもあまり上手くいかなかったので結局どうすればいいのだろう.

K,M,Nが変数条件を満たす時, unshift 1 (λ shift 1 (shift 1 M)) $ N ⟶ Mが示せればいいと思うけれど.

原始帰納的函数で素数を計算

原始帰納函数; prf

原始帰納函数(primitive recursive function; prf)とは, 以下で定義される函数 $ \mathbb{N}^n \to \mathbb{N} $ のこと:

  • 定数函数 $zero : \mathbb{N}^0 \to \mathbb{N}$, 後続函数 $suc : \mathbb{N} \to \mathbb{N}$, 射影函数 $p^n_i : \mathbb{N}^n \to \mathbb{N}$ はprfである
  • prf $f : \mathbb{N}^m \to \mathbb{N}, f_i : \mathbb{N}^n \to \mathbb{N} (1 < i \leq m)$に対し, 合成 $g(\vec{x}) = f(f_1(\vec{x}), \cdots ,f_m(\vec{x}));\; g : \mathbb{N}^n \to \mathbb{N}$ はprfである
  • prf $f : \mathbb{N}^{n} \to \mathbb{N}, h : \mathbb{N}^{n+2} \to \mathbb{N}$に対し, 原始帰納法で定まる函数 $g(0,\vec{x}) = f(\vec{x});\;\; g(y+1,\vec{x})=h(g(y,\vec{x}),y,\vec{x});\;\; g : \mathbb{N}^{n+1} \to \mathbb{N}$ はprfである

n番目の素数を計算するprf $pr : \mathbb{N} \to \mathbb{N}$が定義できることを見る.

準備(自然数, n項直積)

自然数とタプル(長さ付きリスト)といくつかの函数を定義しておく.

data Nat = Z | S Nat deriving (Show)
type N0 = Z
type N1 = S N0
type N2 = S N1
type N3 = S N2
 
type family (+) a b where
  Z + m = m
  (S n) + m = S (n + m)
 
data Tuple n a where
  TEmpty :: Tuple N0 a
  TSuc :: a -> Tuple n a -> Tuple (S n) a
 
instance Functor (Tuple n) where
  fmap f TEmpty = TEmpty
  fmap f (TSuc x xs) = TSuc (f x) (fmap f xs)
 
lengthT :: Tuple n a -> Nat
lengthT TEmpty = Z
lengthT (TSuc _ t) = S (lengthT t)
 
infixr 4 %:
(%:) = TSuc

`Tuple n a`は函数合成を定義するために函数の組を与えるときに使う.
どうせ同じ型の組しか作らないのでこれで十分. 長さ付きにしたのは型チェックのため(後述).

prfの定義

上の定義をそのままGADTで書き下す.
ただし, `PrimRec n`で函数 $\mathbb{N}^n \to \mathbb{N}$を表すことにする. `PrimRec`に「引数の数」という情報を型で持たせることで, 引数の数が合わない場合はコンパイルエラーで弾くことができる.
(e.g. `suc(1,2)`や`add(2)`みたいなものはコンパイルエラーになる)

-- datatype of `p.r. function :: N^p -> N`
data PrimRec p where
  Zero :: PrimRec N0
  Suc :: PrimRec N1
  Proj :: Nat -> PrimRec n
  Comp :: PrimRec m -> Tuple m (PrimRec n) -> PrimRec n
  Induct :: PrimRec n -> PrimRec (S (S n)) -> PrimRec (S n)

infixl 6 %.
(%.) = Comp

prfの実行

さて`PrimRec`を定義したのでprfを宣言できるようになった.
函数を定義したので評価(実行)できるようにしたいということで`evalPR`を定義する.
evalPRは以下のような函数として定義したい.

\[ \begin{eqnarray}
eval(zero,*) &=& 0 \\
eval(suc,n) &=& n+1 \\
eval(p^n_i,\vec{x}) &=& x_i \\
eval(f \circ \{f_i\}_i,\vec{x}) &=& f(f_1(\vec{x}), \cdots f_m(\vec{x})) \\
eval(induction(f,h),0,\vec{x}) &=& f(\vec{x}) \\
eval(induction(f,h),y+1,\vec{x}) &=& h(eval(induction(f,h),y,\vec{x}),y,\vec{x})
\end{eqnarray} \]

ここで, Tupleのi番目を抜き出す操作が必要になるのでそれは別途`(!) :: Tuple n a -> Nat -> a`を定義しておく.

(!) :: Show a => Tuple n a -> Nat -> a
(TSuc n _) ! Z = n
(TSuc _ t) ! (S n) = t ! n
k ! a = error $ show k ++ " ! " ++ show a

evalPR :: PrimRec n -> Tuple n Nat -> Nat
evalPR Zero _ = Z
evalPR Suc (TSuc n TEmpty) = S n
evalPR (Proj k) t = t ! k
evalPR (Comp f fs) xs = evalPR f (go fs xs)
 where
   go :: Tuple m (PrimRec n) -> Tuple n Nat -> Tuple m Nat
   go (TSuc f fs) xs = TSuc (evalPR f xs) (go fs xs)
   go TEmpty xs = TEmpty
evalPR (Induct f h) (TSuc Z xs) = evalPR f xs
evalPR g@(Induct f h) (TSuc (S y) xs) = evalPR h (evalPR g (y %: xs) %: y %: xs)
evalPR _ z = case z of {}

prの構成

素数を計算する函数prを構成する.
まずは準備.

infixr 8 %..
(%..) f x = f %. (x %: TEmpty)
 
extend :: PrimRec N0 -> PrimRec n
extend n = n %. TEmpty
 
arg :: Int -> PrimRec n
arg n = Proj (nat n)
 
one :: PrimRec N0
one = Suc %.. Zero
 
two :: PrimRec N0
two = Suc %.. one

add(x,y)=x+y, pred(x)=x-1, sub(x,y)=if x < y then 0 else x-y, mult(x,y)=x*y
が原始帰納法で定義できる.
(例えば, add(x,0)=x; add(x,y+1)=suc(add(x,y)))

add :: PrimRec N2
add = Induct (arg 0) (Suc %. (arg 0 %: TEmpty))
 
predpr :: PrimRec N1
predpr = Induct Zero (Proj (nat 1))
 
sub :: PrimRec N2
sub = sub' %. (arg 1 %: arg 0 %: TEmpty) where
  sub' = Induct (arg 0) (predpr %. (arg 0 %: TEmpty))
 
mult :: PrimRec N2
mult = Induct (extend Zero) (add %. (arg 2 %: arg 0 %: TEmpty))

次に述語(Bool値をとる函数)を定義する.
今は便宜上, 原始帰納的述語を, {0,1}へのprfとして話を進めることにする.

述語$x=0$,$x=y$,$x \leq y$,$x < y$,not,and,orが定義できる.
(例えば, eq0(x)=sub(1,sub(1,x))とすれば, x=0のときeq0(x)=0, それ以外はeq0(x)=1となる)

equalZ :: PrimRec N1
equalZ = sub %. (extend one %: sub %. (extend one %: arg 0 %: TEmpty) %: TEmpty)
 
equal :: PrimRec N2
equal = equalZ %.. add %. (sub %. (arg 0 %: arg 1 %: TEmpty) %: sub %. (arg 1 %: arg 0 %: TEmpty) %: TEmpty)
 
leq :: PrimRec N2
leq = equalZ %.. sub %. (arg 0 %: arg 1 %: TEmpty)
 
lneq :: PrimRec N2
lneq = andp %. (leq %. (arg 0 %: arg 1 %: TEmpty) %: notp %.. equal %. (arg 0 %: arg 1 %: TEmpty) %: TEmpty)
 
notp :: PrimRec N1
notp = sub %. (extend one %: arg 0 %: TEmpty)
 
andp :: PrimRec N2
andp = equalZ %.. add %. (arg 0 %: arg 1 %: TEmpty)
 
orp :: PrimRec N2
orp = mult %. (arg 0 %: arg 1 %: TEmpty)

次に, $prod_r(y,\vec{x}) = \Pi_{z < y} r(z,\vec{x})$なるprodを定義する.
これを用いて, 原始帰納的述語pに対し, 原始帰納的述語$\exists z < y: p(z,\vec{x})$が定義できる.
同様にして函数$\sum_{z < y} r(z,\vec{x})$や述語$\forall z < y: p(z,\vec{x})$などが定義できる.

さて, ここで初めて可変長の引数を扱う必要が出てくる.
函数prodはprod(r,y,x1..xs)に対しr(z,x1..xs)を掛け算していく函数なので, 受け取った引数の「2番目から最後まで」(=x1..xs)を受け取り, rの0番目にz, 1番目以降にx1..xsを指定して合成する必要がある.
これをするために, $\text{varTuple} :: \text{PrimRec(n)} \to \mathbb{N}^n; \;\; \text{varTuple} _ = (p^n_0, \cdots p^n_{n-1})$なる函数を定義してやれば良さそうに見えるのでこれを定義する.

(varTupleの定義のために手間のかかるハックをしているけれど, もっと綺麗に書けない?)

class VarLenTuple n where
  varLenTuple :: Tuple n Nat
 
instance VarLenTuple Z where
  varLenTuple = TEmpty
 
instance (VarLenTuple n) => VarLenTuple (S n) where
  varLenTuple = Z %: varLenTuple
 
dimDom :: (VarLenTuple n) => PrimRec n -> Nat
dimDom = lengthT . go where
  go :: (VarLenTuple n) => PrimRec n -> Tuple n Nat
  go _ = varLenTuple
 
class VarTuple n where
  varTuple' :: Nat -> Tuple n (PrimRec m)
 
instance VarTuple Z where
  varTuple' _ = TEmpty
 
instance (VarTuple n) => VarTuple (S n) where
  varTuple' (S n) = attach (Proj n) $ varTuple' n

varTuple :: (VarLenTuple n, VarTuple n) => PrimRec n -> Tuple n (PrimRec m)
varTuple = varTuple' . dimDom 

prod :: (VarTuple n, VarLenTuple n) => PrimRec (S n) -> PrimRec (S n)
prod r = Induct (extend one) (mult %. (arg 0 %: r' %: TEmpty)) where
  r' = r %. fmap (mapProj S) (varTuple r)
 
existL :: (VarTuple n, VarLenTuple n) => PrimRec (S n) -> PrimRec (S n)
existL r = equalZ %.. prod r
 
sumpr :: (VarTuple n, VarLenTuple n) => PrimRec (S n) -> PrimRec (S n)
sumpr r = Induct (extend Zero) (add %. (arg 0 %: r' %: TEmpty)) where
  r' = r %. fmap (mapProj S) (varTuple r)
 
forallL :: PrimRec N1 -> PrimRec N1
forallL r = equalZ %.. sumpr r

$prime(x) = x > 1 \land \neg (\exists u < x: \exists v < x: u v = x)$はxが素数であるという述語を表す.
これを定義する.

primep :: PrimRec N1
primep = andp %. (lneq %. (extend one %: arg 0 %: TEmpty) %:
                  notp %.. (existL r' %. (arg 0 %: arg 0 %: TEmpty)) %: TEmpty)
  where
    r' :: PrimRec N2
    r' = existL r'' %. (arg 1 %: arg 0 %: arg 1 %: TEmpty)
 
    r'' :: PrimRec N3
    r'' = equal %. (mult %. (arg 0 %: arg 1 %: TEmpty) %: arg 2 %: TEmpty)

さて, $\mu_{z < y}p(z,\vec{x}) = \min\{z\in \mathbb{N}; z < y \land p(z,\vec{x})\} \cup \{y\}$とおく. この函数はy未満のzで$p(z,\vec{x})$を満たすものの中で最小の数(そのようなものが存在しなければy)を返す函数である.
この函数はprfとして, $\mu_{z < y}p(z,\vec{x}) = \sum_{v < y} g(v,\vec{x}) ; \;\; g(v,\vec{x}) = \exists z\leq v: p(z,\vec{x})$で定義できる.

mu :: (VarTuple n, VarLenTuple n) => PrimRec (S n) -> PrimRec (S n)
mu p = sumpr g %. (varTuple p) where
  f :: Tuple (S n) (PrimRec m) -> Tuple (S n) (PrimRec m)
  f (TSuc x xs) = TSuc (Suc %.. x) xs
 
  g :: (VarTuple n, VarLenTuple n) => PrimRec (S n)
  g = existL p %. (f $ varTuple p)

最後に, 目標だった函数prは,
$pr(0) = 2; pr(n+1) = \mu_{z < pr(x)!+2}(pr(x) < z \land prime(z))$
とすれば得られる.

factorial :: PrimRec N1
factorial = Induct (extend one) (mult %. (arg 0 %: Suc %.. arg 1 %: TEmpty))
 
pr :: PrimRec N1
pr = Induct (extend two) (mu p %. (add %. (extend two %: factorial %.. arg 0 %: TEmpty) %: arg 1 %: TEmpty)) where
  p :: PrimRec N2
  p = andp %. (lneq %. (pr %.. arg 1 %: arg 0 %: TEmpty) %: (primep %.. arg 0) %: TEmpty)

実行

上の函数が正しい挙動をするかどうかを手で検証するのは大変なのでQuickCheckを使った.
なお, 直接実行するとprやfactorialなどはすごく時間がかかる(私のPCだとpr(4)ですでに数分では終わらない)ので本当にこれで正しい出力を返すかどうかは怪しい.

コード全体


参考文献

計算論 計算可能性とラムダ計算 (コンピュータサイエンス大学講座)

計算論 計算可能性とラムダ計算 (コンピュータサイエンス大学講座)