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

Julia Language Programming

Julia で遊ぼう

[ Home | Light | Julia ]

数で遊ぼう

今回は簡単な数理パズルを出題します。プログラムを作って解いてもかまいませんが、なかには筆算 (電卓) で解くことができる問題もあるので、興味のある方は挑戦してみてください。

  1. 1000000 以下の自然数で、3 の倍数になっている数字の和を求めてください。
  2. 10000! の末尾に付く 0 の個数を求めてください。
  3. 7654321 の末尾の数字を求めてください。
  4. 将棋盤の1ずつのマスに米粒を置きます。最初のマスへは1粒、次のマスへは2粒、そのつぎのマスへは4粒というようにして、つぎつぎに倍増していきます。最後のマス (81 マス) まで置き終わったときの米粒の総数を求めてください。
  5. 7 以上の素数で割り切れない N 以下の正の整数を求めるプログラムを作ってください (ハミングの問題)。
-- 参考文献 --------
1. Steven G. krantz (著), 関沢正躬 (訳), 『問題解決への数学』, 丸善, 2001
2. 中村義作, 『どこまで解ける日本の算法』, ブルーバックス, 1994
3. 大村平, 『数学公式のはなし』, 日科技連, 1996
4. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991

●解答1

今のパソコンは高性能なので、次のようにプログラムしても瞬時に答えを求めることができます。

julia> function sum_of_multiples(n, m)
       a = zero(n)
       for x in m : m : n
       a += x
       end
       a
       end
sum_of_multiples (generic function with 1 method)

julia> sum_of_multiples(big(1000000), 3)
166666833333

ところが、数列の和を求める公式を使うと、もっと簡単に答えを求めることができます。

1 + 2 + 3 + ... + n = n(n + 1)/ 2

上記公式より n 個の 3 の倍数の和は 3 + 6 + 9 + ... + 3n = 3n(n + 1) / 2 となります。したがって、1000000 以下の 3 の倍数の和は 1 から div(1000000, 3) = 333333 までの和を 3 倍することで求めることができます。

3 * 333333 * (333333 + 1) / 2 = 166666833333

これをプログラムすると次のようになります。

julia> function sum_of_multiples1(n, m)
       x = div(n, m)
       div(m * x * (x + 1), 2)
       end
sum_of_multiples1 (generic function with 1 method)

julia> sum_of_multiples1(big(1000000), 3)
166666833333

●等差数列の和

次のように、一定の差で並んだ数列を「等差数列」といいます。

a, a + d, a + 2d, a + 3d, ..., a + (n - 1)d

a を「初項」、d を「公差」といいます。等差数列の一般項は次の式で表すことができます。

an = a + (n - 1)d

初項から an までの和 Sn は次の式で求めることができます。

Sn = n(2a + (n - 1)d) / 2

初項を 1, 公差 を 1 とすると、1 から n までの和は n(n + 1)/ 2 となります。

この公式は次のように導出することができます。

Sn = a              + (a + d)        + ,,, + (a + (n - 2)d) + (a + (n - 1)d)
Sn = (a + (n - 1)d) + (a + (n - 2)d) + ... + (a + d)        + a

足し算すると

2Sn = (2a + (n - 1)d) + (2a + (n - 1)d) + ... (2a + (n - 1)d) + (2a + (n - 1)d)
2Sn = n(2a + (n - 1)d)
 Sn = n(2a + (n - 1)d)/2

このように、右辺を逆順に並べ替えて足し算すると、2a + (n - 1)d が n 個並ぶことになります。あとは、これを 2 で割り算すればいいわけです。


●解答2

10000! であれば、次のようなプログラムでも瞬時に答えを求めることができます。

julia> a = factorial(big(10000)); c = 0; while a % 10 == 0; c += 1; a = div(a, 10); end; c
2499

単純に n! を求めて、10 で割れる回数を求めているだけです。ところが、この方法では n が大きくなると極端に遅くなります。多倍長整数の場合、除算や余りを求める処理は乗算よりもはるかに時間がかかります。たとえば、1 桁増やした 100000! の場合、階乗の値は短時間で求めることができても、div(a, 10) の回数が増えることにより実行時間が極端に遅くなるのです。

