M.Hiroi's Home Page

Julia Language Programming

お気楽 Julia プログラミング超入門


Copyright (C) 2018-2021 Makoto Hiroi
All rights reserved.

Julia の基礎知識

●高階関数の使い方

julia> map(x -> x * x, 1 : 10)
10-element Vector{Int64}:
   1
   4
   9
  16
  25
  36
  49
  64
  81
 100

julia> a = zeros(Int, 10)
10-element Vector{Int64}:
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0

julia> map!(x -> x * x, a, 1 : 10)
10-element Vector{Int64}:
   1
   4
   9
  16
  25
  36
  49
  64
  81
 100

julia> a
10-element Vector{Int64}:
   1
   4
   9
  16
  25
  36
  49
  64
  81
 100

julia> filter(iseven, a)
5-element Vector{Int64}:
   4
  16
  36
  64
 100

julia> filter!(iseven, a)
5-element Vector{Int64}:
   4
  16
  36
  64
 100

julia> a
5-element Vector{Int64}:
   4
  16
  36
  64
 100

julia> d = Dict("foo" => 1, "bar" => 2, "baz" => 3, "oops" => 4)
Dict{String,Int64} with 4 entries:
  "bar"  => 2
  "baz"  => 3
  "oops" => 4
  "foo"  => 1

julia> filter(p -> iseven(p.second), d)
Dict{String,Int64} with 2 entries:
  "bar"  => 2
  "oops" => 4

julia> filter!(p -> iseven(p.second), d)
Dict{String,Int64} with 2 entries:
  "bar"  => 2
  "oops" => 4

julia> d
Dict{String,Int64} with 2 entries:
  "bar"  => 2
  "oops" => 4

julia> reduce((a, x) -> (x, a), 1 : 10)
(10, (9, (8, (7, (6, (5, (4, (3, (2, 1)))))))))

julia> reduce((a, x) -> (x, a), 1 : 10, init=0)
(10, (9, (8, (7, (6, (5, (4, (3, (2, (1, 0))))))))))

julia> foldl((a, x) -> (x, a), 1 : 10)
(10, (9, (8, (7, (6, (5, (4, (3, (2, 1)))))))))

julia> foldl((a, x) -> (x, a), 1 : 10, init=0)
(10, (9, (8, (7, (6, (5, (4, (3, (2, (1, 0))))))))))

julia> foldr((x, a) -> (x, a), 1 : 10)
(1, (2, (3, (4, (5, (6, (7, (8, (9, 10)))))))))

julia> foldr((x, a) -> (x, a), 1 : 10, init=11)
(1, (2, (3, (4, (5, (6, (7, (8, (9, (10, 11))))))))))

julia> m = reshape(collect(1 : 16), 4, 4)
4×4 Matrix{Int64}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> reduce(+, m, dims=1)
1×4 Matrix{Int64}:
 10  26  42  58

julia> reduce(+, m, dims=2)
4×1 Matrix{Int64}:
 28
 32
 36
 40
julia> all(iseven, [2,4,6,8,10])
true

julia> all(iseven, [2,4,6,8,11])
false

julia> any(isodd, [2,4,6,8,11])
true

julia> any(isodd, [2,4,6,8,10])
false

julia> all([true, true, true])
true

julia> all([true, true, false])
false

julia> any([true, true, false])
true

julia> any([false, false, false])
false

julia> count(iseven, 1 : 10)
5

julia> count(iseven, 1 : 2 : 10)
0

julia> count(map(iseven, 1 : 10))
5

julia> a = [1, 2, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 5]
14-element Vector{Int64}:
 1
 2
 1
 2
 3
 1
 2
 3
 4
 1
 2
 3
 4
 5

julia> findfirst(isequal(4), a)
9

julia> findnext(isequal(4), a, 10)
13

julia> findnext(isequal(4), a, 14)

julia> findlast(isequal(4), a)
13

julia> findprev(isequal(4), a, 12)
9

julia> findprev(isequal(4), a, 8)

julia> findall(isequal(4), a)
2-element Vector{Int64}:
  9
 13

julia> findall(isequal(10), a)
Int64[]

julia> foreach(x -> print("$x "), a)
1 2 1 2 3 1 2 3 4 1 2 3 4 5

簡単なプログラム

●約数の個数

