関数プログラミングのアプローチ (7)

モンテカルロ法の関数的プログラム (2)

前回、遅延リストを作成しました。これを用いて前に作成したモンテカルロ法のプログラムを書き換えたいと思います。しかし、プログラムを書き換える前に、次の２つのコードを見てください。

x, ls = next(ls)

x, rs = rand(rs)

上のコードは遅延リストの要素にアクセスするコードで、下は乱数状態オブジェクトから次の乱数を生成するコードです。並べて比較すると、この２つは全く同じであることが分かります。つまり、乱数は遅延リストとして表現することができるのです。下のコードは、乱数を遅延リストとして実装したものです。

def rand ():
    def f (x,y,z,w,t):
        t=(x^(x<<11)) & 0xffffffff
        x=y & 0xffffffff
        y=z & 0xffffffff
        z=w & 0xffffffff
        w=(w^(w>>19))^(t^(t>>8)) & 0xffffffff
        return w, LazyList(lambda: f(x,y,z,w,t))
    return LazyList(lambda: f(123456789, 362436069, 521288629, 88675123, 1))

これを用いてモンテカルロ法のプログラムを書くと次のようになります。

def main ():
    print math.sqrt(6.0 / estimate(cesaro_test, 10000, groupL(2, rand())))

def estimate (test, times, rands):
    x, y = 0.0, 0.0
    for b in forceList(takeL(times, mapL(test, rands))):
        if b: x = x + 1
        y = y + 1
    return x / y

def cesaro_test ((x,y)):
    return 1 == gcd(x, y)

def gcd (a, b):
    while b > 0: a, b = b, a % b
    return a

最後に、プログラムの全文をまとめて掲載しておきます。また、比較のために同じものを関数型言語 Clean で書いたものを併記します。専用の言語を用いることで、同じプログラムがより簡潔に書けることを見てください。

Cleanによるモンテカルロ法の関数的プログラム

module MonteCarlo

import Int, Real, Misc

Start = sqrt (6.0 / estimate cesaro_test 10000 (group 2 rand))

estimate test times rands
    # (x,y) = foldl accm (0.0,0.0) $ take times $ map test rands
    = x / y
  where
    accm (x,y) b = if b (x+1.0,y+1.0) (x,y+1.0)

cesaro_test [x,y] = 1 == gcd x y

//gcd a 0 = a
//gcd a b = gcd b (a rem b)

rand = f 123456789 362436069 521288629 88675123 1
  where
    f x y z w t
        # t=(x bitxor (x << 11))
          x=y
          y=z
          z=w
          w=(w bitxor (w >>> 19)) bitxor (t bitxor (t >>> 8))
        = [w: f x y z w t]

map f [] = []
map f [a:ls] = [f a: map f ls]

take 0 _  = []
take _ [] = []
take n [a:ls] = [a: take (n - 1) ls]

group n ls
    # (x,ls) = split n ls
    = [x: group n ls]
  where
    split 0 ls = ([],ls)
    split _ [] = ([],[])
    split n [a:ls] # (x,ls) = split (n - 1) ls
                   = ([a:x],ls)

foldl f x [] = x
foldl f x [a:ls] = foldl f (f x a) ls

Pythonによるモンテカルロ法の関数的プログラム

# -*- coding: utf-8 -*-
import math

def main ():
    print math.sqrt(6.0 / estimate(cesaro_test, 10000, groupL(2, rand())))

def estimate (test, times, rands):
    x, y = 0.0, 0.0
    for b in forceList(takeL(times, mapL(test, rands))):
        if b: x = x + 1
        y = y + 1
    return x / y

def cesaro_test ((x,y)):
    return 1 == gcd(x, y)

def gcd (a, b):
    while b > 0: a, b = b, a % b
    return a

def rand ():
    def f (x,y,z,w,t):
        t=(x^(x<<11)) & 0xffffffff
        x=y & 0xffffffff
        y=z & 0xffffffff
        z=w & 0xffffffff
        w=(w^(w>>19))^(t^(t>>8)) & 0xffffffff
        return w, LazyList(lambda: f(x,y,z,w,t))
    return LazyList(lambda: f(123456789, 362436069, 521288629, 88675123, 1))

class LazyList:
    def __init__ (self, fun):
        self.undef = True
        self.fun = fun
    def __next__ (self):
        if self.undef:
            self.value = self.fun()
            self.undef = False
        return self.value

def next (x):
    return x.__next__()

def mapL (fun, ls):
    def f (ls):
        v = next(ls)
        if v:
            x,ls = v
            return fun(x), LazyList(lambda: f(ls))
        else:
            return v
    return LazyList(lambda: f(ls))

def takeL (n, ls):
    def f (n,ls):
        if n <= 0:
            return None
        v = next(ls)
        if v:
            x,ls = v
            return x, LazyList(lambda: f(n-1,ls))
        else:
            return v
    return LazyList(lambda: f(n,ls))

def groupL (n, ls):
    def f (ls):
        ret = []
        for i in xrange(0,n):
            v = next(ls)
            if v:
                x, ls = v
                ret.append(x)
            elif ret:
                return ret, LazyList(lambda: v)
            else:
                return v
        return ret, LazyList(lambda: f(ls))
    return LazyList(lambda: f(ls))

def forceList (ls):
    ret = []
    v = next(ls)
    while v:
        x,ls = v
        ret.append(x)
        v = next(ls)
    return ret