そこで、他の方法を考えてみましょう。階乗を計算するとき、末尾に 0 が付くのは値を 10 倍したときです。これは数字 10 や 100 を乗算するときだけではありません。次の例を見てください。

1 = 1
1 * 2 = 2
1 * 2 * 3 = 6
1 * 2 * 3 * 4 = 24
1 * 2 * 3 * 4 * 5 = 120
1 * 2 * 3 * 4 * 5 * 6 = 720
1 * 2 * 3 * 4 * 5 * 6 * 7 = 5040
1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 = 40320 
1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 = 362880 
1 * 2 * 3 * 4 * 5 * 6 * 7  * 8 * 9 * 10 = 3628800 

10 は 2 * 5 に素因数分解することができます。つまり、2 と 5 の組があれば、末尾に 0 がひとつ追加されるわけです。また、0 が複数追加されることもあります。次の例を見てください。

24! = 620448401733239439360000
25! = 15511210043330985984000000

25 は 5 * 5 と素因数分解することができます。このとき、2 * 5 の組が 2 つできるので、末尾に 0 が 2 つ付くわけです。階乗を素因数分解したとき、因数 2 の個数は因数 5 の個数よりも多くなるので、2 と 5 は必ず組にすることができます。つまり、因数 5 の個数が末尾に付く 0 の 個数になるわけです。

階乗の場合、因数の個数を求めるのは簡単です。10000! の場合、10000 / 5 で 5 の倍数の個数 2000 を求めることができます。次に、25 (= 5 * 5) の倍数の個数を 10000 / 25 で求めます。さらに、125 (= 5 * 5 * 5) の倍数の個数を 10000 / 125 で求めます、これを 10000 > 5m が成立する m まで繰り返し、その総和が 5 の因子の個数になります。

10000 / 5    = 2000
10000 / 25   =  400
10000 / 125  =   80
10000 / 625  =   16
10000 / 3125 =    3.2 (小数点切捨て)
----------------------
        合計 = 2499

プログラムと実行結果を示します。

julia> function solver(n)
       a = 0
       m = big(5)
       while m <= n
       a += div(n, m)
       m *= 5
       end
       a
       end
solver (generic function with 1 method)

julia> for x in 1 : 10
       println(solver(big(10) ^ x))
       end
2
24
249
2499
24999
249998
2499999
24999999
249999998
2499999997

●解答3

Julia の場合、次のように簡単に求めることができます。

julia> (big(7) ^ 654321) % 10
7

この問題は筆算で簡単に求めることができます。7n (n > 0) の末尾の数字は次のように 7, 9, 3, 1, ... と巡回します。

7^1 = 7
7^2 = 49
7^3 = 343
7^4 = 2401
7^5 = 16807
7^6 = 117649
7^7 = 823543
7^8 = 5764801

ここで、74 の末尾は 1 になることに注目してください。末尾が 1 の数字を何回乗算しても、その結果の末尾は 1 になります。7654321 は (74)163580 * 7 なので、末尾の数字は 7 と求めることができます。


●解答4

米粒の合計値は 1 + 2 + 22 + 23 + ... + 280 になります。これを素直にプログラムすると次のようになります。

julia> sum([big(2) ^ x for x in 0 : 80])
2417851639229258349412351

とても大きな数になるので、普通の電卓では計算できません。Windows の電卓を使用するときは関数電卓に切り替えてください。もちろん、数学の公式を使うともっと簡単に求めることができます。

●等比数列の和

次のように、一定の比で並んだ数列を「等比数列」といいます。

a, ar, ar2, ..., arn-1, ...

a を「初項」、d を「公比」といいます。等比数列の一般項は次の式で表すことができます。

an = arn-1

初項から an までの和 Sn は次の式で求めることができます。

Sn = a(1 - rn) / (1 - r)

問題は初項 1 で公比 2 なので、米粒の合計は次のようになります。

(1 - 281) / (1 - 2) = 281 - 1
julia> big(2) ^ 81 - 1
2417851639229258349412351