自然数 n の約数の個数を求めるプログラムを作りましょう。n の素因数分解ができると、約数の個数を求めるのは簡単です。\(n = p^{a} \times q^{b} \times r^{c}\) とすると、約数の個数は \((a + 1)(b + 1)(c + 1)\) になります。たとえば、12 は \(2^{2} \times 3^{1}\) になるので、約数の個数は 3 * 2 = 6 になります。実際、12 の約数は 1, 2, 3, 4, 6, 12 の 6 個です。

拙作のページ「素因数分解」で作成した関数 factorization() を使うと、プログラムは次のようになります。

リスト : 約数の個数

divisor_num(n) = foldl((a, p) -> a * (p.second + 1), factorization(n), init=1)

関数 divisor_num() は高階関数 foldl() で配列の要素 (Pair) を順番に取り出し、x + 1 を a に掛け算していくだけです。Pair の最初の要素は first で、二番目の要素は second でアクセスすることができます。

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

julia> divisor_num(24)
8

julia> divisor_num(12345678)
24

julia> divisor_num(123456789)
12

julia> divisor_num(1234567890)
48

julia> divisor_num(1111111111)
16

●約数の和

n の素因数分解ができると、約数の合計値を求めるのは簡単です。素因数分解した結果が \(p^a\) の場合、その約数は \(1, p, p^2, \cdots, p^a\) となり、初項が 1 で公比が p の等比数列になります。約数の合計値を \(s(p, a)\) で表すことにすると、等比数列の和の公式から \(s(p, a)\) は以下の式になります。

\( s(p, a) = p^a + p^{a-1} + \cdots + p^2 + p + 1 = \dfrac{p^{a+1} - 1}{p - 1} \)

たとえば、8 の素因数分解は \(2^3\) になり、素数の合計値は 8 + 4 + 2 + 1 = 15 になります。整数 n の素因数分解が \(p^a \times q^b \times r^c\) だとすると、\(s(p, a) \times s(q, b) \times s(r, c)\) は約数の合計値になります。実際に式を展開すると次のようになります。

\(\begin{array}{ll} & (p^{a} + p^{a-1} + \cdots + p^{2} + p + 1) \\ \times & (q^{b} + q^{b-1} + \cdots + q^{2} + q + 1) \\ \times & (r^{c} + r^{c-1} + \cdots + r^{2} + r + 1) \\ = & p^{a} * q^{b} * r^{c} + \cdots + p + q + r + 1 \end{array}\)

このように、各項は約数と一対一に対応しているので、その和は約数の合計値になるわけです。たとえば、12 は \(2^2 \times 3\) に素因数分解できますが、その合計値は (4 + 2 + 1) * (3 + 1) = 28 となります。12 の約数は 1, 2, 3, 4, 6, 12 なので、その合計値は確かに 28 になります。

拙作のページ「素因数分解」で作成した関数 factorization() を使うと、プログラムは次のようになります。

リスト : 約数の合計値

function divisor_sum(n)
    sigma(p, n) = div(p ^ (n + 1) - 1, p - 1)    # s(p, n) の計算
    foldl((a, p) -> a * sigma(p...), factorization(n), init = 1)
end

局所関数 sigma() は \(s(p, n)\) を計算します。あとは foldl() で sigma() の返り値を累積変数 a に掛け算していくだけです。

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

julia> divisor_sum(24)
60

julia> divisor_sum(12345678)
27319968

julia> divisor_sum(123456789)
178422816

julia> divisor_sum(1234567890)
3211610688

julia> divisor_sum(1111111111)
1246404096

●約数

p が素数の場合、pa の約数は次のように簡単に求めることができます。

\( p^{a}, \ p^{a-1}, \ \cdots, \ p^{2}, \ p, \ 1 \)

n の素因数分解が \(p^a \times q^b\) だったとすると、その約数は次のようになります。

\(\begin{array}{l} (p^{a}, \ p^{a-1}, \ \cdots, \ p^{2}, \ p, \ 1) \times q^b, \\ (p^{a}, \ p^{a-1}, \ \cdots, \ p^{2}, \ p, \ 1) \times q^{b-1}, \\ \quad \quad \cdots \cdots \\ (p^{a}, \ p^{a-1}, \ \cdots, \ p^{2}, \ p, \ 1) \times q^2, \\ (p^{a}, \ p^{a-1}, \ \cdots, \ p^{2}, \ p, \ 1) \times q^1, \\ (p^{a}, \ p^{a-1}, \ \cdots, \ p^{2}, \ p, \ 1) \times 1 \end{array}\)

