M.Hiroi's Home Page

Julia Language Programming

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

Puzzle DE Julia!!

[ Home | Light | Julia ]

N 色ライツアウトの解法

次は N 色ライツアウトの解法プログラムを作りましょう。N が素数のとき、その剰余は有限体になるので、とりあえず N は素数に限定します。なお、それ以外の数でも工夫次第で解くことができるようです。興味のある方は下記 URL をお読みくださいませ。有益な情報を公開されている高橋謙一郎さんと deepgreen さんに深く感謝いたします。

GF(N) の場合、足し算や掛け算は簡単です。計算結果に mod N を適用するだけです。ちょっと面倒なのが引き算と割り算です。これらの演算は足し算と掛け算に直して行うことにします。x の負数 (加法の逆元) は N - x で簡単に求めることができます。逆数 (乗法の逆元) は「拡張ユークリッドの互除法」や「フェルマーの小定理」を使うと簡単です。今回は拡張ユークリッドの互除法を使うことにしましょう。詳しい説明は拙作のページ Julia 入門 番外編 拡張ユークリッドの互除法 をお読みください。

●逆数表の作成

それではプログラムを作りましょう。最初に逆数をリストに格納して返す関数 make_inv() を作ります。 

リスト : 逆数表の作成

# 拡張ユークリッドの互除法
function ext_euclid(a::Int, b::Int)
    xs = a, 1, 0
    ys = b, 0, 1
    while ys[1] != 0
        q, z = divrem(xs[1], ys[1])
        xs, ys = ys, (z, xs[2] - q * ys[2], xs[3] - q * ys[3])
    end
    xs
end

# 逆数表の生成
function make_inv(n)
    inv = zeros(Int, n)
    for i = 1 : n - 1
        g, x, _ = ext_euclid(i, n)
        if g != 1 error("色数は素数にしてください") end
        inv[i] = (x + n) % n
    end
    inv
end

関数 ext_euclid() は拡張ユークリッドの互除法で a * x + b * y = gcd(a, b) を満たす (x, y) の組を求めます。返り値は gcd(a, b), x, y です。make_inv(n) の引数は色の個数で、1 から n - 1 までの逆数をリストに格納して返します。n が素数ではない場合、ext_euclid() の返り値 g が 1 にならない数があります。この数には逆数がないのでエラーを送出します。逆数は x にセットされますが、負数が返ってくる場合があるので、(x + n) % n を計算して inv[i] にセットします。

●枢軸の選択

次は枢軸を選択する関数 select_pivot() を作ります。 

リスト : 枢軸の選択

function select_pivot(xs::Matrix{Int}, i::Int)
    # 0 以外のものと交換する
    for j = i : size(xs, 1)
        if xs[j, i] != 0
            if j != i
                temp = xs[i, :]
                xs[i, :] = xs[j, :]
                xs[j, :] = temp
            end
            break
        end
    end
end

今回は整数の演算なので、最大値を選択する必要はありません。xs[j, i] が 0 ではない行 j と行 i を交換するだけです。

●解法プログラムの作成

最後に、ガウスの消去法で解を求める関数 nlo() を作ります。 

リスト : N 色ライツアウトの解法 (N は素数)

function nlo(ys::Matrix{Int}, c::Int)
    ans = []
    inv = make_inv(c)
    zs = make_matrix((c .- ys) .% c)  # 負数は c - x
    n = size(zs, 1)
    m = n + 1
    for i = 1 : n
        select_pivot(zs, i)
        if zs[i, i] == 0
            if all(isequal(0), zs[i:end, i:end])
                m = i
                break    # 複数の解がある
            end
            return ans   # 解なし
        end
        for j = i + 1 : n
            if zs[j, i] != 0
                temp = zs[j, i] * inv[zs[i, i]]
                zs[j, i:end] += c .- ((temp * zs[i, i:end]) .% c)
                zs[j, i:end] .%= c
                if zs[j, i] != 0 error("係数を 0 にできないよ") end
            end
        end
    end
    print_matrix(zs)    # debug

    temp = zs[:, end]
    for j = 0 : c ^ (n + 1 - m) - 1
        if n + 1 - m > 0
            zs[m:end, end] = num_conv(j, n + 1 - m, c)
        end
        # 後退代入
        for k = m - 1 : -1 : 1
            zs[k, end] += c - (dot(zs[k, k+1:n], zs[k+1:end, end]) % c)
            zs[k, end] *= inv[zs[k, k]]
            zs[k, end] %= c
        end
        push!(ans, reshape(zs[:, end], size(ys)))
        zs[:, end] = temp
    end
    ans
