M.Hiroi's Home Page
http://www.geocities.jp/m_hiroi/

Julia Language Programming

簡単なプログラム (2)

[ Home | Light | Julia ]

●ファイルの連結

#
# cat.jl : ファイルの連結
#
#          Copyright (C) 2016 Makoto Hiroi
#

function cat(fin)
    for line in eachline(fin)
        print(line)
    end
end

#
if length(ARGS) == 0
    cat(STDIN)
else
    for name in ARGS
        open(name, "r") do fin
            cat(fin)
        end
    end
end

●単語のカウント

#
# wc.jl : 単語のカウント
#
#         Copyright (C) 2016 Makoto Hiroi
#
function wc(fin)
    word = 0
    inword = false
    while (c = read(fin, Char)) != '\0'
        if (isspace(c))
            inword = false
        elseif !inword
            inword = true
            word += 1
        end
    end
    word
end

if length(ARGS) == 0
    wc(STDIN)
else
    open(ARGS[1], "r") do fin
        println(wc(fin))
    end
end

●文字列の検索

#
# grep.jl : 文字列の検索
#
#           Copyright (C) Makoto Hiroi
#
function grep(fin, key)
    for line in eachline(fin)
        if search(line, key) != 0:-1
            print(line)
        end
    end
end

if length(ARGS) == 1
    # 検索文字列のみ
    grep(STDIN, ARGS[1])
elseif length(ARGS) >= 2
    open(ARGS[2], "r") do fin
        grep(fin, ARGS[1])
    end
else
    println("usage: julia grep.jl key [filename]")
end

●文字の置換

#
# tr.jl : 文字の置換
#
#         Copyright (C) 2016 Makoto Hiroi
#
function tr(fin, s1, s2)
    while (c = read(fin, Char)) != '\0'
        n = search(s1, c)
        if n > 0; c = s2[n]; end
        print(c)
    end
end

if length(ARGS) >= 2
    if length(ARGS[1]) == length(ARGS[2])
        if length(ARGS) == 2
            tr(STDIN, ARGS[1], ARGS[2])
        else
            open(ARGS[3], "r") do fin
                tr(fin, ARGS[1], ARGS[2])
            end
        end
    else
        println("第 1 引数と第 2 引数の長さが異なる")
    end
else
    println("usage: julia tr.jl str1 str2 [filename]")
end

●文字列の置換

#
# gres.jl : 文字列の置換
#
#           Copyright (C) 2015 Makoto Hiroi
#
function gres(fin, pat, r)
    for line in eachline(fin)
        print(replace(line, pat, r))
    end
end

if length(ARGS) == 2
    gres(STDIN, ARGS[1], ARGS[2])
elseif length(ARGS) == 3
    open(ARGS[3], "r") do fin
        gres(fin, ARGS[1], ARGS[2])
    end
else
    println("usage: julia gres.jl pattern replace [filename]")
end

●ソート

#
# sort.jl : ソート
#
#           Copyright 2016 Makoto Hiroi
#

# 単純挿入ソート
function insert_sort{T}(buff::Vector{T})
    for i in 2 : length(buff)
        temp = buff[i]
        j = i - 1
        while j > 0 && temp < buff[j]
            buff[j + 1] = buff[j]
            j -= 1
        end
        buff[j + 1] = temp
    end
    buff
end

# 中央値を返す
function median3(a, b, c)
    if a > b
        if b > c
            b
        elseif a < c
            a
        else
            c
        end
    else
        if b < c
            b
        elseif a < c
            c
        else
            a
        end
    end
end

# 9 つの中から中央値を選ぶ
function median9(buff, low, high)
    m2 = div(high - low, 2)
    m4 = div(m2, 2)
    m8 = div(m4, 2)
    a = buff[low]
    b = buff[low + m8]
    c = buff[low + m4]
    d = buff[low + m2 - m8]
    e = buff[low + m2]
    f = buff[low + m2 + m8]
    g = buff[high - m4]
    h = buff[high - m8]
    i = buff[high]
    median3(median3(a,b,c), median3(d,e,f), median3(g,h,i))
end

# クイックソート
function quick_sort{T}(buff::Vector{T})
    function qsort(low::Int, high::Int)
        if high - low < 16; return; end
        pivot = median9(buff, low, high)
        i = low
        j = high
        while true
            while pivot > buff[i]; i += 1; end
            while pivot < buff[j]; j -= 1; end
            if i >= j; break; end
            buff[i], buff[j] = buff[j], buff[i]
            i += 1
            j -= 1
        end
        qsort(low, i - 1)
        qsort(j + 1, high)
    end
    #
    qsort(1, length(buff))
    insert_sort(buff)
