読者です 読者をやめる 読者になる 読者になる

Just $ A sandbox

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

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

原始帰納函数; 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)ですでに数分では終わらない)ので本当にこれで正しい出力を返すかどうかは怪しい.

コード全体


参考文献

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

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