end

引数 ys が盤面の状態、c が色数を表します。最初に make_inv() で逆数表を求めて、それを変数 inv にセットします。make_matrix() で拡大係数行列を生成するとき、引数 ys の負数を渡すことに注意してください。

前進消去の計算で、浮動小数点数の場合は temp = zs[j, i] / zs[i, i] としましたが、GF(C) の世界では temp = zs[j, i] * inv[zs[i, i]] とします。zs[j, i:end] から引き算する処理も、(temp * zs[i, i:end]) .% c の負数を求めて、それを加算する処理になります。加算と乗算のあとは剰余を計算することをお忘れなく。

m < n + 1 の場合、c ** (n + 1 - m) 個の解があります。for ループですべての解を出力します。関数 num_conv() は整数 j を n + 1 - m 桁の c 進数に変換し、その結果をリストに格納して返します。その値を zs[m:end, end] にセットしてから後退代入の処理を行います。ここでも、引き算を足し算に、割り算を掛け算に直して計算しています。

あとのプログラムは簡単なので説明は割愛します。詳細は プログラムリスト１ をお読みください。

●実行結果

簡単なテストと実行結果を示します。ご参考までに、行の基本変形後の拡大係数行列を表示しています。 

julia> printanswer(nlo(ones(Int, 5, 5), 2))
1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 1 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 1 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 1 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 1 1 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
No 1, 手数 = 15
0 0 0 1 1
1 1 0 1 1
1 1 1 0 0
0 1 1 1 0
1 0 1 1 0

No 2, 手数 = 15
0 1 1 0 1
0 1 1 1 0
0 0 1 1 1
1 1 0 1 1
1 1 0 0 0

No 3, 手数 = 15
1 1 0 0 0
1 1 0 1 1
0 0 1 1 1
0 1 1 1 0
0 1 1 0 1

No 4, 手数 = 15
1 0 1 1 0
0 1 1 1 0
1 1 1 0 0
1 1 0 1 1
0 0 0 1 1

julia> printanswer(nlo(ones(Int, 5, 5), 3))
1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2
0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2
0 0 1 0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 1 1 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2
0 0 0 0 2 0 1 1 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 2 1 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 2
0 0 0 0 0 0 0 1 2 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 2
0 0 0 0 0 0 0 0 1 1 2 2 2 0 0 1 0 0 0 0 0 0 0 0 0 2
0 0 0 0 0 0 0 0 0 1 1 1 2 0 0 1 1 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 2 1 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 1 2 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 2 1 1 1 2 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 2 0 1 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 2 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 2 2 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 2 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 0 2 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 2 0 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
No 1, 手数 = 24
1 1 1 0 1
0 2 0 0 1
2 2 2 1 0
1 0 0 2 0
2 2 1 2 0

No 2, 手数 = 27
2 2 0 2 2
1 1 1 1 1
1 0 2 0 1
0 1 2 1 0
0 2 2 2 0

No 3, 手数 = 27
0 0 2 1 0
2 0 2 2 1
0 1 2 2 2
2 2 1 0 0
1 2 0 2 0

・・・ 省略 ・・・

No 25, 手数 = 27
0 0 1 1 2
2 1 0 1 2
2 2 2 1 0
2 1 0 1 2
0 0 1 1 2

No 26, 手数 = 24
1 1 0 0 0
0 0 1 2 2
1 0 2 0 1
1 2 2 0 2
1 0 2 1 2

No 27, 手数 = 33
2 2 2 2 1
1 2 2 0 2
0 1 2 2 2
0 0 1 2 2
2 0 0 1 2