end

if length(ARGS) == 0
    write(STDOUT, quick_sort(readlines(STDIN)))
elseif length(ARGS) == 1
    open(ARGS[1], "r") do fin
        write(STDOUT, quick_sort(readlines(fin)))
    end
else
    println("usage: julia sort.jl [filename]")
end

●ファイルのコピー

#
# cp.jl : ファイルのコピー
#
#         Copyright (C) Makoto Hiroi
#

# バイト単位のコピー
function cp(fin, fout)
    while !eof(fin)
        write(fout, read(fin, UInt8))
    end
end

# バッファに読み込む
function cp1(fin, fout)
    local n
    buff = Array(UInt8, 16)
    while (n = readbytes!(fin, buff)) == 16
        write(fout, buff)
    end
    write(fout, buff[1:n])
end

# ファイルを全部読み込む
function cp2(fin, fout)
    write(fout, readbytes(fin))
end

if length(ARGS) == 2
    open(ARGS[1], "r") do fin
        open(ARGS[2], "w") do fout
            cp1(fin, fout)
        end
    end
else
    println("usage: julia cp.jl src dst")
end

●ファイルのエントロピー

#
# entoropy.jl : ファイルのエントロピーを求める
#
#               Copyright (C) 2016 Makoto Hiroi
#
function entoropy(fin)
    count = fill(0, 256)
    while !eof(fin)
        c = read(fin, UInt8)
        count[c + 1] += 1
    end
    total = sum(count)
    e = 0.0
    for x in count
        if x == 0; continue; end
        p = x / total
        e += - p * log(2, p)
    end
    e, e * total / 8
end

if length(ARGS) == 1
    open(ARGS[1], "r") do fin
        println(entoropy(fin))
    end
else
    println("usage: julia entoropy.jl filename")
end

●鶴亀算

  1. 鶴と亀、合わせて 100 匹いる。足の合計が 272 本のとき、鶴と亀はそれぞれ何匹ずついるか。
  2. 鶴と亀とトンボが合わせて 10 匹いる。足の合計が 38 本で羽の合計が 14 枚であるとき、鶴と亀とトンボはそれぞれ何匹ずついるか。(トンボの足は 6 本で羽は 4 枚)

問題 1 は次の連立方程式を解けば求めることができます。

x + y = 100
2x + 4y = 272
julia> inv([1 1; 2 4]) * [100, 272]
2-element Array{Float64,1}:
 64.0
 36.0

問題 2 は次の連立方程式を解けば求めることができます。

 x +  y +  z = 10
2x + 4y + 6z = 38
2x +      4z = 14
julia> inv([1 1 1; 2 4 6; 2 0 4]) * [10, 38, 14]
3-element Array{Float64,1}:
 3.0
 5.0
 2.0

●フィボナッチ数

フィボナッチ数 Fn の漸化式を行列で表すと次のようになります。

F0 = 0
F1 = 1
[Fn Fn-1] = [1 1; 1 0] * [Fn-1 Fn-2]

したがって、Fn を求める問題は行列 [1 1; 1 0] の n 乗を求める問題に帰着します。xn は log2 n 程度の手間で求めることが可能です。Fn を繰り返しで求めると n に比例する時間がかかるので、n が大きくなると行列を使ったほうが速くなると思われます。

リスト : フィボナッチ数

function fibo(n)
    if n == 0
        0
    elseif n == 1
        1
    else
        f = ([1 1; 1 0] ^ (n - 1)) * [1, 0]
        f[1]
    end
end

for x in 0:40
    @printf "%d " fibo(x)
end
0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 
28657 46368 75025 121393 196418 317811 514229 832040 1346269 2178309 3524578 
5702887 9227465 14930352 24157817 39088169 63245986 102334155

フィボナッチ数列の最初の 2 項を 2, 1 に置き換えた数列の項を「リュカ数 (Lucas number)」といいます。

リスト : リュカ数

function lucas(n)
    if n == 0
        2
    elseif n == 1
        1
    else
        f = ([1 1; 1 0] ^ (n - 1)) * [1, 2]
        f[1]
    end
end

