次は 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 が素数ではない場合の対処方法は、高橋謙一郎さんと 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