julia> printanswer(nlo(ones(Int, 5, 5), 5))
1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4
0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4
0 0 1 0 0 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 1 1 4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4
0 0 0 0 4 0 1 1 4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 4 1 0 4 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 3 1 4 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 4 2 4 3 3 1 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 4 0 2 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 2 4 2 2 0 1 0 0 0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 0 0 0 1 0 4 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 2 4 3 1 0 1 0 0 0 0 0 0 0 0 4
0 0 0 0 0 0 0 0 0 0 0 0 4 4 4 1 3 1 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 4 3 4 1 0 0 0 0 0 0 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 1 2 0 1 0 1 0 0 0 0 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 1 0 4 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 2 0 3 1 1 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 2 3 1 1 1 0 0 3
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 4 0 1 1 1 0 3
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 2 2 3 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 4 0 1 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 4 0 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
No 1, 手数 = 51
2 4 4 4 2
3 4 2 4 3
0 1 0 1 0
0 4 0 4 0
0 4 1 4 0

No 2, 手数 = 56
2 0 3 3 4
2 4 3 4 2
1 0 0 2 4
1 4 4 4 1
3 0 2 3 0

No 3, 手数 = 56
2 1 2 2 1
1 4 4 4 1
2 4 0 3 3
2 4 3 4 2
1 1 3 2 0

・・・ 省略 ・・・

No 23, 手数 = 56
3 0 2 3 0
1 4 4 4 1
1 0 0 2 4
2 4 3 4 2
2 0 3 3 4

No 24, 手数 = 56
3 1 1 2 2
0 4 0 4 0
2 4 0 3 3
3 4 2 4 3
0 1 4 2 4

No 25, 手数 = 66
3 2 0 1 4
4 4 1 4 4
3 3 0 4 2
4 4 1 4 4
3 2 0 1 4

5 * 5 盤で連立方程式が解ける場合、3 色では 27 通り、5 色では 25 通りの解が表示されました。行の基本変形後の係数行列の自由度 (行数 - 階数) を n, 色数を c とすると、解の総数は cn 通りになります。自由度を調べたところ、次のようになりました。 

         表 : 自由度

       |  2  3  5  7 11 13 17 19
-------+-------------------------
(3, 3) |  0  0  0  1  0  0  0  0
(4, 4) |  4  2  2  2  2  2  2  2
(5, 5) |  2  3  2  2  3  2  2  2
(6, 6) |  0  0  0  0  0  3  0  0
(7, 7) |  0  0  0  1  0  0  3  0
(8, 8) |  0  4  0  0  0  0  2  1
(9, 9) |  8  2  2  2  2  2  2  4

自由度 0 はどんな点灯パターンでも解けることを意味します。たとえば、3 * 3 盤は 7 色以外の自由度は 0 なので、ボタンが全部同じ色のパターンは解くことができます。ところが 7 色の場合、ボタンが全部同じ色のパターンは解くことができません。次に示す点灯パターンは解くことができて、7 通りの解が表示されます。 

julia> printanswer(nlo([6 6 6; 6 6 6; 6 0 0], 7))
1 1 0 1 0 0 0 0 0 1
0 1 1 0 0 1 0 0 0 1
0 0 1 6 1 0 0 0 0 0
0 0 0 1 0 1 1 0 0 1
0 0 0 0 2 0 0 1 0 0
0 0 0 0 0 6 0 1 0 0
0 0 0 0 0 0 6 0 1 6
0 0 0 0 0 0 0 4 2 0
0 0 0 0 0 0 0 0 0 0
No 1, 手数 = 2
0 0 1
1 0 0
0 0 0

No 2, 手数 = 20
1 3 2
4 2 3
1 3 1

No 3, 手数 = 31
2 6 3
0 4 6
2 6 2

No 4, 手数 = 28
3 2 4
3 6 2
3 2 3

No 5, 手数 = 39
4 5 5
6 1 5
4 5 4

No 6, 手数 = 29
5 1 6
2 3 1
5 1 5

No 7, 手数 = 40
6 4 0
5 5 4
6 4 6

盤面を大きくしたり色数を増やすと、ボタンを押す回数がとても多くなります。ゲームとして遊ぶには、5 * 5 盤の 2, 3, 5 色程度が手ごろなのかもしれませんね。

●N が素数以外の場合

N が素数ではない場合の対処方法は、高橋謙一郎さんと deepgreen さんが考案されていますが、今回は deepgreen さんの方法を使わせてもらうことにします。基本的な考え方はとてもシンプルです。係数行列 zs の対角成分 zs[i, i] において、その値の逆数が存在しないとき、逆数が存在する行 j と交換する、という方法です。つまり、割り算できる値を選択するわけです。ただし、この方法ですべて上手くいくわけではなく、対角成分に逆数が存在しない値が残る場合があります。簡単な例を示しましょう。 

