例のテキスト

  「自己言及の論理と計算」
  http://www.kurims.kyoto-u.ac.jp/~hassei/selfref.pdf

の話の続きです。前回の復習をすると、

-------------------------------------------------------^
【Xn 問題】
不等号 <= が定義されている集合をPos と呼ぶ。 a, b が 
Pos のとき、power(a, b) で b から a への単調増加関数〔つ
まり x <= y なら f(x) <= f(y) を満たす〕f 全体を表す。更
に、power(a, b) の要素 f, g に対して、「全ての x に対して 
f(x) <= g(x)」となるとき「f <= g」と不等号を定義する。こ
のとき power(a, b) も Pos になる。さて、Pos である a に
対して Pos である X(a)(n), n = 0, 1, 2,.. を

  X(a)(0) = 1点 I からなる自明な Pos
  X(a)(n+1) = power(a, X(a)(n)), n = 0, 1, 2, ...

と帰納的に定義する。

例えば、X(a)(1) は a と同一視でき、X(a)(2) は単調増加写像 
a -> a 全体と同一視できる。

(問題1)X([0, 1, 2])(3) を図示せよ。

-- 出力は graphviz で画像にすることにする。
-- http://www.research.att.com/sw/tools/graphviz/
-------------------------------------------------------$

という問題でした。これは、単調増加写像 f : a -> a の不動
点を求める問題と関連しています。


そして、今回はテキストの  p.17 にあるような、各 Xn 間の写
像を表してみようという問題を考えてみます。

-------------------------------------------------------^
【X(n) <-> X(n+1) 問題】
a を Pos とするとき、n = 0, 1, 2,.. に対して 2 つの写像:

  e(a)(n) : X(a)(n) -> X(a)(n+1)
  p(a)(n) : X(a)(n+1) -> X(a)(n)

を次のように帰納的に定義する。

  e(a)(0) = (a の最小値への定数関数)
  p(a)(0) = (I への定数関数)
  e(a)(n+1) = (f に合成関数 f . p(a)(n) を対応させる写像)
  p(a)(n+1) = (f に合成関数 f . e(a)(n) を対応させる写像)

(問題2) e([0, 1, 2])(2), p([0, 1, 2])(2) を図示せよ。
-------------------------------------------------------$

答えは、

  http://blade.nagaokaut.ac.jp/~sinara/haskell/xn/2-2-1.jpg

