Hatena::ブログ(Diary)

Life Goes On このページをアンテナに追加 RSSフィード

2012-02-09

[][] Haskell vs F# その後

mkotha さんに直してもらったりして、Haskellのコードはだいぶ速くなりました。どうも2重ループの内側がボトルネックのようなので、そこを展開して、データ構造も変えて、UNPACKプラグマは効くので残して、正格評価を1ヶ所だけ。性能と可読性のバランスがそこそことれたかなと思ってます。C++ や F# のコードにも同じような改修を加えたら、Haskell はまた抜かれてしまいました。まぁでも、目くじら立てるほどの差でもないので、そのままにしています。

実行環境が Windows というアドバンテージがあるとはいえ、C++ も超える F# の健闘が光ります。明示的な副作用がない関数プログラミングでこれだけ速いとうれしい。コード書いてても気持ちがいいし、Microsoft でなければもっと流行っていいはず。

最終形のコードを以下に載せておきます。

ついでに Scala でも書いてみました。C++ の2倍くらいの時間でしたが、書くだけで疲れてチューニングする元気もないので、処理時間はなしでコードだけ。

f:id:rst76:20120209232247p:image

Haskell のコード

Data.Vector を早く標準モジュールに採用してほしい。cabal install するくらいはいいけど、プロファイラが使えないのが痛い。cabal install するときは -p オプションを忘れずに。

import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U

