次は 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] にセットしてから後退代入の処理を行います。ここでも、引き算を足し算に、割り算を掛け算に直して計算しています。
あとのプログラムは簡単なので説明は割愛します。詳細はプログラムリスト1をお読みください。
簡単なテストと実行結果を示します。ご参考までに、行の基本変形後の拡大係数行列を表示しています。
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 でループの先頭に戻り、次の値を試します。
あとのプログラムは簡単なので説明は割愛します。詳細はプログラムリスト2をお読みくださいませ。
それでは実際に試してみましょう。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