たとえば、12 の約数は \(2^2\) = (1, 2, 4) と 3 = (1, 3) から、(1, 2, 4) * 1 と (1, 2, 4) * 3 のすべての要素 (1, 2, 4, 3, 6, 12) になります。

プログラムは次のようになります。

リスト : 約数をすべて求める

function divisor(n)
    # pq の約数を求める
    divisor_sub(p, q) = [p ^ i for i = 0 : q]

    xs = factorization(n)
    ys = divisor_sub(xs[1]...)
    for p = xs[2 : end]
        ys = [x * y for x = divisor_sub(p...) for y = ys]
    end
    sort(ys)
end

局所関数 divisor_sub() は \(p^q\) の約数を配列に格納して返します。引数 n を factorization() で素因数分解して変数 xs にセットします。xs の先頭要素を divisor_sub() に渡して配列に変換し、それを変数 ys にセットします。あとは for ループで xs の 1 番目から要素を順番に取り出し、\(p^q\) を divisor_sub() でリストに変換して、それを内包表記で累積変数 ys のリストの要素と掛け合わせていくだけです。

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

julia> foreach(x -> print("$x "), divisor(24))
1 2 3 4 6 8 12 24
julia> foreach(x -> print("$x "), divisor(123456789))
1 3 9 3607 3803 10821 11409 32463 34227 13717421 41152263 123456789
julia> foreach(x -> print("$x "), divisor(1234567890))
1 2 3 5 6 9 10 15 18 30 45 90 3607 3803 7214 7606 10821 11409 18035 19015 21642 
22818 32463 34227 36070 38030 54105 57045 64926 68454 108210 114090 162315 171135 
324630 342270 13717421 27434842 41152263 68587105 82304526 123456789 137174210 
205761315 246913578 411522630 617283945 1234567890
julia> foreach(x -> print("$x "), divisor(1111111111))
1 11 41 271 451 2981 9091 11111 100001 122221 372731 2463661 4100041 27100271 
101010101 1111111111

●完全数

完全数 - Wikipedia によると、『完全数(かんぜんすう,perfect number)とは、その数自身を除く約数の和が、その数自身と等しい自然数のことである。』 とのことです。完全数を求めるプログラムを作りましょう。

リスト : 完全数

function perfect_number(n)
    [x for x = 2 : n if divisor_sum(x) - x == x]
end

完全数を求める perfect_number() は簡単です。x の約数の合計値を divisor_sum() で求め、その値から x を引いた値が x と等しければ完全数です。内包表記で完全数を格納した配列を生成して返します。

実行結果を示します。

julia> perfect_number(10000)
4-element Vector{Int64}:
    6
   28
  496
 8128

ところで、完全数 - Wikipedia によると、メルセンヌ素数を Mn とすると、偶数の完全数は 2n-1 * Mn で表すことができるそうです。この式を使うと偶数の完全数は次のようになります。

 n : メルセンヌ素数 : 完全数
---+----------------+----------------------
 2 : 3              : 6
 3 : 7              : 28
 5 : 31             : 496
 7 : 127            : 8128
13 : 8191           : 33550336
17 : 131071         : 8589869056
19 : 524287         : 137438691328
31 : 2147483647     : 2305843008139952128

なお、奇数の完全数はまだ発見されておらず、偶数の完全数 (つまりメルセンヌ素数) が無数に存在するか否かも未解決な問題だそうです。


●友愛数

友愛数 - Wikipedia によると、『友愛数(ゆうあいすう)とは、異なる2つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。』 とのことです。友愛数を求めるプログラムを作りましょう。

リスト : 友愛数

function yuuai_number(n)
    check(x, m) = m < x && x == divisor_sum(m) - m
    [p for p = [(x, divisor_sum(x) - x) for x = 2 : n] if check(p...)]
end

友愛数を求める関数 yuuai_number() も簡単です。局所関数 check(x, m) は x と m が友愛数かチェックします。m の約数の合計値から m を引いた値が x と等しければ、x と m は友愛数です。同じ組を表示しないようにするため、m < x を条件に入れています。あとは内包表記でタプル (x, m) を格納した配列を生成し、その要素が友愛数か check() でチェックするだけです。この処理も内包表記を使えば簡単です。

実行結果を示します。