この公式は次のように導出することができます。

Sn = a + ar + ar2 + ... + arn-1
両辺を r 倍すると
rSn =    ar + ar2 + ... + arn-1 + arn
これを引き算すると
Sn - rSn = a - arn => Sn = a(1 - rn) / (1 - r)

右辺を引き算すると ar から arn-1 の項がなくなって、a - arn だけになります。あとは、1 - r で割り算すればいいわけです。


●解答5

7 以上の素数で割り切れない正の整数は、素因子が 2, 3, 5 しかない自然数のことです。これを「ハミング数 (Hamming Numbers)」といいます。ハミング数は素因数分解したとき、2i * 3j * 5k (i, j, k >= 0) の形式になります。たとえば、100 以下のハミング数は次のようになります。

1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 
54, 60, 64, 72, 75, 80, 81, 90, 96, 100

それではプログラムを作りましょう。一番簡単な方法は、1 から n までの整数列を生成して、そこからハミング数を取り出していくことです。これを julia でプログラムすると次のようになります。

julia> function check(n)
       while n % 2 == 0; n = div(n, 2); end
       while n % 3 == 0; n = div(n, 3); end
       while n % 5 == 0; n = div(n, 5); end
       n == 1
       end
check (generic function with 1 method)

julia> hamming(n) = filter(check, 1 : n)
hamming (generic function with 1 method)

julia> for x in 2 : 8
       n = 10 ^ x
       m = length(hamming(n))
       println("$n, $m")
       end
100, 34
1000, 86
10000, 175
100000, 313
1000000, 507
10000000, 768
100000000, 1105

関数 check(n) は n がハミング数かチェックします。これは 2, 3, 5 だけで割り切れるか試しているだけです。プログラムはとても簡単ですが、引数 n の値が大きくなると時間がかかるようになります。n に比べてハミング数の個数は少ないようなので、式 2i * 3j * 5k (i, j, k >= 0) を使ってハミング数を生成したほうがよさそうです。引数 n に対して i, j, k の上限値は log2 n, log3 n, log5 n で求めることができます。たとえば、100000000 の場合は次のようになります。

i : 0 - 26
j : 0 - 16
k : 0 - 11

全体で 27 * 17 * 12 = 5508 個しかありません。この中から 100000000 以下の数を選べばいいわけです。プログラムは次のようになります。

リスト : ハミングの問題 (2)

function hamming(n)
    zs = []
    for x in [big(2) ^ i for i in 0 : floor(Int, log(2, n))]
        for y in [big(3) ^ j for j in 0 : floor(Int, log(3, n))]
            for z in [big(5) ^ k for k in 0 : floor(Int, log(5, n))]
                if x * y * z <= n
                    push!(zs, x * y * z)
                end
            end
        end
    end
    sort(zs)
end

2, 3, 5 のべき乗の集合を内包表記で生成し、その要素を掛け合わせて、条件を満たす数値を選択していくだけです。実行結果は次のようになりました。

julia> for x in 8 : 12
       n = big(10) ^ x
       m = length(hamming(n))
       println("$n, $m")
       end
100000000, 1105
1000000000, 1530
10000000000, 2053
100000000000, 2683
1000000000000, 3429

この方法だと短時間で答えを求めることができます。


関数で遊ぼう

●クロージャによる連結リストの実装

クロージャをサポートしているプログラミング言語では、効率を考慮しないでよければ、クロージャを使って「連結リスト」を実装しすることができます。連結リストの操作関数 cons, car, cdr には次の関係が成り立ちます。

x == car(cons(x, y)) => true
y == cdr(cons(x, y)) => true

ここで cons(x, y) で生成したオブジェクトがセルではない場合を考えてみましょう。もし、そのオブジェクトに car を適用すれば cons の第 1 引数 x を返し、cdr を適用すれば第 2 引数を返すことができれば、セルと同じことが実現できます。

そこで、cons はセルではなくクロージャを返すことにしましょう。クロージャは引数 x, y の値を保持することができます。そして、このクロージャは引数に関数を受け取ることにします。あとは、この関数に引数 x, y を渡して評価すれば car と cdr を実現することができます。Julia で cons, car, cdr をプログラムすると次のようになります。

