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

Julia Language Programming

簡単なプログラム (5)

[ Home | Light | Julia ]

●カッコ列

カッコ列は ( と ) からなる列のことで、バランスが取れているカッコ列は、右カッコで閉じることができる、つまり右カッコに対応する左カッコがある状態のことをいいます。たとえば n = 1 の場合、( ) はバランスの取れたカッコ列ですが、) ( はバランスが取れていません。今回はカッコ列を生成する高階関数 kakko() を作ります。

リスト : カッコ列

function kakko(f, n)
    function iter(x, y, a)
        if x == y == n
            f(a)
        else
            if x < n
                iter(x + 1, y, a * "(")
            end
            if y < x
                iter(x, y + 1, a * ")")
            end
        end
    end
    iter(0, 0, "")
end

kakko(println, 3)
kakko(println, 4)

カッコ列の生成は簡単です。局所関数 iter() の引数 x が左カッコの個数、引数 y が右カッコの個数を表します。引数 a は累積変数でカッコ列を表す文字列です。

バランスの取れたカッコ列の場合、x, y, n には y <= x <= n の関係が成り立ちます。x == y == n の場合、カッコ列がひとつ完成しました。引数の関数 f を呼び出します。そうでなければ、iter() を再帰呼び出しします。x < n であれば左カッコを追加し、y < x であれば右カッコを追加します。これでカッコ列を生成することができます。

((()))
(()())
(())()
()(())
()()()
(((())))
((()()))
((())())
((()))()
(()(()))
(()()())
(()())()
(())(())
(())()()
()((()))
()(()())
()(())()
()()(())
()()()()

●カッコ列の総数

カタラン数 - Wikipedia によると、カッコ列の総数は「カタラン数 (Catalan number)」になるとのことです。カタラン数は次に示す公式で求めることができます。

         (2n)!
Cn = ----------
       (n+1)!n!

これをそのままプログラムしてもいいのですが、それではちょっと面白くないので別な方法でプログラムを作ってみましょう。カタラン数は次に示す経路図において、A から B までの最短距離の道順を求めるとき、対角線を超えないものの総数に一致します。

                    B                      B  
  ┌─┬─┬─┬─┐      ┌─┬─┬─0─14    
  │  │  │  │  │      │  │  │  │  │    
  ├─┼─┼─┼─┤      ├─┼─0─5─14    
  │  │  │  │  │      │  │  │  │  │    
  ├─┼─┼─┼─┤      ├─0─2─5─9    
  │  │  │  │  │      │  │  │  │  │    
  ├─┼─┼─┼─┤      0─1─2─3─4    
  │  │  │  │  │      │  │  │  │  │    
  └─┴─┴─┴─┘      1─1─1─1─1    
A                      A                      

            図 : 道順の総数の求め方

A からある地点にいたる最短距離の道順の総数は、左隣と真下の地点の値を足したものになります。一番下の地点は 1 で、対角線を越えた地点は 0 になります。あとは下から順番に足し算していけば、A から B までの道順の総数を求めることができます。上図の場合はカラタン数 C4 に相当し、その値は 14 となります。

プログラムは配列を使うと簡単です。次の図を見てください。

0 : [1, 1, 1, 1, 1]

1 : [1, 1, 1, 1, 1,]

2 : [1, 1, 1+1=2, 2+1=3, 3+1=4]
 => [1, 1, 2, 3, 4]

3 : [1, 1, 2, 3+2=5, 5+4=9]
 => [1, 1, 2, 5, 9]

4 : [1, 1, 2, 5, 5+9=14]
 => [1, 1, 2, 5, 14]

上図は Cn (n = 4) を求める場合です。大きさが n + 1, 要素の値が 1 のベクタを用意します。n = 0, 1 の場合は n 番目の要素をそのまま返します。n が 2 よりも大きい場合、変数 i を 2 に初期化して、i - 1 番目以降の要素の累積和を求めます。

たとえば i = 2 の場合、2 番目の要素は 1 番目の要素と自分自身を加算した値 2 になります。3 番目の要素は 2 番目の要素と自分自身を足した値 3 になり、4 番目の要素は 3 + 1 = 4 になります。次に i を +1 して同じことを繰り返します。3 番目の要素は 2 + 3 = 5 になり、4 番目の要素は 5 + 4 = 9 になります。i = 4 のとき、4 番目の要素は 5 + 9 = 14 となり、C4 の値を求めることができました。

リスト : カッコ列の総数

function kakko_num(n)
    table = Array(BigInt, n)
    fill!(table, big(1))
    for i in  2 : n
        for j in i : n
            table[j] += table[j - 1]
        end
    end
    table[n]
end

for i in 1 : 10
    println(kakko_num(i))
end

println(kakko_num(50))
println(kakko_num(100))
1
2
5
14
42
132
429
1430
4862
16796
1978261657756160653623774456
896519947090131496687170070074100632420837521538745909320

●カッコ列と二分木

バランスの取れたカッコ列と二分木は 1 対 1 に対応します。二分木を行きがけ順で巡回するとき、途中の節では左カッコ ( を出力して左右の枝をたどり、葉に到達したら右カッコ ) を出力すると、カッコ列を生成することができます。

リスト : カッコ列と二分木

# 二分木
abstract Tree

# 節
immutable Node <: Tree
    left::Tree
    right::Tree
end

# 葉
immutable Leaf <: Tree
end

# 二分木をカッコ列に変換
function tree_to_kakko(tree::Tree)
    function iter(tree::Tree)
        if typeof(tree) == Leaf
            ")"
        else
            "(" * iter(tree.left) * iter(tree.right)
        end
    end
    s = iter(tree)
    s[1:end-1]
end

# カッコ列を二分木に変換
function kakko_to_tree(s)
    function iter(i)
        if length(s) < i
            (Leaf(), i)
        elseif s[i] == ')'
            (Leaf(), i + 1)
        elseif s[i] == '('
            l, j = iter(i + 1)
            r, k = iter(j)
            (Node(l, r), k)
        else
            error("illegal char")
        end
    end
    iter(1)[1]
end

function test(s)
    println(s)
    t = kakko_to_tree(s)
    println(t)
    u = tree_to_kakko(t)
    println(u)
end

kakko(test, 3)

関数 tree_to_kakko() の実際の処理は局所関数 iter() で行います。引数が葉 (Leaf) の場合は右カッコを返します。引数が節 (Node) の場合は、iter() を再帰呼び出しして左部分木 left をたどり、それから右部分木 right をたどります。その結果と左カッコを演算子 * で連結すればいいわけです。ただし、最後に余分な右カッコが付いてくるので、最後の要素を削除します。

二分木の場合、葉 (Leaf) の個数を n とすると、節 (Node) の個数は n - 1 になります。カッコ列と違って Leaf の個数が一つ多くなることに注意してください。

関数 kakko_to_tree() も局所関数 iter() で変換処理を行います。文字列 s の i 番目の要素が右カッコの場合は (Leaf(), i + 1) を返します。左カッコの場合は、iter() を再帰呼び出しして左部分木 l を生成し、それから右部分木 r を生成します。あとは (Node(l, r), k) を返すだけです。ただし、葉に対応する右カッコがひとつ少ないので、引数 i が文字列の終端に到達した場合は Leaf() を返すようにします。

((()))
Node(Node(Node(Leaf(),Leaf()),Leaf()),Leaf())
((()))
(()())
Node(Node(Leaf(),Node(Leaf(),Leaf())),Leaf())
(()())
(())()
Node(Node(Leaf(),Leaf()),Node(Leaf(),Leaf()))
(())()
()(())
Node(Leaf(),Node(Node(Leaf(),Leaf()),Leaf()))
()(())
()()()
Node(Leaf(),Node(Leaf(),Node(Leaf(),Leaf())))
()()()

●累乗

累乗は下図のように x の n 乗という x を n 回掛ける計算です。Julia には累乗を計算する演算子 ^ がありますが、繰り返しでも再帰定義でも簡単にプログラムすることができます。

x ^ 3 = x * x * x
x ^ 4 = x * x * x * x
x ^ 5 = x * x * x * x * x

  図 : 累乗の計算

累乗を単純な繰り返しで実装すると、次のようになります。

リスト : 累乗の計算 (1)

function power0(x::BigInt, n::Int)
    a = big(1)
    for _ in 1 : n; a *= x; end
    a
end

この場合、n 回の乗算が必要になります。ところが、式を変形するともっと少ない回数で求めることができます。

x ^ 4  = (x ^ 2) ^ 2 -> 2 回
x ^ 8  = (x ^ 4) ^ 2 -> 3 回
x ^ 16 = (x ^ 8) ^ 2 -> 4 回

  一般化すると

x ^ n = (x ^ (n / 2)) ^ 2; (n は偶数)

x ^ n = ((x ^ (n / 2)) ^ 2) * x; (n は奇数)

  図 : 累乗を求める式の変形

階乗計算では n を n - 1 の計算に置き換えていきますが、累乗の場合は n を n / 2 に置き換えていくことができます。n が半分になっていくので、ひとつずつ n を減らすよりも減少の度合いは大きくなります。その分だけ計算回数が少なくなるわけです。

それでは、この考え方をプログラムしてみましょう。x ^ (n / 2) を計算する部分は、再帰を使えば簡単です。

リスト : 累乗の計算 (2)

function power1(x::BigInt, n::Int)
    if n == 0
        big(1)
    elseif n % 2 == 1
        x * power1(x, n - 1)
    else
        v = power1(x, div(n, 2))
        v * v
    end
end

関数 power1() は一般化した累乗の定義をそのままプログラムしただけです。最初の if 文が再帰呼び出しの停止条件です。次に、n が奇数ならば x と power1(x, n - 1) を掛け算します。n が偶数ならば power1() で x ^ (n/2) を求めて、 それを二乗した値を返します。

ところで、関数 power1() も繰り返しに変換することができます。整数 n は下図のように 2 ^ i の和の形で表せることを利用します。

11 = 1 + 2 + 8     (2^0 + 2^1 + 2^3)
15 = 1 + 2 + 4 + 8 (2^0 + 2^1 + 2^2 + 2^3)

        図 : 整数の 2 進数表現

これは整数 n の下位ビットから順番にビットをチェックしていけば簡単に求めることができます。これを利用すると x ^ n は下図のように求めることができます。

x ^ 11 = x * x ^ 2 * x ^ 8
x ^ 15 = x * x ^ 2 * x ^ 4 * x ^ 8

x ^ 2 = x * x
x ^ 4 = x ^ 2 * x ^ 2
x ^ 8 = x ^ 4 * x ^ 4
  .....

        図 : 累乗の計算

x ^ 2, x ^ 4, x ^ 8, ... は値を次々に 2 乗していけば簡単に求めることができます。これをプログラムすると次のようになります。

リスト : 累乗の計算 (3)

function power2(x::BigInt, n::Int)
    v = big(1)
    while n > 0
        if n & 1 != 0; v *= x; end
        n >>= 1
        x *= x
    end
    v
end

最初に局所変数 v を big(1) に初期化します。そして、n が 0 よりも大きいあいだ、つまり 1 のビットがあるあいだ処理を繰り返します。最下位ビットが 1 の場合は v に x を掛け算します。そして、n を 1 ビット右へシフトします。x の値は x *= x で自乗されていくので、x, x ^ 2, x ^ 4, x ^ 8, x ^ 16 と増えていきます。繰り返しを終了したら v を返します。

それでは、実際に実行速度を比較してみましょう。

リスト : 簡単なテスト

function test(f)
    for _ in 1 : 100
        f(big(2), 100000)
    end
end

@time test(power0)
@time test(power1)
@time test(power2)
@time test(^)
 25.185582 seconds (40.01 M allocations: 58.487 GB, 12.18% gc time)
  0.057852 seconds (10.99 k allocations: 2.554 MB)
  0.142687 seconds (11.28 k allocations: 4.888 MB, 1.22% gc time)
  0.006189 seconds (1.15 k allocations: 1.213 MB)

組み込み演算子 ^ が一番速くなりました。power0() と比べると、ower1() と power2() はとても速いですね。そして、power2() よりも power1() の方が速くなったのには驚きました。power1() のような再帰定義でも十分実用的といえるでしょう。

●平方根

実数 a の平方根 √a の値を求める場合、方程式 x2 - a = 0 を Newton (ニュートン) 法で解くことが多いと思います。方程式を f(x), その導関数を f'(x) とすると、ニュートン法は次の漸化式の値が収束するまで繰り返す方法です。

xn+1 = xn - f(xn) / f'(xn)

平方根を求める場合、導関数は f'(x) = 2x になるので、漸化式は次のようになります。

xn+1 = (xn + a / xn) / 2

参考文献 1 によると、√a より大きめの初期値から出発し、置き換え x <- (x + a / x) / 2 を減少が止まるまで繰り返すことで √a の正確な値を求めることができるそうです。

Julia でプログラムすると、次のようになります。

リスト : 平方根を求める

function sqrt1(x::Float64)
    function init(x::Float64, s = 1.0)
        while s < x
            s *= 2.0
            x /= 2.0
        end
        s
    end
    if x < 0; error("sqrt1: domain error"); end
    p = x > 1.0 ? init(x) : 1.0
    while true
        q = (p + x / p) / 2
        if q >= p; break; end
        p = q
    end
    p
end

println(sqrt1(2.0))
println(sqrt(2))
println(sqrt1(123456.0))
println(sqrt(123456))
1.414213562373095
1.4142135623730951
351.363060095964
351.363060095964

局所関数 init は √x よりも大きめの初期値を求めます。たとえば、√123456 を求める場合、初期値の計算は次のようになります。

   s         x
-------------------
  1.0  123456.0
  2.0   61728.0
  4.0   30864.0
  8.0   15432.0
 16.0    7716.0
 32.0    3858.0
 64.0    1929.0
128.0     964.5
256.0     482.25
512.0     241.125

√123456 = 351.363060095964 

s を 2 倍、x を 1 / 2 していき、s >= x となったときの s が初期値 (512) となります。4, 16, 64, 256, ... 22n の平方根はこれだけで求めることができます。

あとは漸化式を計算して変数 q にセットし、q がひとつ前の値 p 以上になったら p を返すだけです。√123456 を求めたときの p と q の値を示します。

   p                  q
--------------------------------------
512.0              376.5625
376.5625           352.20622925311204
352.20622925311204 351.3640693544162
351.3640693544162  351.3630600974135
351.3630600974135  351.363060095964
351.363060095964   351.363060095964

√123456 = 351.363060095964 

6 回の繰り返しで √123456 を求めることができます。

●めのこ平方

平方根の整数部もニュートン法を使って求めることができますが、次の公式を使って平方根の整数部分を求めることもできます。

(1) 1 + 3 + 5 + ... + (2n - 1) = n2
(2) 1 + 3 + 5 + ... + (2n - 1) = n2 < m < 1 + 3 + ... (2n - 1) + (2n + 1) = (n + 1)2

式 (1) は、奇数 1 から 2n - 1 の総和は n2 になることを表しています。式 (2) のように、整数 m の値が n2 より大きくて (n + 1)2 より小さいのであれば、m の平方根の整数部分は n であることがわかります。これは m から奇数 1, 3, 5 ... (2n - 1), (2n + 1) を順番に引き算していき、引き算できなくなった時点の (2n + 1) / 2 = n が m の平方根になります。参考文献 2 によると、この方法を「めのこ平方」と呼ぶそうです。

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

リスト : めのこ平方

# めのこ平方
function isqrt(n, m = 1)
    while n >= m
        n -= m
        m += 2
    end
    div(m, 2)
end

println(isqrt(4))
println(isqrt(16))
println(isqrt(64))
println(isqrt(80))
println(isqrt(81))
println(isqrt(82))
println(isqrt(100))
2
4
8
8
9
9
10

この方法はとても簡単ですが、数が大きくなると時間がかかるようになります。そこで、整数を 2 桁ずつ分けて計算していくことにします。次の図を見てください。

整数 6789 を 67 と 89 に分ける

1 + 3 + ... + 15 = 82 < 67

両辺を 100 倍すると 802 < 6700 < 6789

802 = 1 + 3 + ... + 159 (= 2 * 80 - 1)

161 + 163 < (6789 - 6400 = 389) < 161 + 163 + 165

整数 6789 を 67 と 89 に分けます。最初に 67 の平方根を求めます。この場合は 8 になり、82 < 67 を満たします。次に、この式を 100 倍します。すると、802 < 6700 になり、6700 に 89 を足した 6789 も 802 より大きくなります。802 は 1 から 159 までの奇数の総和であることはすぐにわかるので、6789 - 6400 = 389 から奇数 161, 163, ... を順番に引き算していけば 6789 の平方根を求めることができます。

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

リスト : めのこ平方 (改良版)

function isqrt1(n)
    if n < 100
        isqrt(n)
    else
        m = 10 * isqrt1(div(n, 100))
        isqrt(n - m * m, 2 * m + 1)
    end
end

println(isqrt1(6789))
println(isqrt1(123456789))
println(isqrt1(1234567890))
82
11111
35136

isqrt1() は n の平方根の整数部分を求めます。n が 100 未満の場合は isqrt() で平方根を求めます。これが再帰呼び出しの停止条件になります。n が 100 以上の場合は、n の下位 2 桁を取り除いた値 div(n, 100) の平方根を isqrt1() で求め、その値を 10 倍して変数 m にセットします。そして、isqrt() で n - m * m から奇数 2 * m + 1, 2 * m + 3 ... を順番に引き算していって n の平方根を求めます。

興味のある方はいろいろ試してみてください。

●参考文献

  1. 奥村晴彦,『C言語による最新アルゴリズム事典』, 技術評論社, 1991
  2. 仙波一郎のページ, 『平方根計算法 (PDF)』

●べき集合

今回は配列 xs のべき集合を求める高階関数 power_set() を作ります。たとえば配列 [1, 2, 3] のべき集合は [], [1], [2], [3], [1, 2], [1, 3], [2, 3], [1, 2, 3] になります。

リスト : べき集合

function power_set(f, xs)
    function iter(i, a)
        if i > length(xs)
            f(a)
        else
            iter(i + 1, a)
            push!(a, xs[i])
            iter(i + 1, a)
            pop!(a)
        end
    end
    a::typeof(xs) = []
    iter(1, a)
end

power_set(println, ["foo", "bar", "baz"])
power_set(println, [1,2,3,4])

power_set() は簡単です。実際の処理は局所関数 iter() で行います。xs の i 番目の要素を選択する場合は、その要素を引数 a に追加して iter() を再帰呼び出しします。選択しない場合は、引数 a に要素を追加せずに iter() を再帰呼び出しするだけです。これでべき集合の要素をすべて求めることができます。

ASCIIString[]
ASCIIString["baz"]
ASCIIString["bar"]
ASCIIString["bar","baz"]
ASCIIString["foo"]
ASCIIString["foo","baz"]
ASCIIString["foo","bar"]
ASCIIString["foo","bar","baz"]
Int32[]
[4]
[3]
[3,4]
[2]
[2,4]
[2,3]
[2,3,4]
[1]
[1,4]
[1,3]
[1,3,4]
[1,2]
[1,2,4]
[1,2,3]
[1,2,3,4]

●マスターマインド

「マスターマインド」は 0 から 9 までの重複しない 4 つの数字からなる隠しコードを当てるゲームです。数字は合っているが位置が間違っている個数を cows で表し、数字も位置も合っている個数を bulls で表します。bulls が 4 になると正解です。

   [6, 2, 8, 1] : 正解
---------------------------------
1: [0, 1, 2, 3] : cows 2 : bulls 0
2: [1, 0, 4, 5] : cows 1 : bulls 0
3: [2, 3, 5, 6] : cows 2 : bulls 0
4: [3, 2, 7, 4] : cows 0 : bulls 1
5: [3, 6, 0, 8] : cows 2 : bulls 0
6: [6, 2, 8, 1] : cows 0 : bulls 4

  図 : マスターマインドの動作例

今回はマスターマインドを解くプログラムを作ることにします。

●推測アルゴリズム

このゲームは 10 個の数字の中から 4 個選ぶわけですから、全体では 10 * 9 * 8 * 7 = 5040 通りのコードしかありません。この中から正解を見つける方法ですが、質問したコードとその結果を覚えておいて、それと矛盾しないコードを作るようにします。具体的には、4 つの数字の順列を生成し、それが今まで質問したコードと矛盾しないことを確かめます。これは生成検定法と同じですね。

矛盾しているかチェックする方法も簡単で、以前に質問したコードと比較して、bulls と cows が等しいときは矛盾していません。たとえば、次の例を考えてみてください。

[6, 2, 8, 1] が正解の場合

[0, 1, 2, 3] => bulls = 0, cows = 2

           [0, 1, 2, 3]  と比較する
     --------------------------------------------------------
           [0, X, X, X]  0 から始まるコードは bulls = 1
                         になるので矛盾する。
           ・・・・

           [1, 0, 3, 4]  cows = 3, bulls = 0 になるので矛盾する

           ・・・・

           [1, 0, 4, 5]  cows = 2, bulls = 0 で矛盾しない。
     --------------------------------------------------------

[1, 0, 4, 5] => bulls = 0, cows = 1

次は、[0, 1, 2, 3] と [1, 0, 4, 5] に矛盾しない数字を選ぶ

        図 : マスターマインドの推測アルゴリズム

[0, 1, 2, 3] で bulls が 0 ですから、その位置にその数字は当てはまりません。したがって、[0, X, X, X] というコードは [0, 1, 2, 3] と比較すると bulls が 1 となるので、矛盾していることがわかります。

次に [1, 0, 3, 4] というコードを考えてみます。[0, 1, 2, 3] の結果は cows が 2 ですから、その中で合っている数字は 2 つしかないわけです。ところが、[1, 0, 3, 4] と [0, 1, 2, 3] と比較すると cows が 3 になります。当たっている数字が 2 つしかないのに、同じ数字を 3 つ使うのでは矛盾していることになりますね。

次に [1, 0, 4, 5] というコードと比較すると、bulls が 0 で cows が 2 となります。これは矛盾していないので、このコードを質問することにします。その結果が bulls = 0, cows = 1 となり、今度は [0, 1, 2, 3] と [1, 0, 4, 5] に矛盾しないコードを選択するのです。

●プログラムの作成

それでは、プログラムを作っていきましょう。まず、質問したコードとその結果を格納するデータ型を定義します。

リスト : データ型の定義

# 定数
const CSIZE = 4

# 質問したコードとその結果
type Query
    bulls::Int
    cows::Int
    code::Array{Int, 1}
end

型名は Query としました。bulls, cows と質問したコード code を格納します。これを大域変数 query の配列に格納します。

次は bulls を数える関数 count_bulls() を作ります。

リスト : bulls を数える

function count_bulls(xs, ys)
    c = 0
    for i in 1 : CSIZE
        if xs[i] == ys[i]; c += 1; end
    end
    c
end

count_bulls() は簡単ですね。配列 xs, ys の要素を順番に比較して、等しい場合は変数 c の値を +1 します。

次は cows を数える処理を作ります。いきなり cows を数えようとすると難しいのですが、2 つのリストに共通の数字を数えることは簡単にできます。この方法では、bulls の個数を含んだ数を求めることになりますが、そこから bulls を引けば cows を求めることができます。関数名は count_same_number() としましょう。プログラムは次のようになります。

リスト : 同じ数字の個数を数える

function count_same_number(xs, ys)
    c = 0
    for x in xs
        for y in ys
            if x == y
                c += 1
                break
            end
        end
    end
    c
end

for ループで xs の要素を順番に取り出して変数 x にセットします。次の for ループで x が ys に含まれているかチェックして、そうであれば c の値を +1 します。

次は、今まで質問したコードと矛盾していないか調べる関数 check を作ります。

リスト : 今まで質問したコードと矛盾していないか

function check(answer, xs)
    global query
    for q in query
        b = count_bulls(q.code, xs)
        c = count_same_number(q.code, xs) - b
        if b != q.bulls || c != q.cows
            return
        end
    end
    b = count_bulls(answer, xs)
    c = count_same_number(answer, xs) - b
    q = Query(b, c, xs)
    push!(query, q)
    n = length(query)
    println("$n: $xs, bulls = $b, cows = $c")
    if b == 4
        throw("Good Job!")
    end
end

引数 answer は正解のコード、xs は生成したコードです。最初に、大域変数 query に格納されたデータをチェックしていきます。count_bulls() と count_same_number() を使って bulls (変数 b) と cows (変数 c) を求めて、質問したときの q.bulls と q.cows に矛盾しないかチェックします。矛盾している場合は return で終了します。

それから、正解のコード answer と xs を比較して bulls と cows を求め、それらを構造体 Query にまとめて query に追加します。あとは関数 permutations() で順列を生成するだけです。詳細は プログラムリスト をお読みください。

●何回で当たるか

これでプログラムは完成です。それでは実行例を示しましょう。

julia> include("mastermind.jl")
solver (generic function with 1 method)

julia> solver([9,8,7,6])
1: [0,1,2,3], bulls = 0, cows = 0
2: [4,5,6,7], bulls = 0, cows = 2
3: [5,4,8,9], bulls = 0, cows = 2
4: [6,7,9,8], bulls = 0, cows = 4
5: [8,9,7,6], bulls = 2, cows = 2
6: [9,8,7,6], bulls = 4, cows = 0
Good Job!

julia> solver([9,4,3,1])
1: [0,1,2,3], bulls = 0, cows = 2
2: [1,0,4,5], bulls = 0, cows = 2
3: [2,3,5,4], bulls = 0, cows = 2
4: [3,4,0,6], bulls = 1, cows = 1
5: [3,5,6,1], bulls = 1, cows = 1
6: [6,5,0,2], bulls = 0, cows = 0
7: [7,4,3,1], bulls = 3, cows = 0
8: [8,4,3,1], bulls = 3, cows = 0
9: [9,4,3,1], bulls = 4, cows = 0
Good Job!

肝心の質問回数ですが、5, 6 回で当たる場合が多いようです。実際に、5040 個のコードをすべて試してみたところ、平均は 5.56 回になりました。これは 参考文献 1 の結果と同じです。質問回数の最大値は 9 回で、そのときのコードは [9 4 3 1], [9 2 4 1], [5 2 9 3], [9 2 0 4], [9 2 1 4] でした。

なお、参考文献 1 には平均質問回数がこれよりも少なくなる方法が紹介されています。単純な数当てゲームだと思っていましたが、その奥はけっこう深いようです。興味のある方はいろいろ試してみてください。

●参考文献

  1. 田中哲郎 「数当てゲーム (MOO, マスターマインド) 」, 松原仁、竹内郁雄 編 『bit 別冊 ゲームプログラミング』 pp150 - 157, 共立出版, 1997

●プログラムリスト

#
# mastermind.jl : マスターマインドの解法
#
#                 Copyright (C) 2016 Makoto Hiroi
#

# 定数
const CSIZE = 4

# 質問したコードとその結果
type Query
    bulls::Int
    cows::Int
    code::Array{Int, 1}
end

# 0 - 9 から 4 個の数字を選ぶ順列を生成
function permutations(f, xs, n = 1)
    if n > CSIZE
        f(xs[1:CSIZE])
    else
        tmp = xs[n]
        for i in n : length(xs)
            xs[n] = xs[i]
            xs[i] = tmp
            permutations(f, xs, n + 1)
            xs[i] = xs[n]
            xs[n] = tmp
        end
    end
end

# bulls を数える
function count_bulls(xs, ys)
    c = 0
    for i in 1 : CSIZE
        if xs[i] == ys[i]; c += 1; end
    end
    c
end

# 同じ数字を数える
function count_same_number(xs, ys)
    c = 0
    for x in xs
        for y in ys
            if x == y
                c += 1
                break
            end
        end
    end
    c
end

function check(answer, xs)
    global query
    for q in query
        b = count_bulls(q.code, xs)
        c = count_same_number(q.code, xs) - b
        if b != q.bulls || c != q.cows
            return
        end
    end
    b = count_bulls(answer, xs)
    c = count_same_number(answer, xs) - b
    q = Query(b, c, xs)
    push!(query, q)
    n = length(query)
    println("$n: $xs, bulls = $b, cows = $c")
    if b == 4
        throw("Good Job!")
    end
end

function solver(answer)
    global query
    query = Query[]
    try
        permutations(xs -> check(answer, xs), collect(0:9))
    catch e
        println(e)
    end
end

Copyright (C) 2016 Makoto Hiroi
All rights reserved.

[ Home | Light | Julia ]