for x in 0:40
    @printf "%d " lucas(x)
end
2 1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364 2207 3571 5778 9349 15127 24476 
39603 64079 103682 167761 271443 439204 710647 1149851 1860498 3010349 4870847 
7881196 12752043 20633239 33385282 54018521 87403803 141422324 228826127

次の漸化式で生成される数列をトリボナッチ数列といいます。

T0 = T1 = 0, T2 = 1
Tn+3 = Tn+2 + Tn+1 + Tn
リスト : トリボナッチ数

function tribo(n)
    if n == 0 || n == 1
        0
    elseif n == 2
        1
    else
        t = ([1 1 1; 1 0 0; 0 1 0] ^ (n - 2)) * [1, 0, 0]
        t[1]
    end
end

for x in 0:30
    @printf "%d " tribo(x)
end
0 0 1 1 2 4 7 13 24 44 81 149 274 504 927 1705 3136 5768 10609 19513 35890 66012
 121415 223317 410744 755476 1389537 2555757 4700770 8646064 15902591

次の漸化式で生成される数列をテトラナッチ数列といいます。

T0 = T1 = T2 = 0, T3 = 1
Tn+4 = Tn+3 + Tn+2 + Tn+1 + Tn
リスト : テトラナッチ数列

function tetra(n)
    if n == 0 || n == 1 || n == 2
        0
    elseif n == 3
        1
    else
        t = ([1 1 1 1; 1 0 0 0 ; 0 1 0 0; 0 0 1 0] ^ (n - 3)) * [1, 0, 0, 0]
        t[1]
    end
end

for x in 0:30
    @printf "%d " tetra(x)
end
0 0 0 1 1 2 4 8 15 29 56 108 208 401 773 1490 2872 5536 10671 20569 39648 76424
147312 283953 547337 1055026 2033628 3919944 7555935 14564533 28074040

●たらいまわし

リスト : たらいまわし関数

function tarai(x, y, z)
  if x <= y
    y
  else
    tarai(tarai(x - 1, y, z), tarai(y - 1, z, x), tarai(z - 1, x, y))
  end
end

function tak(x, y, z)
  if x <= y
    z
  else
    tak(tak(x - 1, y, z), tak(y - 1, z, x), tak(z - 1, x, y))
  end
end

@time print(tarai(14, 7, 0))
@time print(tak(22, 11, 0))

関数 tarai() や tak() は「たらいまわし関数」といい、再帰的に定義されています。これらの関数は、引数の与え方によっては実行に時間がかかるため、Lisp などのベンチマークに利用されることがあります。tarai() は通称「竹内関数」と呼ばれていて、日本の代表的な Lisper である竹内郁雄氏によって考案されたそうです。そして、tak() は tarai() のバリエーションで、John Macarthy 氏によって作成されたそうです。

それでは、さっそく実行してみましょう。

14  2.209532 seconds (12.30 k allocations: 428.349 KB)
11  2.283424 seconds (823 allocations: 28.012 KB)

実行環境 : Windows 7, Core i7-2670QM 2.20GHz

このように、たらいまわし関数は引数の値が小さくても実行に時間がかかります。

●遅延評価による高速化

tarai() は「遅延評価」を行う処理系、たとえば関数型言語の Haskell では高速に実行することができます。また、Scheme でも delay と force を使って遅延評価を行うことができます。tarai() のプログラムを見てください。x <= y のときに y を返しますが、このとき引数 z の値は必要ありませんね。引数 z の値は x > y のときに計算するようにすれば、無駄な計算を省略することができます。なお、tak() は x <= y のときに z を返しているため、遅延評価で高速化することはできません。ご注意ください。

完全ではありませんが、Julia でもクロージャを使って遅延評価を行うことができます。

リスト : 遅延評価

function tarai_lazy(x, y, z)
  if x <= y
    y
  else
    zz = z()
    tarai_lazy(tarai_lazy(x - 1, y, () -> zz),
               tarai_lazy(y - 1, zz, () -> x),
               () -> tarai_lazy(zz - 1, x, () -> y))
  end
end

@time print(tarai(14, 7, 0))
@time print(tarai_lazy(14, 7, () -> 0))

遅延評価したい処理をクロージャに包んで引数 z に渡します。そして、x > y のときに引数 z を評価 (関数呼び出し) します。すると、クロージャ内の処理が実行されて z の値を求めることができます。たとえば、() -> 0 を z に渡す場合、z() とすると返り値は 0 になります。() -> x を渡せば、x に格納されている値が返されます。() -> tarai( ... ) を渡せば、tarai() が実行されてその値が返されるわけです。