----- 4 -----
[ 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 1 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;

        ・・・ 省略 ・・・

  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 3 3 2;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 1;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
----- 6 -----
[ 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 1 0 0 5 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;

        ・・・ 省略 ・・・

  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 3 1 0 2;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 3 2 3;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 3 3 ]

5 * 5 盤 4 色の係数行列は今までと同様に変形できましたが、6 色の係数行列は対角成分に 3 が残ってしまいます。これらの方程式は解の条件を示しています。右辺値を d1, d2, d3 とすると、3 * (x23 + x24 + x25) = di (mod 6, i = 1, 2, 3) を満たす x23, x24, x25 を求めればいいわけです。

これらの方程式を解く賢い方法があるのかもしれませんが、今回は単純な力任せで探索することにします。具体的には、[0, 0, 0] から [5, 5, 5] までのベクタ (63 = 216 通り) を生成します。あとは、左辺式を計算して、その値が右辺値と等しいことをチェックするだけです。

●逆数表とピボット選択の修正

それではプログラムを作りましょう。まずは枢軸を選択する関数 select_pivot() を修正します。 

リスト : 逆数表とピボット選択

# 逆数表の生成
function make_inv(n)
    inv = zeros(Int, n)
    for i = 1 : n - 1
        g, x, _ = ext_euclid(i, n)
        if g == 1
            inv[i] = (x + n) % n
        end
    end
    inv
end

# ピボット選択
function select_pivot(xs::Matrix{Int}, i::Int, inv::Vector{Int})
    for j = i : size(xs, 1)
        if xs[j, i] != 0 && inv[xs[j, i]] != 0
            if j != i
                temp = xs[i, :]
                xs[i, :] = xs[j, :]
                xs[j, :] = temp
            end
            break
        end
    end
end

make_inv() では、g が 1 のときだけ逆数をセットします。逆数表が 0 の場合、その数の逆数は存在しません。select_pivot() は引数 inv に逆数表を渡します。inv[xs[j, i]] が 0 でなければ逆数が存在するので行 i と j を交換します。

●解法プログラムの修正

最後に、解を求める関数 nlo() を修正します。 

リスト : N 色ライツアウトの解法

function nlo(ys::Matrix{Int}, c::Int)
    ans = []
    inv = make_inv(c)
    zs = make_matrix((c .- ys) .% c)  # 負数は c - x
    n = size(zs, 1)
    m = n + 1
    check = false
    for i = 1 : n
        select_pivot(zs, i, inv)
        if zs[i, i] == 0 || inv[zs[i, i]] == 0
            if all(isequal(0), zs[i:end, i:n])
                if any(x -> x > 0, zs[i:end, end])
                    return ans  # 解なし
                end
            else
                check = true
            end
            m = i
            break    # 複数の解がある
        end
        for j = i + 1 : n
            if zs[j, i] != 0
                temp = zs[j, i] * inv[zs[i, i]]
                zs[j, i:end] += c .- ((temp * zs[i, i:end]) .% c)
                zs[j, i:end] .%= c
                if zs[j, i] != 0 error("係数を 0 にできないよ") end
            end
        end
    end
    # print_matrix(zs)    # debug

    temp = zs[:, end]
    for j = 0 : c ^ (n + 1 - m) - 1
        if n + 1 - m > 0
            vs = num_conv(j, n + 1 - m, c)
            if check && !((zs[m:end, m:n] * vs) .% c == zs[m:end, end]) continue end
            zs[m:end, end] = vs
        end
        # 後退代入
        for k = m - 1 : -1 : 1
            zs[k, end] += c - (dot(zs[k, k+1:n], zs[k+1:end, end]) % c)
            zs[k, end] *= inv[zs[k, k]]
            zs[k, end] %= c
        end
        push!(ans, reshape(zs[:, end], size(ys)))
        zs[:, end] = temp
    end
    ans
end

select_pivot() で行を選択したあと zs[i, i] と inv[zs[i, i]] をチェックし、値が 0 ならば前進消去を終了します。係数行列 zs[i:end, i:n] の要素が全部 0 ならば、今までと同じパターンです。右辺値 zs[i:end, end] をチェックして 0 以外の値がある場合、連立方程式に解はありません。

zs[i:end, i:n] の要素に 0 以外の値がある場合、残っている方程式を満たす解を探索します。変数 check を true にセットして for ループを脱出します。なお、変数 check は false に初期化しておきます。

後退代入では、num_conv() でベクトル vs を生成するところは同じです。check が true ならば、左辺式を計算して右辺値と等しいかチェックします。左辺式の計算は簡単で、zs[m:end, m:n] と vs の積を計算し、その値の剰余を求めるだけです。この値が右辺値 zs[m:end, end] と等しければ条件を満たしています。そうでなければ continue でループの先頭に戻り、次の値を試します。

あとのプログラムは簡単なので説明は割愛します。詳細は プログラムリスト２ をお読みくださいませ。

●実行結果

それでは実際に試してみましょう。5 * 5 盤で全点灯パターンの最小手数を求めます。 

リスト : 簡単なテスト

function test()
    for n = 2 : 9
        println("----- $n -----")
        m = n * 25
        xs = nlo(fill(n - 1, 5, 5), n)
        ys = map(zs -> sum(zs), xs)
        z, i = findmin(ys)
        print_matrix(xs[i])
        println(z, ", ", length(xs))
    end
end

最小手数だけではなく、解の総数も表示します。結果は次のようになりました。 

julia> test()
----- 2 -----
0 0 0 1 1
1 1 0 1 1
1 1 1 0 0
0 1 1 1 0
1 0 1 1 0
15, 4
----- 3 -----
0 0 1 2 0
1 0 1 1 2
0 2 1 1 1
1 1 2 0 0
2 1 0 1 0
21, 27
----- 4 -----
0 0 2 1 3
1 3 2 3 1
1 3 3 2 2
0 3 3 3 0
1 0 1 1 0
39, 16
----- 5 -----
3 1 1 1 3
2 1 3 1 2
0 4 0 4 0
0 1 0 1 0
0 1 4 1 0
34, 25
----- 6 -----
0 3 1 2 3
4 3 1 1 2
0 2 1 1 1
1 1 2 3 3
5 1 0 4 0
45, 108
----- 7 -----
1 0 5 0 1
0 2 3 2 0
5 3 3 3 5
0 2 3 2 0
1 0 5 0 1
47, 49
----- 8 -----
0 4 2 1 3
5 3 2 3 5
1 3 7 6 6
0 3 7 3 0
5 4 5 1 0
79, 64
----- 9 -----
1 4 6 4 7
5 8 5 2 8
5 6 7 0 2
3 2 1 8 0
0 7 1 1 0
93, 243

ところで、N 色ライツアウトで全点灯パターンが必ず解けるわけではありません。実際、3 * 3 盤 7 色, 5 * 5 盤 11 色, 8 * 8 盤 3, 6, 9 色の全点灯パターンは解くことができませんでした。なお、8 * 8 盤 6 色の場合、8 本の方程式が残るので 0 から 6 ** 8 - 1 = 1679615 までの値をチェックすることになります。Python3 で試したときはけっこう時間がかかったのですが、Julia では 1 秒かかりませんでした。やっぱり Julia は速いですね。

そうはいっても、単純な力任せの探索では限界があります。残った方程式の数が多くなると、実用的な時間では解けないことがあります。Julia の場合、7 * 7 盤 21 色までは高速に解くことができました。8 * 8 盤で色数を増やすと、このプログラムでは時間がかかる場合があります。盤の大きさや色数を増やすときはご注意ください。ご参考までに、全点灯パターンの解の総数を下表に示します。 

                表 : 全点灯パターンの解の総数

       |  2  3  4  5   6  7   8   9  10  11  12  13  14  15   16  17  18   19   20   21
-------+--------------------------------------------------------------------------------  
(3, 3) |  1  1  1  1   1  0   1   1   1   1   1   1   0   1    1   1   1    1    1    0
(4, 4) | 16  9 64 25 144 49 256  81 400 121 576 169 784 225 1024 289 1296 361 1600  441
(5, 5) |  4 27 16 25 108 49  64 243 100   0 432 169 196 675  256 289  972 361  400 1323
(6, 6) |  1  1  1  1   1  1   1   1   1   1   1   0   1   1    1   1   1    1    1    1
(7, 7) |  1  1  1  1   1  7   1   1   1   1   1   1   7   1    1   0   1    1    1    7
(8, 8) |  1  0  1  1   0  1   1   0   1   1  NG   1   1   0    1 289  NG   19    1    0

8 * 8 盤の NG は、「完全ピボット選択」を行うと実用的な時間で解くことが可能です。興味のある方は拙作のページ NumPy 超入門 N 色ライツアウトの解法 の追記をお読みくださいませ。

●プログラムリスト１

#
# nlo.jl : N 色ライツアウトの解法 (N は素数)
#
#          Copyright (C) 2018-2021 Makoto Hiroi
#
using LinearAlgebra

# 拡張ユークリッドの互除法
function ext_euclid(a::Int, b::Int)
    xs = a, 1, 0
    ys = b, 0, 1
    while ys[1] != 0
        q, z = divrem(xs[1], ys[1])
        xs, ys = ys, (z, xs[2] - q * ys[2], xs[3] - q * ys[3])
    end
    xs
end

# 逆数表の生成
function make_inv(n)
    inv = zeros(Int, n)
    for i = 1 : n - 1
        g, x, _ = ext_euclid(i, n)
        if g != 1 error("色数は素数にしてください") end
        inv[i] = (x + n) % n
    end
    inv
end

# 係数行列の作成
function make_matrix(ys::Matrix{Int})
    n, m = size(ys)
    k = n * m
    xs = zeros(Int, k, k + 1)
    for y = 1 : n, x = 1 : m
        z = (x - 1) * n + y
        for (dx, dy) = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
            y1 = y + dy
            x1 = x + dx
            if 0 < x1 <= m && 0 < y1 <= n
                xs[z, (x1 - 1) * n + y1] = 1
            end
        end
    end
    xs[:, end] = reshape(ys, k, 1)
    xs
end

# ピボット選択
function select_pivot(xs::Matrix{Int}, i::Int)
    # 0 以外のものと交換する
    for j = i : size(xs, 1)
        if xs[j, i] != 0
            if j != i
                temp = xs[i, :]
                xs[i, :] = xs[j, :]
                xs[j, :] = temp
            end
            break
        end
    end
end

# 数値 m を n 桁 c 進数に変換
function num_conv(m::Int, n::Int, c::Int)
    buff = zeros(Int, n)
    i = 1
    while m > 0
        buff[i] = m % c
        m = div(m, c)
        i += 1
    end
    buff
end

# 行列の表示
function print_matrix(xs::Matrix{Int})
    n, m = size(xs)
    for i = 1 : n
        for j = 1 : m
            print("$(xs[i, j]) ")
        end
        println("")
    end
end

# 解の表示
function printanswer(ans)
    for (i, a) = enumerate(ans)
        println("No $i, 手数 = $(sum(a))")
        print_matrix(a)
        println("")
    end
end

# ライツアウトの解法
function nlo(ys::Matrix{Int}, c::Int)
    ans = []
    inv = make_inv(c)
    zs = make_matrix((c .- ys) .% c)  # 負数は c - x
    n = size(zs, 1)
    m = n + 1
    for i = 1 : n
        select_pivot(zs, i)
        if zs[i, i] == 0
            if all(isequal(0), zs[i:end, i:end])
                m = i
                break    # 複数の解がある
            end
            return ans   # 解なし
        end
        for j = i + 1 : n
            if zs[j, i] != 0
                temp = zs[j, i] * inv[zs[i, i]]
                zs[j, i:end] += c .- ((temp * zs[i, i:end]) .% c)
                zs[j, i:end] .%= c
                if zs[j, i] != 0 error("係数を 0 にできないよ") end
            end
        end
    end
    print_matrix(zs)    # debug

    temp = zs[:, end]
    for j = 0 : c ^ (n + 1 - m) - 1
        if n + 1 - m > 0
            zs[m:end, end] = num_conv(j, n + 1 - m, c)
        end
        # 後退代入
        for k = m - 1 : -1 : 1
            zs[k, end] += c - (dot(zs[k, k+1:n], zs[k+1:end, end]) % c)
            zs[k, end] *= inv[zs[k, k]]
            zs[k, end] %= c
        end
        push!(ans, reshape(zs[:, end], size(ys)))
        zs[:, end] = temp
    end
    ans
end

●プログラムリスト２

#
# nlo1.jl : N 色ライツアウトの解法 (N は素数でなくてもよい)
#
#           Copyright (C) 2018-2021 Makoto Hiroi
#
using LinearAlgebra

# 拡張ユークリッドの互除法
function ext_euclid(a::Int, b::Int)
    xs = a, 1, 0
    ys = b, 0, 1
    while ys[1] != 0
        q, z = divrem(xs[1], ys[1])
        xs, ys = ys, (z, xs[2] - q * ys[2], xs[3] - q * ys[3])
    end
    xs
end

# 逆数表の生成
function make_inv(n)
    inv = zeros(Int, n)
    for i = 1 : n - 1
        g, x, _ = ext_euclid(i, n)
        if g == 1
            inv[i] = (x + n) % n
        end
    end
    inv
end

# 係数行列の作成
function make_matrix(ys::Matrix{Int})
    n, m = size(ys)
    k = n * m
    xs = zeros(Int, k, k + 1)
    for y = 1 : n, x = 1 : m
        z = (x - 1) * n + y
        for (dx, dy) = [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]
            y1 = y + dy
            x1 = x + dx
            if 0 < x1 <= m && 0 < y1 <= n
                xs[z, (x1 - 1) * n + y1] = 1
            end
        end
    end
    xs[:, end] = reshape(ys, k, 1)
    xs
end

# ピボット選択
function select_pivot(xs::Matrix{Int}, i::Int, inv::Vector{Int})
    # 0 以外のものと交換する
    for j = i : size(xs, 1)
        if xs[j, i] != 0 && inv[xs[j, i]] != 0
            if j != i
                temp = xs[i, :]
                xs[i, :] = xs[j, :]
                xs[j, :] = temp
            end
            break
        end
    end
end

# 数値 m を n 桁 c 進数に変換
function num_conv(m::Int, n::Int, c::Int)
    buff = zeros(Int, n)
    i = 1
    while m > 0
        buff[i] = m % c
        m = div(m, c)
        i += 1
    end
    buff
end

# 行列の表示
function print_matrix(xs::Matrix{Int})
    n, m = size(xs)
    for i = 1 : n
        for j = 1 : m
            print("$(xs[i, j]) ")
        end
        println("")
    end
end


# 解の表示
function printanswer(ans)
    for (i, a) = enumerate(ans)
        println("No $i, 手数 = $(sum(a))")
        print_matrix(a)
        println("")
    end
end

# ライツアウトの解法
function nlo(ys::Matrix{Int}, c::Int)
    ans = []
    inv = make_inv(c)
    zs = make_matrix((c .- ys) .% c)  # 負数は c - x
    n = size(zs, 1)
    m = n + 1
    check = false
    for i = 1 : n
        select_pivot(zs, i, inv)
        if zs[i, i] == 0 || inv[zs[i, i]] == 0
            if all(isequal(0), zs[i:end, i:n])
                if any(x -> x > 0, zs[i:end, end])
                    return ans  # 解なし
                end
            else
                check = true
            end
            m = i
            break    # 複数の解がある
        end
        for j = i + 1 : n
            if zs[j, i] != 0
                temp = zs[j, i] * inv[zs[i, i]]
                zs[j, i:end] += c .- ((temp * zs[i, i:end]) .% c)
                zs[j, i:end] .%= c
                if zs[j, i] != 0 error("係数を 0 にできないよ") end
            end
        end
    end
    # print_matrix(zs)    # debug

    temp = zs[:, end]
    for j = 0 : c ^ (n + 1 - m) - 1
        if n + 1 - m > 0
            vs = num_conv(j, n + 1 - m, c)
            if check && !((zs[m:end, m:n] * vs) .% c == zs[m:end, end]) continue end
            zs[m:end, end] = vs
        end
        # 後退代入
        for k = m - 1 : -1 : 1
            zs[k, end] += c - (dot(zs[k, k+1:n], zs[k+1:end, end]) % c)
            zs[k, end] *= inv[zs[k, k]]
            zs[k, end] %= c
        end
        push!(ans, reshape(zs[:, end], size(ys)))
        zs[:, end] = temp
    end
    ans
end
2021/12/05 改訂: Julia のバージョンを 1.6 に変更
Copyright (C) 2018-2021 Makoto Hiroi
All rights reserved.
[ Home | Light | Julia ]