data Node = Node {
  df :: {-# UNPACK #-} !Double,
  pu :: {-# UNPACK #-} !Double,
  pm :: {-# UNPACK #-} !Double,
  pd :: {-# UNPACK #-} !Double,
  iu :: {-# UNPACK #-} !Int,
  im :: {-# UNPACK #-} !Int,
  id :: {-# UNPACK #-} !Int
  }

induceBackward :: V.Vector Node -> U.Vector Double -> U.Vector Double
induceBackward nodes values = values `seq` V.convert (V.map value nodes)
  where
  next k = U.unsafeIndex values (U.length values `div` 2 + k)
  value (Node df pu pm pd iu im id) = (pu * next iu + pm * next im + pd * next id) * df

iteration = 10000

main :: IO()
main =  print (maximum [value i | i <- [0..iteration-1]])
  where
  value i = (foldr induceBackward (lastValues i) testTree) U.! 0
  lastValues i = U.replicate 201 (fromIntegral i)
  createNode j = Node 1.0 (1.0/6.0) (2.0/3.0) (1.0/6.0) (j+1) j (j-1)
  testTree = map (\i -> V.fromList (map createNode [-i..i])) [0..99]

F# のコード

自分の中で No.2 言語の座に躍り出ました。はてな記法も対応してほしい。

[<Struct>]
type ShortRateTreeNode = 
  val DF : double
  val PU : double
  val PM : double
  val PD : double
  val IU : int
  val IM : int
  val ID : int

  new (df, pu, pm, pd, iu, im, id) = {
    DF = df
    PU = pu
    PM = pm
    PD = pd
    IU = iu
    IM = im
    ID = id
    }

let induceBackward (nodes : ShortRateTreeNode []) (values : double []) : double [] =
  let next k = values.[values.Length / 2 + k]
  let value (node : ShortRateTreeNode) =
    (node.PU * next node.IU + node.PM * next node.IM + node.PD * next node.ID) * node.DF
  Array.map value nodes

let ITERATION = 10000

let maturityValues i = Array.create 201 (double i)
let createNode j = ShortRateTreeNode(1.0, 1.0/6.0, 2.0/3.0, 1.0/6.0, j+1, j, j-1)
let testTree = Array.init 100 (fun i -> Array.map createNode [|-i..i|])
let value i = (Array.foldBack induceBackward testTree (maturityValues i)).[0]
stdout.WriteLine (Array.max (Array.init ITERATION value))

C++ のコード

性能の話になると、けっきょく引き合いに出されるんですよね。

# include <vector>
# include <utility>
# include <memory>
# include <cstdio>

using namespace std;

struct node
{
  double df;
  double pu;
  double pm;
  double pd;
  int iu;
  int im;
  int id;
};

auto_ptr<vector<double> > induce_backward(const vector<node> &nodes, const vector<double> &values)
{
  const int n_nodes = nodes.size();
  auto_ptr<vector<double> > ret(new vector<double>());
  ret->reserve(n_nodes);
  const int n = values.size() / 2;
  for(int j = 0; j < n_nodes; j++)
  {
    const node &node = nodes[j];
    double sum = (node.pu * values[n + node.iu] + node.pm * values[n + node.im] + node.pd * values[n + node.im]) * node.df;
    ret->push_back(sum);
  }
  return ret;
}

const int iteration = 10000;

int main()
{
  vector<vector<node> > test_tree;
  for(int i = 0; i < 100; i++)
  {
    test_tree.push_back(vector<node>());
    vector<node> &nodes = test_tree.back();
    for(int j = -i; j <= i; j++)
    {
      node node = struct node();
      node.df = 1.0;
      node.pu = 1.0/6.0;
      node.pm = 2.0/3.0;
      node.pd = 1.0/6.0;
      node.iu = j+1;
      node.im = j;
      node.id = j-1;
      nodes.push_back(node);
    }
  }
  double maxval = 0;
  for(int i = 1; i <= iteration; i++)
  {
    auto_ptr<vector<double> > values(new vector<double>(201, (double)i));
    for(int j = 99; j >= 0; j--)
      values = induce_backward(test_tree[j], *values);
    maxval = max(maxval, (*values)[0]);
  }
  printf("%f\n", maxval);
  return 0;
}

Scala のコード

3年くらい前に触ったときと比べると、API が格段に充実しています。でも疲れた・・・

object ScalaBackwardInductionTest {
  class Node(j : Int) {
    val df = 1.0
    val pu = 1.0 / 6.0
    val pm = 2.0 / 3.0
    val pd = 1.0 / 6.0
    val iu = j + 1
    val im = j
    val id = j - 1
  }

  def induceBackward(nodes : Array[Node], values : Array[Double]) : Array[Double] = {
    def next(k : Int) : Double = values(values.length / 2 + k)
    def value (node : Node) : Double =
      (node.pu * next(node.iu) + node.pm * next(node.im) + node.pd * next(node.id)) * node.df
    return nodes.map(value)
  }

  val iteration = 10000

  def main(args : Array[String]) {
    def createStep(i : Int) : Array[Node] = Array.range(-i, i+1).map(j => new Node(j))
    val testTree : Array[Array[Node]] = Array.range(0, 99+1).map(createStep)
    def lastValues(i : Int) : Array[Double] = Array.fill(201)(i.toDouble)
    def value (i : Int) : Double = testTree.foldRight(lastValues(i))(induceBackward)(0)
    println((0 to iteration-1).map(value).max)
  }
}

2012-02-04

[][] Haskell vs F#

VB .NET で書かれたプログラムを速くしろと言われて、Haskell と F# で書き換えたりしています。僕の目論見ではHaskellの方が速くなるはずで、そしたら『F# もいい言語ですけどね、やはり速度を求めるならネイティブアプリケーションでないと。』とか何とか言って、Haskell 化してしまう予定でした。だって僕は F# 初心者だし、HaskellC++ にも負けないわけだし、Haskell が速いに決まってるじゃないですか。

ところが実際に試してみたら、F# の方が速くなってしまいました。それも2倍以上。というわけで、図らずも自分自身の Haskell 力の無さを露呈してしまったわけですが、それはともかく、Haskell が遅いとかいう汚名を着せられてしまうのは、避けなくてはなりません。

でも、ごめんなさい、僕の能力では限界です。

世の Haskeller のみなさま、なんとかこのコードを高速化あるいは、この類のコードが遅くなる理由を教えていただけないでしょうか。STUArray は面倒とか、unsafeなんちゃらは避けたいとか、無いこともないですが、まずは速くなることを確認したい。

それから F#er のみなさま、もっと速くなるよとか、こう書いた方がいいよとかあれば、忌憚なきご意見を承りたく。代入使って速いのはいいとしても、ちょっと書き方変えるだけでこれだけ結果が変わるのは意外です。

実行時間

Haskell、F# それぞれのコードの実行時間はおよそ以下の通りで、F# は何種類か書いてみました。グラフが切れてますが、F# の一番遅いやつは、一番速いやつの240倍くらいの時間がかかります。

あと構造体を使うと、F# のコードはさらに1割くらい速くなったりするのですが、気づかなかったことにしてます。

f:id:rst76:20120204013937j:image

Haskell のコード

import Data.Array.Unboxed

data Node = Node {
  df :: Double,
  branch :: [(Int, Double)]
  }

induceBackward :: Array Int Node -> UArray Int Double -> UArray Int Double
induceBackward nodes values = accumArray (+) 0 (bounds nodes)
  [(j, p * values ! k * df) | (j, Node df branch) <- assocs nodes, (k, p) <- branch]

iteration = 1000

main :: IO()
main = print (maximum [value i | i <- [0..iteration-1]])
  where
  value i = foldr induceBackward (lastValues i) testTree ! 0
  lastValues i = listArray (-100, 100) (repeat (fromIntegral i))
  testTree = [listArray (-i, i)
    [Node 1.0 [(j-1, 1.0/6.0), (j, 2.0/3.0), (j+1, 1.0/6.0)] | j <- [-i..i]]
    | i <- [0..99]]

F# のコード

type Node = 
  val DF : double
  val Branch : (int * double) []
  
  new (df, branch) = {
    DF = df
    Branch = branch
    }

let induceBackwardVerySlow (nodes : Node []) (values : double []) : double [] =
    let n = values.Length / 2
    [| for node in nodes ->
        Array.sum [| for (k, p) in node.Branch -> p * values.[n + k] * node.DF |] |]

let induceBackwardSlow (nodes : Node []) (values : double []) : double [] =
    let n = values.Length / 2
    [| for node in nodes ->
        Array.sum (Array.map (fun (k, p) -> p * values.[n + k] * node.DF) node.Branch) |]

let induceBackwardMedium (nodes : Node []) (values : double []) : double [] =
    let n = values.Length / 2
    let value (node : Node) =
        Array.sum (Array.map (fun (k, p) -> p * values.[n + k] * node.DF) node.Branch)
    Array.map value nodes

let induceBackwardFast (nodes : Node []) (values : double []) : double [] =
    let n = values.Length / 2
    let arr = Array.zeroCreate nodes.Length
    for j in 0 .. nodes.Length-1 do
        let node = nodes.[j]
        for k in 0 .. 2 do
            let next = node.Branch.[k]
            arr.[j] <- arr.[j] + snd next * values.[n + fst next] * node.DF
    arr

let ITERATION = 1000

let lastValues i = Array.create 201 (double i)
let testTree =
    [| for i in 0..99 ->
        [| for j in -i..i ->
            Node(1.0, [|(j-1, 1.0/6.0); (j, 2.0/3.0); (j+1, 1.0/6.0)|]) |] |]
let value i = (Array.foldBack induceBackwardFast testTree (lastValues i)).[0]
stdout.WriteLine (Array.max (Array.init ITERATION value))

2011-06-11

[] 冪集合

2011-06-08 - HaHaHa!にあったお題に挑戦。

まずはごくごく素直に。nobsunがなぜNatとか定義しているのか、よく分かっていません。

data Set = S [Set]

instance Show Set where
  show (S xs)
    | null xs = "0"
    | otherwise = "{" ++ drop 2 (concat [", " ++ show x | x <- xs]) ++ "}"

main :: IO ()
main = interact $ show . powerSet . read

powerSet :: Int -> Set
powerSet 0 = S []
powerSet n = S $ power xs where
  S xs = powerSet $ n-1

power :: [Set] -> [Set]
power [] = [S []]
power (x:xs) = [S z | S y <- power xs, z <- [y, x:y]]
ghci> powerSet 0
0
ghci> powerSet 1
{0}
ghci> powerSet 2
{0, {0}}
ghci> powerSet 3
{0, {0}, {{0}}, {0, {0}}}

いろいろな解があると面白そうだったので、あなごるにも出題しました。

anarchy golf - Power Set

プログラマのみなさん、ゴルファーも非ゴルファーもぜひぜひご参加を。

2011-05-23

[] Function call expression 参戦記

公私ともにバタバタしていて、blogtwitter も放置していますが、どうにか生きてます。

anarchy golf - Function call expression に関して、youzさんが解説を書けと言ってるので、2位という立場で僭越ながら書いてみる。

2週間前に問題を見て、なんて好みの問題だろうと思い、誰も参加してないけど submit。

こういう処理系ちっくな問題のときはだいたい、Parsec で考えを整理してから実装してる。今回もまずParsec で下書き。

import Text.ParserCombinators.Parsec

m @ main = getLine >>= putStrLn . either show id . parse expr "" >> m

expr = do
  f <- char 'f'
  xs <- many $ do
    char '('
    e <- expr
    char ')'
    return e
  return $ foldl apply [f] xs

apply f x = '(' : f ++ ' ' : x ++ ")"

これを Parsec 使わずに書きなおす。Programming in Haskell にもある通り、パーサーは文字列を受け取って、パースした結果と残りの文字列を返せばいい。いろいろ小細工して 98B。

m@main=getLine>>=putStrLn.(!1)>>m
(_:'(':a)?b=a!0?('(':b++' ':a!1++")")
(_:a)?b=[a,b]
a!i=a?"f"!!i

もう少し縮みそうな気がするも、競争相手がいないので放置していたところ、1週間前に notogawaさんが以下のツィート。

@Lost_dog_ さんが嘆いていたのでFunction call expressionを処理

そうですか。こちらが身を削って書いたコードを易々と“処理”されたのですね。しかも 91B。7B 差。

悔しさに歯噛みしつつ、アイデアも時間もなくそのまま deadline 当日を迎える。

前夜に GCJ Round1B があり、寝不足の頭と残念な結果をお供に、子供と鉄道博物館へ。

半分諦めていたけど、運転シミュレータに並びながらあれこれ考える。

今の方針で縮めるとしたら getLine と putStrLn で1行ずつ処理してるところくらいだけど、局所最適化してるので、このまま縮めるのは無理だろう。スタックマシンの処理系実装みたいなことをやればもう少しシンプルに書けるのでは?

そんなアイデアをベースに、プラレールに興じる子供の横で、鉄道手帳にメモしながら検討する。文字 'f' が来たらスタックに積む。'(' はたぶん無視しても大丈夫。')' が来たらスタックのトップの2要素を関数適用して "(f x)" みたいな文字列にすればよさそう。さらに '\n' も '(' と同じように無視して、入力をまとめて処理すれば、入力1行に対して解析結果がスタックに1行積み上がるから、反転して出力すれば求める結果が得られるんじゃないか?!

"f(f)(f)" とか "f(f(f))" みたいな入力を当てはめて上手くいきそうだったので、コーディングにとりかかる。といってもPCなど持ってないので、携帯からあなごるのサイトにアクセスして、フォームで入力。記号の入力が激しく面倒でした。でもやればできるもんです。コンパイルエラーもなく一発で success !92B!惜しい。

main=interact$unlines.reverse.foldl(&)[]
s&'f'="f":s                       -- 'f'が来たらスタックに積む
(b:a:s)&')'=('(':a++' ':b++")"):s -- ')'が来たらスタックの先頭2要素を関数適用
s&_=s                             -- 残り('(', '\n')は無視

あと 1B をどうすればいいか、残り30分でつらつら考えましたが、時間切れ。

notogawaさんの解を見て納得。'\n' もスタックに積んじゃえば、そのままスタックの文字を結合して答えが出るんだ。

あと、改めて考えると、')' は後置(逆ポーランド記法関数適用演算子と捉えられるんですよね。面白いなぁ。

負けたことは悔しいですが、とても楽しい問題でした。infix to postfix も同じように縮まないだろうかと考えたり。

2010-08-17

[] ポケットキューブ(2×2×2キューブ)の最短手数

ルービックキューブを揃えるための最短手数は最大でも20手だそうです。3×3×3を解くには相当なリソースが必要ですが、2×2×2なら何とかなるんじゃないかと思い、プログラムを組んでみました。

結果は↓の通りで最大で11手でした。これくらい誰か既に計算してそうですが、見つかりません。我こそはという方は追試をお願いします。

距離
01
19
254
3321
41847
59992
650136
7227536
8870072
91887748
10623800
112644

プログラムはこんな感じです。

キューブの状態を表すのに、最初は展開図っぽいデータ構造にしていたのですが、けっこうメモリを食うので、8個のコーナーキューブ(a〜h)がどの場所にあってどっちを向いているかを保持するようにしました。

キューブの向きは厳密には各々24通りありますが、キューブの場所で向きが3通りに限定されるので、0〜2の整数で表現しています。

また左下奥のキューブ(h)を位置、向きとも固定することで、重複を排除しています。鏡像とかは特に考慮していません。

import Data.List (foldl')
import Data.Set (Set, singleton, member, insert)

data Axis = F | R | U
type Move = (Axis, Int)
type Cube = [Int]

rotate :: Move -> Cube -> Cube
rotate move [a,b,c,d,e,f,g,h] = case move of
  (F,1) -> [c%1,b,g%2,d,a%2,f,e%1,h]
  (F,2) -> [g,b,e,d,c,f,a,h]
  (F,3) -> [e%1,b,a%2,d,g%2,f,c%1,h]
  (R,1) -> [e%2,a%1,c,d,f%1,b%2,g,h]
  (R,2) -> [f,e,c,d,b,a,g,h]
  (R,3) -> [b%2,f%1,c,d,a%1,e%2,g,h]
  (U,1) -> [b,d,a,c,e,f,g,h]
  (U,2) -> [d,c,b,a,e,f,g,h]
  (U,3) -> [c,a,d,b,e,f,g,h]
  where
  x % i = let d=rem x 3 in x - d + rem(d+i)3

next :: (Set Cube, [Cube]) -> (Set Cube, [Cube])
next (vs0,xs) = foldl' add (vs0,[]) [rotate (a,r) x | x<-xs, a<-[F,R,U], r<-[1..3]]
  where
  add (vs,ys) y
    | member y vs = (vs, ys)
    | otherwise   = (insert y vs, y:ys)

ini :: Cube
ini = [0,3..21]

countAll :: [Int]
countAll = takeWhile (/=0) $ map (length.snd) $ iterate next (singleton ini, [ini])

main :: IO ()
main = print $ countAll

2010-07-06

[] Haskellと確率と限定継続

会社の勉強会で、Haskell による確率プログラミングについて話しました。

元ネタは Oleg さんの論文HANSEI: Embedded Probabilistic Programming)です。OCaml で書かれていたプログラムHaskell に翻訳しながら解説しました。

もともとの論文の主張は、木探索として定義したモデルを継続渡し方式(Continuation-Passing Style / CPS)にすることで効率化できる、さらに限定継続を利用することで、継続を意識せず直接方式(Direct Style)で書けるようになるというものでした。

木探索から CPS にすることで効率化という流れは Haskell でも同じなのですが、限定継続を利用しても直接方式では書けないため、あまり嬉しくありません。でも Haskell には do 記法があるので、モナドとして定義すれば擬似的に直接方式で書けてハッピーというのが、今回の発表の主旨です。

なにぶん継続初心者が書いていますので、誤りなどあれば指摘いただけると助かります。

ソースコード

ソースコードgithub に上げました。まだ使い方がよく分かっていません。

rst76/probability - GitHub

限定継続について

限定継続を学ぶにあたっては、以下の資料がとても分かりやすく参考になりました。また私の発表資料中で、面識もない浅井先生のお名前を勝手に拝借しています。すみません。

継続を使った Printf の型付け

Haskell での限定継続の実装

Haskell での限定継続の実装は以下のようにいくつかありました。今回は 2 つ目の実装(継続モナドを拡張したもの)を簡略化して、利用しています。

2010-05-09

[] Google Code Jam 2010 - Qualification Round

今年も参加してます。予選だしあまり言うこともないのですが、せっかくなのでコードを晒しときます。

去年のコードと比べると、ゴルフの影響を受けすぎ(悪い意味で)。もうちょっと可読性を大事にする人間だったはずなのですが・・・。

ともあれ Haskell で解けるというのはウレシイ限りです。

A

最初、問題の意味がぜんぜん分かりませんでした(snapper って何だよ!)。分かれば簡単。2 進表記で指定された桁以下が全部 1 かどうか答える。

main=interact$unlines.zipWith(++)["Case #"++shows i": "|i<-[1..]].map(f.map read.words).tail.lines
f[n,k]=["OFF","ON"]!!(0^mod(k+1)(2^n))

B

与えられた数列の周期を求めて、0 以上の最小値を答える。これもまぁ簡単。

main=interact$unlines.zipWith(++)["Case #"++shows i": "|i<-[1..]].map(f.map read.words).tail.lines
f(_:x:ys)=show$mod(-x)$foldl1 gcd[abs$x-y|y<-ys,y/=x]

C

グループ単位でジェットコースターに乗るのを繰り返すとき、何人乗れるか答える。Small はともかく Large で 10^8*1000 ってフツウにやったら終わらないよなぁ、面倒だなぁとか思いながら書いてみたら、やっぱりものすごく面倒でした。なぜこの問題だけ?

人数のリストが途中から循環するので、循環が始まる点を適当に求めて、それ以前/以後で場合分けして求めてます(6 〜 7 行目の a というやつがそれです)。さすがにこれはもっと簡単に書けるのでは。他の人のコードを読みたい。

import Data.List
main = interact $ unlines.zipWith(++)["Case #"++shows i": "|i<-[1..]].f.map(map read.words).tail.lines
f([r,k,n]:g:x) = show a : f x
  where
  -- 求める人数は循環前の部分の人数+循環部の人数(完全に一周する分)+循環部の人数(途中まで辿る分)
  a | r>d0 = sum g0 + sum g1*((r-d0)`div`d1) + sum(genericTake((r-d0)`mod`d1)g1)
    | otherwise = sum $ genericTake r g0
  -- d0, d1 は g0, g1 の長さ
  [d0,d1] = map genericLength[g0,g1]
  -- g0, g1 は人数のリストの、循環前の部分と循環部
  [g0,g1] = map(snd.unzip)[g0'',g1'']
  g1'' = b : takeWhile(/=b)g2''
  (g0'',_:g2'') = span(/=b)g''
  -- b は循環部に含まれる点。ここから循環部が始まるとする。
  b = g''?tail g''
  -- リストの循環部に含まれる点を見つける。ウサギとカメのロジック
  (x:xs)?(y:_:ys)
    | x==y = x
    | otherwise = xs?ys
  -- g'' の要素は一度にコースターに乗れる人数(ID つき)
  g'' = (0!0)g'
  -- 一度にコースターに乗れる人数を求めてリストにする。
  (s!m)((i,x):y)
    | s+x>k||m+1>n = (i,s) : (0!0)((i,x):y)
    | otherwise = (s+x)!(m+1) $ y
  -- グループのリスト g に ID を振る
  g' = cycle $ zip[0..]g
f _ = []

2010-04-16

[] Haskellers Meeting 2010 Spring

行ってきました。

Haskell の神さまにたくさん会いました。

サインももらいました。オレオレ関数型言語を実装したいって言ったら、『今はこういう実装じゃないけど最初に読むにはいい本だよ』みたいなことを言って、↓のメッセージを書いてくれました。

家宝にしようかと思ったけど、それじゃ意味がないので、ちゃんと読もうと思います。

f:id:rst76:20100416225221j:image:w320

2010-03-02

[] SPOJ の時間制限が厳しい件について

Lost Dog さんの記事が面白そうだったのと、あなごるサーバが落ちてて何もできなかったのとで、SPOJ に初挑戦。

Sphere Online Judge (SPOJ) - Problem TRT

問題は Problem 67 - Project Euler に近いと思います。O(n2) の DP(?)で解きました。

ただ、何も工夫しないと 1s で終わらず TLE になってしまうので、

  1. インデックスアクセスしてるところはリストから配列
  2. 正格評価を使う
  3. foldl' とか zipWith とか使ってたのをやめてベタに再帰

などで高速化を図りました。

なんとか accept されましたが、もうちょいスマートな解き方がないもんでしょうか?

以下にソースコードを載せておきます。可読性は著しく低いので、ネタバレにはならないかと。

続きを読む

[] Haskell(GHC)でのプロファイルのとりかた

いつも忘れて、その度にぐぐってるので、メモ。

コンパイル時に -prof -auto-all オプションをつける

ghc --make -prof -auto-all hoge.hs

実行時に、+RTS -p オプションをつけて、出力される 〜.prof ファイルを読む

./hoge.exe +RTS -p -RTS

特定の式のコストを見たいときは、ソースコード中でコスト集約点を指定(Set Cost Center)する

hoge = {-# SCC "hoge" #-} <expression>

詳しくは、第5章プロファイルを取る を参照

2010-02-03

[] 鳥たちと戯れる

To Mock a Mockingbird を読んでいます。言わずと知れたコンビネータ論理の教科書(?)です。

そこで算術や論理を関数で表すという、ラムダ計算でよくある話が登場するのですが、自然数の定義が一般的なチャーチ数と違っていて分かりにくかったので、検算用のコードを書きました。

ふつうのチャーチ数では、zero f x = x = KIfx、one f x = f x = Ifx、two f x = f (f x) = SBIfx と定義しますが、zero = I、one = V(KI)I、two = V(KI)(V(KI)I) という定義になってます。

こんな感じ(↓)で遊べます。

Main> toInt $ exp five (multi three two)     
15625
Main> toBool $ greater (multi two five) (multi three three)
True
Main> toInt $ concat (multi three four) three              
123

コンビネータの本なので、下にあげた plus とか greater みたいな再帰関数コンビネータで(ポイントフリーで)書く方法なんかも載っていて、勉強するにはとてもいい本だと思います。お薦め。

ちょっと確認する程度ならこんなコードで十分だけど、unsafeCoerceも鬱陶しいし、やはりここはちゃんとした処理系の実装に挑戦しようかな。

import Prelude hiding (succ,pred,exp,even,not,and,or,length,concat)
import Unsafe.Coerce

u = unsafeCoerce

-- 事前準備 各種コンビネータ
i x = x
k x y = x
t x y = y x
l x y = x (y (u y))
b x y z = x (y z)
c x y z = x z y
v x y z = z x y
r x y z = y z x
s x y z = x z (y (u z))

-- 論理
true  = k
false = k i
not   = v false true
and   = r false
or    = t true
imply = r true
equiv = c s not

-- 算術
zero   = i
succ   = v false
one    = succ zero
two    = succ one
three  = succ two
four   = succ three
five   = succ four
ten    = multi two five
pred   = t false
isZero = t true

-- 各種関数
plus    m n = isZero m n     $ plus (pred $ u m) (succ $ u n)
multi   m n = isZero m zero  $ plus  n (multi (pred $ u m) n)
exp     m n = isZero n one   $ multi m (u exp m (pred $ u n))
even      n = isZero n true  $ not $ u even $ pred $ u n 
greater m n = isZero m false $ isZero n true $ greater (pred $ u m) (pred $ u n)
length    n = a1 zero n
  where a1 l n = a l n (u l) $ u a1 (succ l) n
        a  l n = greater (exp ten l) n
concat  m n = plus (multi m (exp ten (length $ u n))) n

-- 表示・デバッグ用関数
toBool b = b True False :: Bool
toInt  n = isZero n 0 $ (+) 1 $ toInt $ pred $ u n :: Int
To Mock a Mockingbird

To Mock a Mockingbird