実行結果は次のようになりました。

14  2.206594 seconds (12.30 k allocations: 428.349 KB)
14  0.007032 seconds (3.54 k allocations: 122.750 KB)

●メモ化による高速化

たらいまわし関数が遅いのは、同じ値を何度も計算しているためです。この場合、表 (table) を使って処理を高速化することができます。同じ値を何度も計算することがないように、計算した値は表に格納しておいて、2 回目以降は表から計算結果を求めるようにします。このような手法を「表計算法」とか「メモ化」といいます。

Julia の場合、辞書を使うと簡単です。次のリストを見てください。

リスト : メモ化

tak_table = Dict()

function tak_memo(x, y, z)
    key = (x, y, z)
    if !haskey(tak_table, key)
        if x <= y
            tak_table[key] = z
        else
            tak_table[key] = tak_memo(tak_memo(x - 1, y, z),
                                      tak_memo(y - 1, z, x),
                                      tak_memo(z - 1, x, y))
        end
    end
    tak_table[key]
end

@time print(tak(22, 11, 0))
@time print(tak_memo(22, 11, 0))

tak_memo() の値を格納する辞書を大域変数 tak_table に用意します。tak_memo() では、引数 x, y, z を要素とするタプルを作り、それをキーとして辞書 tak_table を検索します。tak_table に key があれば、その値を返します。そうでなければ、値を計算して tak_table にセットして、その値を返します。

実行結果は次のようになりました。

11  2.426078 seconds (12.30 k allocations: 428.341 KB)
11  0.112966 seconds (101.83 k allocations: 2.945 MB, 1.86% gc time)

メモ化により実行速度は速くなりましたが、遅延評価のような劇的な効果はありませんでした。Julia の場合、辞書のアクセスに時間がかかるのかもしれません。興味のある方はいろいろ試してみてください。

●メモ化関数

ところで、プログラミング言語によっては、関数をメモ化する「メモ化関数 (memoize)」を定義することができます。たとえば Python の場合、メモ化関数を memoize() とすると、次のように使います。

tak = memoize(tak)

Python の場合、tak の値を書き換えないと、tak() の中で再帰呼び出しするとき、メモ化した関数を呼び出すことはできません。Julia の場合、メモ化関数は簡単に定義することができますが、funciton などで定義した関数名は定数として扱われるため、値を書き換えることができません。次の例を見てください。

julia> foo(a, b) = a + b
foo (generic function with 1 method)

julia> foo = (a, b) -> a * b
ERROR: invalid redefinition of constant foo

julia> foo = 10
ERROR: invalid redefinition of constant foo

これでは memoize() で再帰関数のメモ化を行うことができません。そこで、今回はちょっと裏技みたいな方法を使ってみましょう。julia の場合、再帰関数は無名関数でも定義することが可能です。次の例を見てください。

julia> fact = n -> if n == 0
       1
       else
       n * fact(n - 1)
       end
(anonymous function)

julia> fact(10)
3628800

julia> fact = 10
10

この場合、fact の値を書き換えることができるので、memoize() でメモ化することができます。memoize() とたらいまわし関数は次のようになります。

リスト : たらいまわし関数

# メモ化関数
function memoize(f)
    table = Dict()
    function func(args...)
        if !haskey(table, args)
            table[args] = f(args...)
        end
        table[args]
    end
    func
end

# たらいまわし関数
tak1 = (x, y, z) ->
    if x <= y
        z
    else
        tak1(tak1(x - 1, y, z), tak1(y - 1, z, x), tak1(z - 1, x, y))
    end

# メモ化
tak1 = memoize(tak1)

@time print(tak(22, 11, 0))
@time print(tak1(22, 11, 0))
11  2.471494 seconds (12.30 k allocations: 428.231 KB)
11  0.131045 seconds (127.11 k allocations: 3.585 MB, 1.47% gc time)

このように、無名関数を使うと再帰関数でも簡単にメモ化することができますが、実はこの方法には重大な欠点があります。tak1 をメモ化しないで実行すると、極端に遅くなるのです。実際に試してみると tak1(22, 11, 0) で約 32 秒もかかってしまいました。Julia では、このような無名関数の使い方は邪道なのかもしれません。