となりました。更に n を一つ増やして e([0, 1, 2)(3) と 
p([0, 1, 2])(3) でやってみると、汚いですが、

  http://blade.nagaokaut.ac.jp/~sinara/haskell/xn/2-3-2.jpg

となります。

ちなみに、テキストの p.17 の図をちょっと広げて X([0, 1])(n) 
を n = 0 から n = 7 まで書かせたのは

  http://blade.nagaokaut.ac.jp/~sinara/haskell/xn/1-7.jpg

で、なかなか graphviz の威力には感心します。

これを眺めて、前に colim と lim が一致するような事を言ったの
が間違いだった事に気づきました。「colim の順序による完備化が 
lim になる」というのが正しいです。

ついでに a を

        D
      / \
     B     C
      \ /
        A

というダイアモンド型にしてみたら、X(a)(2) は

  http://blade.nagaokaut.ac.jp/~sinara/haskell/xn/diamond-2-2.jpg

で、X(a)(1) <-> X(a)(2) は

  http://blade.nagaokaut.ac.jp/~sinara/haskell/xn/diamond-2-1.jpg

でした。やはり 3 次元表示が欲しいところです。

プログラムは以下の通り。酒井さんのものを元にしています。

-----------------------------------^ xn.hs
import System

data Funct x y = F [(x, y)] deriving Show
instance Eq y => Eq (Funct x y) where
  F f1 == F f2 = and [b == d | ((_, b), (_, d)) <- zip f1 f2]

type Map x y = x -> y
type Map0 x y = (x -> y, [x])
type Map00 x y = (x -> y, [x], [y])

map2funct :: Map0 x y -> Funct x y
map2funct (f, u) = F [(x, f x) | x <- u]

apply :: Eq x => Funct x y -> Map x y
apply (F f) a = head [b | (a', b) <- f, a == a']

data Pos a = P [a] (a -> a -> Bool)
frgtP :: Pos a -> [a]
frgtP (P ax _) = ax

power :: Pos b -> Pos a -> Pos (Funct a b)
power (P bx le_b) (P ax le_a) = P set le where
  set = map F (foldr extend [[]] ax)
  extend a funcs = [(a,b):func | func <- funcs, b <- bx, check a b func]
  check a b func = and [(le_b b d) || not(le_a a c) | (c, d) <- func ]
  le (F f1) (F f2) = and [le_b b d | ((_, b), (_, d)) <- zip f1 f2]

data Univ a = I | U (Funct (Univ a) a)
instance Eq a => Eq (Univ a) where
  I == I     = True
  U f == U g = f == g
  _ == _     = False

posU :: Pos (Funct (Univ a) a) -> Pos (Univ a)
posU (P xs le) = P (map U xs) le0 where le0 (U x1) (U x2) = le x1 x2

powerU :: Pos a -> Map (Pos (Univ a)) (Pos (Univ a))
powerU a = posU . power a

-- Xn over a
-- xU :: Eq a => Pos a -> [Pos (Univ a)]
-- xU a = iterate (powerU a) initialPos

initialPos :: Pos (Univ a)
initialPos = P [I] le where le I I = True

-- pairs of embeddings and projections
type EPpair x = (Map x x, Map x x, Pos x, Pos x)

powerPE :: Eq a => Pos a -> Map (EPpair (Univ a)) (EPpair (Univ a))
powerPE a (e, p, x, y) = (star (p, frgtP y), star (e, frgtP x),
                           powerU a x, powerU a y)
  where
  star :: Eq a => Map0 (Univ a) (Univ a) -> Map (Univ a) (Univ a)
  star = liftU . m2f . starOperation
    where
    starOperation :: Map0 x y -> Map (Map y a) (Map0 x a)
    starOperation (f, ux) g  = (g.f, ux)

    m2f :: Eq y => Map (Map y a) (Map0 x a) -> Map (Funct y a) (Funct x a)
    m2f f = map2funct . f . apply

    liftU :: Map (Funct (Univ a) a) (Funct (Univ b) b) -> Map (Univ a) (Univ b)
    liftU f (U g) = U (f g)

epU :: Eq a => Pos a -> [EPpair (Univ a)]
epU a = iterate (powerPE a) (initialPE a) where
  initialPE :: Eq a => Pos a -> EPpair (Univ a)
  initialPE a = (const bottom1, const bottom0, x0, x1) where
    x0 = initialPos
    x1 = powerU a x0
    bottom0 = head $ frgtP $ x0
    bottom1 = head $ frgtP $ x1 -- provisory bottom

-- graphs of indices
withIdx xs = zip xs [0..]

pointshow xs no = concat [no i ++ ";\n" | (_, i) <- withIdx xs]

graphshow (P u le) no = concat [no i ++ " -> " ++ no j
  ++ " [style = dotted];\n" | ((_, i), (_, j)) <- graphij le (withIdx u)]
  where
  graphij le ui = [(fi, gj)| fi <- ui, gj <- foldr (extend_graph le fi) [] ui]
  extend_graph le (f, i) (g, j) fs
    | i == j = fs
    | (or [le g h | (h, _) <- fs]) = fs
    | (le g f) = (filter (\(h, _) -> not (le h g)) fs) ++ [(g, j)]
    | True = fs

mapshow ::(Eq a,Eq b)=> Map00 a b -> (Int->String) -> (Int->String) -> String
mapshow (f, xs, ys) xno yno = 
  concat[xno i ++" -> "++ yno (jf x) ++ ";\n" |(x, i) <- withIdx xs]
  where jf x = apply (F (withIdx ys)) $ f x

main :: IO ()
main =
  do o:n:n0:_ <- fmap (map read . update ["1", "2", "0"]) getArgs
     printGraph (epU (P [0..o] (<=))) [n0..n] where
  update (d:ds) (a:as) = a : update ds as
  update [] as = as
  update ds [] = ds

  printGraph epxy (k0:range) = do
    putStr "digraph G {\nrankdir=TB;\n"
    putStr (showSubGraph (x k0) k0)
    mapM_ (\k -> do
      putStr (mapshow (e (k-1)) (node (k-1)) (node k))
      putStr (mapshow (p (k-1)) (node k) (node (k-1)))
      putStr (showSubGraph (x k) k)) range
    putStr "}\n"
    where
      e k = (e, frgtP x, frgtP y) where (e, p, x, y) = epxy !! k
      p k = (p, frgtP y, frgtP x) where (e, p, x, y) = epxy !! k
      x k = x                     where (e, p, x, y) = epxy !! k

  showSubGraph x n = "subgraph cluster" ++ (label n) ++ "{\n"
    ++ "label = " ++ (label n) ++ "\n"
    ++ pointshow (frgtP x) (node n)
    ++ graphshow x (node n) ++ "}\n"

  node n i = "\"x" ++ show n ++ ":" ++ show i ++ "\""
  label n = "X" ++ show n
-----------------------------------$ xn.hs

このプログラムは runhugs xn.hs m n n0 で X([0..m])(n0) から 
X([0..m])(n) までのグラフを graphviz の dot ファイル形式で
出力します。内部的には任意の有限 Pos に対応してます。(外か
ら Pos を指定する方法がまだわからない。)

敢えて forall を排除した、再帰的な型を用いた方法によりまし
た。なんとなくなんですが。


ここまでやって、この手の問題が Haskell は得意かどうかという
事について、つたない感想を言うと:

・多分 Ruby で書いたほうが 5 倍速く書けた。(^^;

・遅延評価により最終的に必要のない計算をしないおかげでアル
  ゴリズムがシンプルになったところがある。

・関数は同じ計算を 2 度しないのかと思ったけど、そうでもない
  なあ。引数が自然数なら配列にすればキャッシュ代わりになる。

・関数は書けるのに関数の型がなかなか書けないことがある。単に
  気力が足りないだけだけど。特に今回のテーマでは「型合わせ」
  がかなり難しかった。

・逆に関数の型が書けるとそれに合う実装ができただけでアルゴ
  リズムも自動的についてきた様な事もあった。ちょっと物理の
  「次元解析」を思い出した。

・やはり中間結果の出力がある意味できないのがデバッグを困難
  にしている。

・写像を操作する方法が豊富なので気持ちいい。エディタ上でプ
  ログラムを推敲するより、ノートに数式を書いてそれをそのま
  まコードに引き写した方が早かったりする。と、いうことは
  やっぱり Haskell って数学向きかも。

というところです。テーマが偏っていたという説もありますが、
いい勉強になりました。でも Haskell はまだまだ謎です。


--
ML: haskell-jp / quickml.com
使い方: http://QuickML.com/