julia> yuuai_number(100000)
13-element Vector{Tuple{Int64, Int64}}:
 (284, 220)
 (1210, 1184)
 (2924, 2620)
 (5564, 5020)
 (6368, 6232)
 (10856, 10744)
 (14595, 12285)
 (18416, 17296)
 (66992, 66928)
 (71145, 67095)
 (76084, 63020)
 (87633, 69615)
 (88730, 79750)

なお、友愛数が無数に存在するか否かは、未解決な問題だそうです。


●別解

perfect_number() と yuuai_number() は、divisor_sum() を呼び出して約数の和を求めていますが、あらかじめ約数の和を計算して配列に格納しておく方法もあります。この場合、約数の和の計算にちょっと時間がかかりますが、完全数と友愛数を求める処理は高速になります。

基本的な考え方は簡単です。約数の和を格納する配列 table を用意します。table の要素は 1 に初期化します。2 から順番に素数 p で割り算して 1 + p1 + ... + pq を求め、それを table の値に掛け算します。素数は「エラトステネスの篩」と同じ方法で求めることができます。素数の倍数であれば、table の値は 1 よりも大きくなっているはずです。つまり table が 1 であれば、その整数は素数であることがわかります。

プログラムは次のようになります。

リスト : n 以下の整数の約数の和を配列に格納して返す

function make_divisor_sum(n)
    function factor_sub(n, p)
        a = 1
        q = 1
        while n % p == 0
            a += p ^ q
            q += 1
            n = div(n, p)
        end
        a
    end
    #
    table = fill(1, n)
    for i = 2 : n
        if table[i] != 1 continue end
        for j = i : i : n
            table[j] *= factor_sub(j, i)
        end
    end
    table
end

局所関数 factor_sub() は n を 素数 p で割り算して、1 + p1 + ... + pq を求めます。make_divisor_sum() は table を 1 で初期化してから、for ループで素数を探します。table[i] が 1 でない場合、i は素数ではありません。contiune で次の数をチェックします。素数の場合は、その倍数に対応する table の値を更新します。i の倍数を 変数 j にセットし、table[j] に factor_sub(j, i) の値を乗算するだけです。最後に table を返します。

make_divisor_sum() を使うと yuuai_number() は次のようになります。

リスト : 友愛数

function yuuai_number1(n)
    table = make_divisor_sum(n)
    check(x, m) = m < x && x == table[m] - m
    [p for p = [(x, table[x] - x) for x = 2 : n] if check(p...)]
end

それでは実際に 1,000,000 以下の友愛数を求めてみましょう。

julia> @time yuuai_number(1000000)
  1.240681 seconds (3.43 M allocations: 352.103 MiB, 16.94% gc time)
40-element Vector{Tuple{Int64, Int64}}:
 (284, 220)
 (1210, 1184)
 (2924, 2620)
 (5564, 5020)
 (6368, 6232)
 (10856, 10744)
 (14595, 12285)
 (18416, 17296)
 (66992, 66928)
 (71145, 67095)
 (76084, 63020)
 (87633, 69615)
 (88730, 79750)
 ⋮
 (514736, 503056)
 (525915, 522405)
 (652664, 643336)
 (669688, 600392)
 (686072, 609928)
 (691256, 624184)
 (712216, 635624)
 (783556, 667964)
 (796696, 726104)
 (863835, 802725)
 (901424, 879712)
 (980984, 898216)

julia> @time yuuai_number1(1000000)
  0.094811 seconds (8 allocations: 22.889 MiB, 4.44% gc time)
40-element Vector{Tuple{Int64, Int64}}:

・・・省略・・・

40 個の友愛数を出力するのに最初のプログラムでは 1.7 秒 (実行環境 : Julia ver 1.6.4, Windows 10, Intel Core i5-6200U 2.30GHz) かかりましたが、make_divisor_sum() を使ったプログラムは約数の和を求める処理を含めて 0.13 秒ですみました。約 13 倍高速化することができました。


●過剰数と不足数

過剰数 - Wikipedia, 不足数 - Wikipedia によると、『その数自身を除く約数の総和が元の数より大きい数』を「過剰数 (abundant number)」といい、『その数自身を除く約数の総和が元の数より小さい数』を「不足数 (deficient number)」というそうです。過剰数と不足数の個数を求めるプログラムも簡単に作ることができます。

リスト : 過剰数と不足数

# 過剰数
abundant_number(n) = count(x -> divisor_sum(x) - x > x, 1 : n)