なお、Julia でメモ化関数を実現する場合、マクロを使う方法があるようです。実際に、GitHub で Memoize.jl が公開されています。興味のある方は試してみてください。

●双方向リスト

双方向リストの詳しい説明は Algorithms with Python 連結リストとキュー をお読みください。

#
# dlist.jl : 双方向リスト
#
#            Copyright (C) 2016 Makoto Hiroi
#
module DList

import Base: push!, pop!, shift!, unshift!, length, isempty, show
export Stack, Queue, Deque, front, back

# 双方向リスト (セル)
type Cell
    item
    prev::Cell
    next::Cell
    Cell(a) = (x = new(a); x.next = x; x.prev = x)
end

# 表示
function show(io::IO, head::Cell)
    xs = head.next
    print(io, "(")
    while xs !== head
        print(io, xs.item)
        xs = xs.next
        if xs !== head; print(io, " "); end
    end
    print(io, ")")
end

# 空か?
isempty(xs::Cell) = xs.next === xs

# 長さ
function length(head::Cell)
    c = 0
    xs = head.next
    while xs !== head
        c += 1
        xs = xs.next
    end
    c
end

# 参照
function front(head::Cell)
    if isempty(head)
        error("front: Empty DList")
    end
    head.next.item
end

function back(head::Cell)
    if isempty(head)
        error("back: Empty DList")
    end
    head.prev.item
end

# 先頭に挿入
function unshift!(xs::Cell, args...)
    for x in args
        ys = Cell(x)
        zs = xs.next
        ys.prev = xs
        ys.next = zs
        zs.prev = ys
        xs.next = ys
    end
    xs
end

# 先頭要素を取り出す
function shift!(xs::Cell)
    if isempty(xs)
        error("shift!: Empty DList")
    end
    ys = xs.next
    zs = ys.next
    xs.next = zs
    zs.prev = xs
    ys.item
end

# 末尾に挿入
function push!(xs::Cell, args...)
    for x in args
        ys = Cell(x)
        zs = xs.prev
        ys.next = xs
        ys.prev = zs
        xs.prev = ys
        zs.next = ys
    end
    xs
end

# 末尾要素を取り出す
function pop!(xs::Cell)
    if isempty(xs)
        error("shift!: Empty DList")
    end
    ys = xs.prev
    zs = ys.prev
    xs.prev = zs
    zs.next = xs
    ys.item
end

# スタック
type Stack
    top::Cell
    Stack() = new(Cell(nothing))
end

isempty(xs::Stack) = isempty(xs.top)
length(xs::Stack) = length(xs.top)
show(io::IO, xs::Stack) = show(io, xs.top)
front(xs::Stack) = front(xs.top)
push!(xs::Stack, args...) = push!(xs.top, args...)
pop!(xs::Stack) = pop!(xs.top)

# キュー
type Queue
    head::Cell
    Queue() = new(Cell(nothing))
end

isempty(xs::Queue) = isempty(xs.head)
length(xs::Queue) = length(xs.head)
show(io::IO, xs::Queue) = show(io, xs.head)
front(xs::Queue) = front(xs.head)
push!(xs::Queue, args...) = push!(xs.head, args...)
shift!(xs::Queue) = shift!(xs.head)

# ディーキュー
type Deque
    head::Cell
    Deque() = new(Cell(nothing))
end

isempty(xs::Deque) = isempty(xs.head)
length(xs::Deque) = length(xs.head)
show(io::IO, xs::Deque) = show(io, xs.head)
front(xs::Deque) = front(xs.head)
back(xs::Deque) = back(xs.head)
push!(xs::Deque, args...) = push!(xs.head, args...)
pop!(xs::Deque) = pop!(xs.head)
shift!(xs::Deque) = shift!(xs.head)
unshift!(xs::Deque, args...) = unshift!(xs.head, args...)

end
julia> using DList

julia> a = Stack()
()

julia> push!(a,1,2,3,4,5)
(1 2 3 4 5)

julia> length(a)
5

julia> while !isempty(a); println(pop!(a)); end
5
4
3
2
1

julia> b = Queue()
()

julia> push!(b, 1,2,3,4,5)
(1 2 3 4 5)

julia> length(b)
5

julia> while !isempty(b); println(shift!(b)); end
1
2
3
4
5

julia> c = Deque()
()

julia> push!(c, 1,2,3,4,5)
(1 2 3 4 5)