julia> cons(x, y) = z -> z(x, y)
cons (generic function with 1 method)

julia> car(z) = z((x, y) -> x)
car (generic function with 1 method)

julia> cdr(z) = z((x, y) -> y)
cdr (generic function with 1 method)

関数 cons() はクロージャを返します。このクロージャは引数 z に関数を受け取り、その関数に x, y を渡して評価します。car() は引数 x にクロージャを渡して評価し、第 1 引数 a を返します。これで Lisp / Scheme の car と同じ動作になります。同様に、cdr() は引数 x にクロージャを渡して評価し、第 2 引数 b を返します。これで Lisp / Scheme の cdr と同じ動作になります。

それでは実際に試してみましょう。

julia> a = cons(1, 0)
(anonymous function)

julia> car(a)
1

julia> cdr(a)
0

julia> b = cons(2, a)
(anonymous function)

julia> car(b)
2

julia> cdr(b)
(anonymous function)

julia> car(cdr(b))
1

このように、クロージャを使って連結リストを作成することができます。ご参考までに簡単なリスト操作関数を下記リストに示します。

#
# list2.jl : クロージャによる連結リストの実装
#
#            Copyright (C) 2016 Makoto Hiroi
#

# 空リスト
const nil = :nil

# 基本関数
cons(x, y) = z -> z(x, y)
car(z) = z((x, y) -> x)
cdr(z) = z((x, y) -> y)

# 述語
null(x) = x == nil
consp(x) = typeof(x) == Function
listp(x) = null(x) || consp(x)
atom(x) = !consp(x)

# リストの表示
function printlist(xs)
    print("(")
    while consp(xs)
        if listp(car(xs))
            printlist(car(xs))
        else
            print(car(xs))
        end
        if consp(cdr(xs)); print(" "); end
        xs = cdr(xs)
    end
    if !null(xs); print(" . $xs"); end
    print(")")
end

# リストの生成
function list(args...)
    xs = nil
    for i in length(args) : -1 : 1
        xs = cons(args[i], xs)
    end
    xs
end

# リストの連結
function append(xs, ys)
    if null(xs)
        ys
    else
        cons(car(xs), append(cdr(xs), ys))
    end
end

# マッピング
function mapcar(f, xs)
    if null(xs)
        nil
    else
        cons(f(car(xs)), mapcar(f, cdr(xs)))
    end
end

# フィルター
function remove(f, xs)
    if null(xs)
        nil
    elseif f(car(xs))
        remove(f, cdr(xs))
    else
        cons(car(xs), remove(f, cdr(xs)))
    end
end

# 畳み込み
function fold_left(f, a, xs)
    while consp(xs)
        a = f(a, car(xs))
        xs = cdr(xs)
    end
    a
end

function fold_right(f, a, xs)
    if null(xs)
        a
    else
        f(car(xs), fold_right(f, a, cdr(xs)))
    end
end
julia> include("list2.jl")
fold_right (generic function with 1 method)

julia> a = list(1,2,3,4,5)
(anonymous function)

julia> printlist(a)
(1 2 3 4 5)
julia> b = list(6.7.8.9.10)
ERROR: syntax: invalid numeric constant "6.7."

julia> b = list(6,7,8,9,10)
(anonymous function)

julia> printlist(append(a, b))
(1 2 3 4 5 6 7 8 9 10)
julia> printlist(mapcar(x -> x * x, a))
(1 4 9 16 25)
julia> printlist(remove(x -> x % 2 == 0, a))
(1 3 5)
julia> printlist(fold_left(+, 0, a))
( . 15)
julia> fold_left(+, 0, a)
15

julia> fold_right(+, 0, a)
15

julia> fold_right(cons, nil, a)
(anonymous function)

julia> printlist(fold_right(cons, nil, a))
(1 2 3 4 5)
julia> printlist(fold_left((x, y) -> cons(y, x), nil, a))
(5 4 3 2 1)

●無名関数とカリー化関数