# 不足数
deficient_number(n) = count(x -> divisor_sum(x) - x < x, 1 : n)
julia> abundant_number(1000000)
247545

julia> deficient_number(1000000)
752451

1000000 以下の過剰数は 247545 個、不足数は 752451 個、完全数は 4 個なので、合計で 1000000 になります。過剰数 - Wikipedia によると、『自然数のうち過剰数が占める割合は0.2474から0.2480の間であると証明されている。』 とのことで、1000000 以下の過剰数の個数は確かにこの範囲内に入っています。


番外編 : たらいまわしと遅延評価

●たらいまわし関数とは?

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

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 氏によって作成されたそうです。

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

$ julia tarai.jl
14  1.704773 seconds (2.52 k allocations: 198.102 KiB, 0.72% compilation time)
11  1.955073 seconds (3 allocations: 104 bytes)

実行環境 : Ubunts 22.04 (WSL2), Julia ver 1.10.5, Intel Core i5-6200U 2.30GHz

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

●遅延評価 (クロージャ) による高速化

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))
@time print(tarai_lazy(14, 7, () -> 0))
@time print(tarai_lazy(14, 7, () -> 0))

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

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

14  1.778538 seconds (2.52 k allocations: 198.102 KiB, 4.71% compilation time)
14  0.054938 seconds (53.82 k allocations: 3.129 MiB, 99.58% compilation time)
14  0.008440 seconds (1.72 k allocations: 115.477 KiB, 97.77% compilation time)
14  0.009211 seconds (1.72 k allocations: 115.477 KiB, 97.85% compilation time)

tarai_lazy() は初回と 2 回目以降だと実行時間が大きく変わりました。JIT の影響かもしれません。

●メモ化による高速化

たらいまわし関数が遅いのは、同じ値を何度も計算しているためです。この場合、表 (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  1.892896 seconds (2.52 k allocations: 198.102 KiB, 2.42% compilation time)
11  0.098352 seconds (58.11 k allocations: 4.053 MiB, 99.05% compilation 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
#5 (generic function with 1 method)

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  1.871978 seconds (2.52 k allocations: 198.102 KiB, 0.62% compilation time)
11  0.122191 seconds (101.61 k allocations: 7.137 MiB, 98.59% compilation time)

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

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

●delay() と force() の作成

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

リスト : 遅延評価 (lazy.jl)

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

macro delay(expr)
    :(Promise(() -> $(esc(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 を返すだけです。

簡単な実行例を示します。上記プログラムを格納したファイル (lazy.jl) を include してください。

julia> a = @delay (println("oops!"); 10 + 20)
Promise(var"#9#10"(), 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! と表示されません。

次は、たらいまわし関数で試してみましょう。

リスト : たらいまわし (遅延評価)

function tarai(x, y, z)
  if x <= y
    y
  else
    zz = force(z)
    tarai(tarai(x - 1, y, @delay zz),
          tarai(y - 1, zz, @delay x),
          @delay tarai(zz - 1, x, @delay y))
  end
end

遅延評価したい処理をプロミスにして引数 z に渡します。そして、x > y のときに引数 z のプロミスを force で評価します。すると、プロミス内の処理が評価されて z の値を求めることができます。たとえば、@delay 0 を z に渡す場合、force(z) とすると返り値は 0 になります。@delay x を渡せば、x に格納されている値が返されます。@delay tarai( ... ) を渡せば tarai() が実行されて、その値を求めることができます。

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

julia> @time tarai(14, 7, @delay 0)
  0.029658 seconds (11.16 k allocations: 746.266 KiB, 99.41% compilation time)
14

julia> @time tarai(14, 7, @delay 0)
  0.000274 seconds (625 allocations: 32.156 KiB, 68.83% compilation time)
14

julia> @time tarai(14, 7, @delay 0)
  0.000169 seconds (625 allocations: 32.156 KiB, 55.05% compilation time)
14

julia> @time tarai(14, 7, @delay 0)
  0.000243 seconds (625 allocations: 32.156 KiB, 54.06% compilation time)
14

実行環境 : Ubunts 22.04 (WSL2), Julia ver 1.10.5, Intel Core i5-6200U 2.30GHz

クロージャを使った遅延評価よりも速くなりました。評価結果をキャッシュしているのが効いているのかもしれません。興味のある方はいろいろ試してみてください。


初版 2018 年 10 月 27 日
改訂 2021 年 11 月 27 日