julia> front(c)
1

julia> back(c)
5

julia> unshift!(c, 1,2,3,4,5)
(5 4 3 2 1 1 2 3 4 5)

julia> pop!(c)
5

julia> c
(5 4 3 2 1 1 2 3 4)

julia> shift!(c)
5

julia> c
(4 3 2 1 1 2 3 4)

●REPL (Read - Eval - Print - Loop)

リスト : Read - Eval - Print - Loop

while true
    print(">>> ")
    s = readline(STDIN)
    try
        println(eval(parse(s)))
    catch e
        println(e)
    end
end
>>> 1 + 2 * 3
7
>>> a = 10
10
>>> a * 10
100
>>> b
UndefVarError(:b)
>>> square(x) = x * x
square
>>> square(1.2345)
1.5239902499999998
>>> exit()

C>

●遅延評価

Scheme の遅延評価は delay と force を使いますが、Julia でもマクロを使って簡単に実装することができます。次のリストを見てください。

リスト : 遅延評価

# プロミス
type Promise
    func
    flag::Bool
    value
end

macro delay(expr)
    :(Promise(() -> $expr, false, nothing))
end

function force(p::Promise)
    if !p.flag
        p.value = p.func()
        p.flag = true
    end
    p.value
end

delay() はマクロで、引数 expr を評価しないでプロミス (Promise) というデータを返します。expr はこのプロミスに保存されていて、force(p) を実行すると、expr を評価してその値を返します。このとき、値がプロミスに保存されることに注意してください。再度 force(p) を実行すると、保存された値が返されます。

プログラムは簡単です。Promis のメンバ変数 func には、評価する式 expr を包んだ無名関数をセットします。格納した式を評価したら、flag を true にセットします。value は式の評価結果を格納します。

delay() は簡単で、引数 expr を無名関数に包んで Promise() に渡します。これで、引数 expr は評価されずにプロミスに格納されます。force() も簡単で、p.flag が false ならば、p.func() を評価して値を求め、その結果を p.value にセットします。あとは、p.value を返すだけです。

簡単な実行例を示します。

julia> include("lazy.jl")

julia> a = @delay (println("oops!"); 10 + 20)
Promise((anonymous function),false,nothing)

julia> force(a)
oops!
30

julia> force(a)
30

delay() の返り値を変数 a にセットします。このとき、式 (println("oops!"); 10 + 20) は評価されていません。force(a) を評価すると、式を評価して値 30 を返します。このとき、println() で oops! と表示されます。また、force(a) を再度実行すると、同じ式を再評価することなく値を求めることができます。この場合は oops! と表示されません。

●遅延ストリーム

「ストリーム (stream)」はデータの流れを抽象化したデータ構造です。たとえば、ファイル入出力はストリームと考えることができます。また、1 次元配列 (ベクタ) を使ってストリームを表すこともできます。ただし、単純なベクタでは有限個のデータの流れしか表すことができません。ところが、遅延評価を用いると擬似的に無限個のデータを表すことができるようになります。これを「遅延ストリーム」と呼びます。

遅延ストリームの基本的な考え方は、必要になったときに新しいデータを生成することです。このときに遅延評価を用います。具体的にはデータを生成する関数を用意し、それを遅延評価してストリームに格納しておきます。そして、必要になった時点で遅延評価しておいた関数を呼び出して値を求めればよいわけです。

遅延ストリームを表すデータ型は次のようになります。

リスト : 遅延ストリームのデータ型

# 終端
immutable Nil
end

# 遅延ストリーム
immutable Cons
    head
    tail::Promise
end

# アクセス関数
stream_head(s::Cons) = s.head
stream_tail(s::Cons) = force(s.tail)

データ型は Cons としました。Nil はストリームの終端を表します。無限ストリームだけを扱うのであれば Nil は必要ありません。Cons の head が現時点での先頭データを表し、tail が遅延ストリームを生成する関数を格納するプロミスになります。

関数 stream_head() は遅延ストリームの先頭データ head を返します。関数 stream_tail() は tail を force() して次の要素を格納した遅延ストリームを返します。

●遅延ストリームの生成

それでは、遅延ストリームを生成する関数を作りましょう。たとえば、low から high までの整数列を生成するストリームは次のようにプログラムすることができます。

リスト : 整数列を生成するストリーム

function intgen(low, high)
    if low > high
        Nil
    else
        Cons(low, @delay intgen(low + 1, high))
    end