Julia の場合、関数をデータ型の一つとして扱うことができるので、関数の返り値として関数を返すことができます。この「関数を返す関数」を使うと、関数の引数が一つでも複数の引数を処理することができます。このような関数を「カリー化関数 (curried function)」といいます。

たとえば、(x, y) -> x + y をカリー化関数にする場合、引数 x を受け取ると「引数 y を受け取って x + y を計算する関数」を返し、その関数に引数 y を渡せば x + y を計算することができます。Julia では、無名関数を使って次のように定義することができます。

julia> foo = x -> y -> x + y
(anonymous function)

無名関数を定義する -> は右結合なので、x -> y -> x + y は x -> (y -> x + y) と同じになります。これで引数をひとつ受け取ったら、関数 y -> x + y を返すことができます。もちろん、引数を 2 つ与えれば、それらを加算した結果を返します。つまり、最初の引数を受け取って関数を生成し、その関数を 2 番目の引数に適用する、という動作になります。

それでは実際に試してみましょう。

julia> foo1 = foo(1)
(anonymous function)

julia> foo1(2)
3

julia> foo(10)(20)
30

引数を一つだけ渡すと「引数 y を受け取って x + y を計算する関数」を返します。返り値を変数 foo1 にセットして、foo1(2) を呼び出せば 1 + 2 = 3 を計算することができます。もちろん、引数を 2 つ渡すと、それらを足し算した値を返します。foo(10) の返り値は関数なので、2 番目の引数を渡すときは foo(10)(20) のようにカッコを続けることに注意してください。

また、局所的な関数を定義して、その関数を返すことでもカリー化関数を実現することができます。次の例を見てください。

julia> function mapcar(f)
       function _map(xs)
       map(f, xs)
       end
       _map
       end
mapcar (generic function with 1 method)

関数 mapcar() はマッピングを 1 引数の関数に直したものです。mapcar() は関数 f を受け取り、その f を呼び出してリストを操作する関数を返します。これでもマッピングの動作ができるのです。簡単な例を示しましょう。

julia> foo2 = mapcar(x -> x * 2)
_map (generic function with 1 method)

julia> foo2([1,2,3,4,5])
5-element Array{Int32,1}:
  2
  4
  6
  8
 10

julia> mapcar(x -> x * x)([1,2,3,4,5])
5-element Array{Int32,1}:
  1
  4
  9
 16
 25

最初の例は mapcar() で生成した関数を変数 foo2 にセットし、それから foo2 を関数呼び出しします。次の例は、mapcar() の返り値を直接関数呼び出ししています。カッコが多くなりますが、2 引数の map() と同じように呼び出すことができます。これでもリストの要素を 2 乗することができます。

Julia では function 文でカリー化関数を定義することはできませんが、関数型言語 (たとえば Haskell や ML (SML/NJ, Ocaml) など) では、簡単にカリー化関数を定義できるようになっています。これらのプログラミング言語では、高階関数はカリー化関数として定義されていて、関数を合成して新しい関数を作ることも簡単にできます。興味のある方は関数型言語にも挑戦してみてください。

●チャーチ数

「ラムダ計算 (lambda calculus)」は、文字λを使って関数を表す「λ記法」という表記法を用いた抽象的な計算モデルで、1930 年代に A. Church 氏によって考案されました。ラムダ計算は Lisp, Scheme, ML, Haskell など多くの関数型言語の基礎理論として、大きな役割を果たしています。ラムダ計算とか計算モデルというと難しい話のように思われるかもしれません。Lisp / Scheme では無名関数のことを「ラムダ式」といいますが、実をいうとこの無名関数の考え方がラムダ計算の基本なのです。

純粋なラムダ計算の定義はとても単純です。ラムダ計算で扱う式は、次に示す 3 通りしかありません。

関数抽象は関数定義、関数適用は関数呼び出し、変数は関数の仮引数のことと考えてください。つまり、純粋なラムダ計算には関数しかないのです。したがって、ラムダ計算では数を表すのにも関数を使います。これを「チャーチ数 (Church numerals)」と呼びます。

ラムダ計算と無名関数はまったく同じではありせんが、今回は難しいことを考えずに Julia の無名関数を使って「チャーチ数」を試してみましょう。

