原始帰納的函数で素数を計算
原始帰納的函数; 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の構成
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)ですでに数分では終わらない)ので本当にこれで正しい出力を返すかどうかは怪しい.
コード全体
参考文献
計算論 計算可能性とラムダ計算 (コンピュータサイエンス大学講座)
- 作者: 高橋正子
- 出版社/メーカー: 近代科学社
- 発売日: 1991/08
- メディア: 単行本
- 購入: 8人 クリック: 160回
- この商品を含むブログ (36件) を見る