end

関数 intgen は遅延ストリームを生成して返します。Cons の第 1 引数が現時点でのデータになります。第 2 引数のプロミスを force すると、intgen(low + 1, high) が実行されて次のデータを格納した遅延ストリームが返されます。そして、その中の遅延オブジェクトを force すると、その次のデータを得ることができます。

簡単な実行例を示しましょう。

julia> s = intgen(1, 4)
Cons(1,Promise((anonymous function),false,nothing))

julia> stream_head(s)
1

julia> s = stream_tail(s)
Cons(2,Promise((anonymous function),false,nothing))

julia> stream_head(s)
2

julia> s = stream_tail(s)
Cons(3,Promise((anonymous function),false,nothing))

julia> stream_head(s)
3

julia> s = stream_tail(s)
Cons(4,Promise((anonymous function),false,nothing))

julia> stream_head(s)
4

julia> s = stream_tail(s)
Nil

もう一つ、簡単な例を示しましょう。フィボナッチ数列を生成する遅延ストリームを作ります。次のリストを見てください。

リスト : フィボナッチ数列を生成する遅延ストリーム

fibgen(a, b) = Cons(a, @delay fibgen(b, a + b))

関数 fibgen() の引数 a がフィボナッチ数列の最初の項で、b が次の項です。したがって、プロミスに fibgen(b, a + b) を格納しておけば、force することでフィボナッチ数列を生成することができます。

julia> f = fibgen(0, 1)
Cons(0,Promise((anonymous function),false,nothing))

julia> stream_head(f)
0

julia> f = stream_tail(f)
Cons(1,Promise((anonymous function),false,nothing))

julia> stream_head(f)
1

julia> f = stream_tail(f)
Cons(1,Promise((anonymous function),false,nothing))

julia> stream_head(f)
1

julia> f = stream_tail(f)
Cons(2,Promise((anonymous function),false,nothing))

julia> stream_head(f)
2

julia> f = stream_tail(f)
Cons(3,Promise((anonymous function),false,nothing))

julia> stream_head(f)
3

julia> f = stream_tail(f)
Cons(5,Promise((anonymous function),false,nothing))

julia> stream_head(f)
5

●遅延ストリームの操作関数

次は遅延ストリームを操作する関数を作りましょう。最初は n 番目の要素を求める関数 stream_ref() です。今回は先頭の要素を 1 番目とします。

リスト  : n 番目の要素を求める

function stream_ref(s::Cons, n::Int)
    while n > 1
        if s === Nil
            error("stream_ref: empty stream")
        end
        s = stream_tail(s)
        n -= 1
    end
    stream_head(s)
end

stream_ref() は stream_tail() を n 回繰り返すことで n 番目の要素を求めます。ストリームから n 個の要素を取り出してリストに格納して返す関数 stream_take() も同様にプログラムすることができます。

リスト : n 個の要素を取り出す

function stream_take(s::Cons, n::Int)
    a = []
    while s !== Nil && n > 0
        push!(a, stream_head(s))
        s = stream_tail(s)
        n -= 1
    end
    a
end

それでは、簡単な実行例を示しましょう。

julia> f = fibgen(big(0), big(1))
Cons(0,Promise((anonymous function),false,nothing))

julia> for x in 1:20
       print(stream_ref(f, x)); print(" ")
       end
0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181

julia> stream_take(f, 20)
20-element Array{Any,1}:
    0
    1
    1
    2
    3
    5
    8
   13
   21
   34
   55
   89
  144
  233
  377
  610
  987
 1597
 2584
 4181

julia> stream_ref(f, 100)
218922995834555169026

変数 f にフィボナッチ数列を生成するストリームをセットします。stream_ref() で順番に要素を 20 個取り出すと、その値はフィボナッチ数列になっていますね。同様に、stream_take() で 20 個の要素を取り出すと、配列の要素はフィボナッチ数列になります。初期値を多倍長整数 (BigInt) で指定しているので、メモリの許す限り大きなフィボナッチ数でも求めることができます。

このほかにも、便利な関数や高階関数などをいろいろ定義することができます。遅延ストリームに興味のある方は拙作のページ お気楽 Scheme プログラミング入門 遅延ストリーム (1), (2) をお読みください。


Copyright (C) 2016 Makoto Hiroi
All rights reserved.

[ Home | Light | Julia ]