●チャーチ数の基本

チャーチ数は関数 f と x を受け取り、x に f を適用した回数で数 (自然数) を表します。たとえば、自然数 n は f(f(f(...(f(x)) ...))) のように f を n 回呼び出すことで表します。簡単な例を示しましょう。

julia> zero = f -> x -> x
(anonymous function)

julia> one = f -> x -> f(x)
(anonymous function)

julia> two = f -> x -> f(f(x))
(anonymous function)

julia> three = f -> x -> f(f(f(x)))
(anonymous function)

変数 zero に格納された関数は、引数 f を受け取ったら関数 x -> x を返す関数に、引数 f と x を 2 つ受け取ったら x をそのまま返す関数になります。このとき、x に f を適用していないことに注意してください。これで 0 を表すことができます。

同様に、one は f を 1 回適用しているので 1 を、two は 2 回適用しているので 2 を、three は 3 回適用しているので 3 を表すことができます。そうはいっても、このままではよくわからないので、実際に引数として関数 n -> n + 1 と 0 を渡して実行してみましょう。実行結果は次のようになります。

julia> zero(n -> n + 1)(0)
0

julia> one(n -> n + 1)(0)
1

julia> two(n -> n + 1)(0)
2

julia> three(n -> n + 1)(0)
3

n -> n + 1 は引数 n に 1 を加算するので、もう一つの引数に 0 を渡せば n -> n + 1 を適用した回数、つまりチャーチ数を Julia の数に変換することができます。

●足し算

次は数 n に 1 を加える succ(f)(n)(x) を定義してみましょう。このとき、n はチャーチ数であることに注意してください。プログラムは次のようになります。

julia> succ = f -> n -> x -> f(n(f)(x))
(anonymous function)

n(f)(x) は数 n を表しているので、それに関数 f を再度適用すれば、n に 1 を加えることができます。簡単な実行例を示します。

julia> succ(n -> n + 1)(zero)(0)
1

julia> succ(n -> n + 1)(one)(0)
2

julia> succ(n -> n + 1)(two)(0)
3

julia> succ(n -> n + 1)(three)(0)
4

このように、succ でチャーチ数 zero, one, two, three に 1 を加算することができます。

次は 2 つのチャーチ数 m, n を足し算する plus(f)(m)(n)(x) を作りましょう。この場合、n(f)(x) で n を表すチャーチ数になるので、この結果に m(f) を適用すれば m + n を実現することができます。プログラムは次のようになります。

julia> plus = f -> m -> n -> x -> m(f)(n(f)(x))
(anonymous function)

julia> plus(n -> n + 1)(one)(two)(0)
3

julia> plus(n -> n + 1)(two)(two)(0)
4

julia> plus(n -> n + 1)(two)(three)(0)
5

正常に動作していますね。

●掛け算

次はチャーチ数 m と n を掛け算する mult(f)(m)(n)(x) を定義してみましょう。m * n は n を m 回足し算すればいいので、関数 n(f) を m に渡して m(n(f))(x) とするだけです。プログラムと実行結果を示します。

julia> mult = f -> m -> n -> x -> m(n(f))(x)
(anonymous function)

julia> mult(n -> n + 1)(two)(three)(0)
6

julia> mult(n -> n + 1)(three)(three)(0)
9

julia> mult(n -> n + 1)(three)(zero)(0)
0

julia> mult(n -> n + 1)(zero)(three)(0)
0

このように、チャーチ数の足し算と掛け算は簡単なのですが、引き算はとても難しくなります。本稿の範囲を超える (M.Hiroi が理解できない) ので、チャーチ数はここまでにしておきましょう。興味のある方は調べてみてください。

●参考文献, URL

  1. Ravi Sethi (著), 神林靖 (訳), 『プログラミング言語の概念と構造』,アジソンウェスレイ, 1995
  2. ラムダ計算入門 (PDF), (住井英二郎さん)
  3. ラムダ計算 - Wikipedia

Copyright (C) 2016 Makoto Hiroi
All rights reserved.

[ Home | Light | Julia ]