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

●連立一次方程式の解法 (2)

●行列式

NumPy には行列 \(A\) の行列式 (determinant) を求める関数 numpy.linalg.det(A) が用意されていますが、\(A\) を LU 分解するときに計算することもできます。次に示す行列式と三角行列の性質を使います。以下では \(A\) の行列式を \(|A|\) と表しています。

n 次の正方行列 \(A, B\) において、\(|AB| = |A||B|\) が成り立つ
行列 \(A\) の行 (または列) を交換した行列 \(A'\) において、\(|A'| = -|A|\) が成り立つ
三角行列の行列式は対角成分の積になる

ここでは数学的な証明は行わずに簡単な例を示すだけにとどめます。

(1)
  [[a, b],  @  [[e, f],  =  [[ae + bg, af + bh],
   [c, d]]      [g, h]]      [ce + dg, cf + dh]]

  右辺の行列式の乗算 (ad - bc) * (eh - fg) = adeh - adfg - bceh + bcfg
  左辺の行列式 (ae + bg)(cf + dh) - (af + bh)(ce + dg)
  = acef + adeh + bcfg + bdgh - acef - adfg - bceh - bdgh
    ^^^^                 ++++   ^^^^                 ++++
  = adeh - adfg - bceh + bcfg

(2)
  [[a, b],    [[c, d],   
   [c, d]]     [a, b]]

   ad - bc     cb - da = - (ad - bc)

(3)
 [[a, b],       [[a, 0],
  [0, d]]        [c, d]]

ad - b0 = ad   ad - 0c = ad

 [[a, b, c],
  [d, e, f],    => aei + bfg + chd - ceg - bdi - ahf
  [g, h, i]],

上三角行列は d, h, g が 0 なので行列式は aei になる
下三角行列は b, c, f が 0 なので行列式は aei になる

前回作成したプログラムでは行列 \(A\) を LU に分解したとき、\(L\) の対角成分が 1 になるので、\(U\) の対角成分の積が \(A\) の行列式になります。ピボット選択を行う場合、行列式の符号を記憶しておいて、行を交換したときは符号を反転すればいいでしょう。

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

リスト : 行列式の計算

import numpy as np

# ピボット選択
def select_pivot(xs, idx, i):
    k = np.abs(xs[i:,i]).argmax() + i
    if k != i:
        temp = xs[i].copy()
        xs[i] = xs[k]
        xs[k] = temp
        idx[i], idx[k] = idx[k], idx[i]
        return -1
    return 1

# ピボット選択付き
def lu_pv(xs):
    n = len(xs)
    zs = xs.astype(np.float_)
    idx = list(range(n))
    det = 1
    for i in range(n):
        det *= select_pivot(zs, idx, i)
        if zs[i, i] == 0: break
        for j in range(i + 1, n):
            temp = zs[j, i] / zs[i, i]
            zs[j, i+1:] -= temp * zs[i, i+1:]
            zs[j, i] = temp
    return det * np.diag(zs).prod(), zs, idx

select_pivot() は行を交換したら -1 を、交換しなければ 1 を返すように変更します。lu_pv() は行列式、LU 分解した行列、行の交換を記録したリストを返します。行列式の符号は変数 det にセットし、select_pivot() の返り値と乗算します。これで行を交換したら符号を反転することができます。対角成分 zs[i, i] が 0 ならば LU 分解できないので、break でループを脱出します。行列式の計算は関数 diag() で zs の対角成分を取り出し、prod() で要素の乗算を求めるだけです。対角成分に 0 の要素があるならば行列式は 0 になります。

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

リスト : 簡単なテスト

a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
a5 = np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]])
a6 = np.array([[2, 4, 2, 2], [4, 10, 3, 3], [2, 6, 1, 1], [3, 7, 1, 4]])
for a in [a1, a2, a3, a4, a5, a6]:
    det, xs, idx = lu_pv(a)
    print(det)
    print(np.linalg.det(a))
    print(xs)
    print(idx)

2.0
2.0
[[ 2.   4. ]
 [ 0.5 -1. ]]
[1, 0]

12.0
12.0
[[ 2.    4.    6.  ]
 [ 1.   -4.   -2.  ]
 [ 0.5   0.25 -1.5 ]]
[1, 2, 0]

72.0
72.0
[[ 8.      4.      2.      1.    ]
 [-1.      8.      0.      2.    ]
 [ 0.125   0.0625  0.75    0.75  ]
 [-0.125   0.1875 -1.      1.5   ]]
[2, 3, 0, 1]

384.0
384.0
[[ 12.          -6.           2.           0.           0.        ]
 [  1.          12.           0.           0.           0.        ]
 [  0.33333333   0.41666667   1.33333333   1.           0.        ]
 [  0.08333333  -0.04166667   0.625       -1.625        1.        ]
 [  0.08333333   0.125        0.625       -0.23076923   1.23076923]]
[1, 3, 4, 0, 2]

-12.0
-12.0
[[ 4.    2.    6.  ]
 [ 0.    2.    4.  ]
 [ 0.25  0.25 -1.5 ]]
[2, 0, 1]

0.0
0.0
[[  4.    10.     3.     3.  ]
 [  0.5   -1.     0.5    0.5 ]
 [  0.75   0.5   -1.5    1.5 ]
 [  0.5   -1.    -0.     0.  ]]
[1, 0, 3, 2]

最後は行列式が 0 になる場合です。この場合、解を一意的に定めることはできません。解は無数にあるか存在しないかのどちらかになります。

●解の存在条件

連立一次方程式の解の存在条件は行列の「階数 (rank)」を使って判定することができます。階数の定義はいくつか方法がありますが、ここでは行列の基本変形を使うことにしましょう。以下に示す操作を行の基本変形といいます。

行列 \(A\) の i 行と j 行を交換する
行列 \(A\) の i 行を c 倍する (\(c \ne 0\))
行列 \(A\) の i 行に j 行の c 倍を加算する (\(c \ne 0\))

同様に列の基本変形を定義することができます。これらの基本変形は、ガウスの消去法やガウス・ジョルダン法で用いる操作と同じです。一般に、任意の行列は基本変形を何回か適用すると「階段行列」に変形することが可能です。階段行列は左下半分に 0 が階段状に並んだ行列のことです。

 [[a, b, c, d],    [[a, b, c, d],    [[a, b, c, d],
  [0, e, f, g],     [0, e, f, g],     [0, e, f, g],
  [0, 0, h, i]]     [0, 0, 0, h]]     [0, 0, 0, 0]]

   (1) rank = 3     (2) rank = 3      (3) rank = 2

                      階段行列

階段行列に変形したあと、零ベクトル (要素がすべて 0 のベクトル) を除いた行ベクトルの本数が階数になります。上図でいうと、(1) と (2) の行列には零ベクトルはないので階数は 3 になりますが、(3) の行列は零ベクトルが一つあるので階数は 2 になります。

ここで、上図の行列が拡大係数行列を表しているとしましょう。(1) の場合、係数行列と拡大係数行列の階数はどちらも 3 になり、拡大係数行列の行数と一致します。この場合、連立一次方程式の解は一意的に定まります。(2) の場合、3 行 3 列の係数行列としてみると、その階数は 2 になり、拡大係数行列の階数よりも少なくなります。この場合は解がありません。係数がすべて 0 なのに定数が残っているので、連立方程式として矛盾しているわけです。

(3) の場合、係数行列と拡大係数行列の階数はともに 2 になりますが、拡大係数行列の行数とは一致しません。このとき、零ベクトルに対応する変数に適当な値を代入すると、残りの変数を決定することができます。つまり、連立一次方程式の解は存在するが、一意的には定めることができないわけです。まとめると次のようになります。

係数行列 \(A\) と拡大係数行列 \((Ab)\) の階数が等しいとき、連立方程式 \(Ax = b\) の解は存在する
\(A\) と \((Ab)\) の階数がその行数と等しいとき、\(Ax = b\) の解は一意的に定まる

簡単な例を示しましょう。

(1)
2a +  4b + 2c + 2d = 8
4a + 10b + 3c + 3d = 17
2a +  6b +  c +  d = 9
3a +  7b +  c + 4d = 11

行の基本変形
[[  4.    10.     3.     3.    17.  ]
 [  0.    -1.     0.5    0.5   -0.5 ]
 [  0.     1.    -0.5   -0.5    0.5 ]
 [  0.    -0.5   -1.25   1.75  -1.75]]
[[  4.   10.    3.    3.   17. ]
 [  0.   -1.    0.5   0.5  -0.5]
 [  0.    0.    0.    0.    0. ]
 [  0.    0.   -1.5   1.5  -1.5]]
[[  4.   10.    3.    3.   17. ]
 [  0.   -1.    0.5   0.5  -0.5]  
 [  0.    0.   -1.5   1.5  -1.5]   
 [  0.    0.    0.    0.    0. ]]  

係数行列と拡大係数行列の階数は 3

d = k とおくと
c = (-1.5 - 1.5d) / -1.5 = k + 1
b = (-0.5 - 0.5c - 0.5d) / -1 = k + 1
a = (17 - 10b - 3c - 3d) / 4 = -4k + 1

k は任意の定数なので解は無数にある

(2)
2a +  4b + 2c + 2d = 8
4a + 10b + 3c + 3d = 17
2a +  6b +  c +  d = 10
3a +  7b +  c + 4d = 11

行の基本変形
[[  4.    10.     3.     3.    17.  ]
 [  0.    -1.     0.5    0.5   -0.5 ]
 [  0.     1.    -0.5   -0.5    1.5 ]
 [  0.    -0.5   -1.25   1.75  -1.75]]
[[  4.   10.    3.    3.   17. ]
 [  0.   -1.    0.5   0.5  -0.5]
 [  0.    0.    0.    0.    1. ]
 [  0.    0.   -1.5   1.5  -1.5]]
[[  4.   10.    3.    3.   17. ]
 [  0.   -1.    0.5   0.5  -0.5]
 [  0.    0.   -1.5   1.5  -1.5]
 [  0.    0.    0.    0.    1. ]]

係数行列の階数は 3 で拡大係数行列の階数は 4 なので解はない

(3)
2a +  4b + 2c + 2d = 8
4a + 10b + 3c + 3d = 17
2a +  6b +  c + 2d = 9
3a +  7b +  c + 4d = 11

行の基本変形
[[  4.    10.     3.     3.    17.  ]
 [  0.    -1.     0.5    0.5   -0.5 ]
 [  0.     1.    -0.5    0.5    0.5 ]
 [  0.    -0.5   -1.25   1.75  -1.75]]
[[  4.   10.    3.    3.   17. ]
 [  0.   -1.    0.5   0.5  -0.5]
 [  0.    0.    0.    1.    0. ]
 [  0.    0.   -1.5   1.5  -1.5]]
[[  4.   10.    3.    3.   17. ]
 [  0.   -1.    0.5   0.5  -0.5]
 [  0.    0.   -1.5   1.5  -1.5]
 [  0.    0.    0.    1.    0. ]]

係数行列と拡大係数行列の階数は 4

d = 0 / 1 = 0
c = (-1.5 - 1.5*0) / -1.5 = 1
b = (-0.5 - 0.5*1 - 0.5*0) / -1 = 1
a = (17 - 10*1 - 3*1 - 3*0) / 4 = 1

解は一意的に定まる

●ライツアウトの解法

それでは簡単な例題としてパズル「ライツアウト」を解いてみましょう。パズルの説明は拙作のページ Puzzle DE Programming: ライツアウトの解法をお読みください。

下記 URL によると、ライツアウトの解法は連立一次方程式を解くことに帰着させることができるそうです。

驚異のライツアウト解法ロジック, (高橋謙一郎さん)
ライツアウト＆８めくりパズル, (deepgreen さん)

具体的にいうと、5 行 5 列盤のライツアウトの解は、次に示す 25 本の連立方程式を解くことで求めることができます。

\( (X_{i,j-1} + X_{i-1,j} + X_{i,j} + X_{i+1,j} + X_{i,j+1} + B_{i,j}) \equiv 0 \pmod 2,\quad (0 \lt i \leq 5, 0 \lt j \leq 5) \)

\(X_{i,j}\) はボタン (i, j) を押す回数、\(B_{i,j}\) はボタン (i, j) の状態 (点灯・消灯) を表します。たとえば、ボタン (i, j) が点灯している場合、自身のボタンと周囲のボタンを押す回数が奇数回であれば、それを消灯することができます。逆に、消灯している場合であれば、ボタンを押す回数が偶数回であれば消灯したままになります。これを行列で表すと \(Ax + B \equiv 0 \pmod 2\) となります。\(\equiv\) は合同式を表す記号です

●有限体

数学の世界では、加減乗除ができる数の体系を「体」といいます。特に、数が n 個しかないものを有限体 (またはガロア体) といい GF(n) と表記します。n が素数のとき、その剰余 (\(0, 1, \ldots, n - 1\)) は有限体になります。2 は素数なので GF(2) であり、四則演算ができるのでガウスの消去法により \(Ax + B \equiv 0 \pmod 2\) を解くことができます。

GF(2) は 0 と 1 だけの世界です。加算と乗算は次のように定義できます。

 + | 0 | 1    * | 0 | 1   
---+---+---  ---+---+---  
 0 | 0 | 1    0 | 0 | 0
---+---+---  ---+---+---
 1 | 1 | 0    1 | 0 | 1

GF(n) の場合、加算と乗算の演算結果 m が n 以上になったら m mod n を計算します。GF(2) の場合、加法は排他的論理和 (XOR) に、乗法は論理積 (AND) と同じになります (0 と 1 の掛け算と考えてもかまいません)。

減算の場合、x - y を x + (- y) と考えます。(- y) の加法の逆元といいます。y + (- y) = 0 が成り立つので、0 の逆元は 0 で 1 の逆元は 1 であることがわかります。減算は次のようになります。

 0 - 0 = 0 + (0) = 0
 0 - 1 = 0 + (1) = 1
 1 - 0 = 1 + (0) = 1
 1 - 1 = 1 + (1) = 0

GF(2) の場合、加算と減算は同じ結果 (排他的論理和) になります。したがって、連立方程式 \(Ax + B \equiv 0 \pmod 2\) は \(Ax \equiv B \pmod 2\) を解けばいいことになります。徐算の場合も同様に、x / y を x * (1 / y) として考えて、逆数 \(y^{-1}\) を求めればいいのですが、GF(2) の徐算は簡単で 0 / 1 = 0 と 1 / 1 = 1 になります。

●プログラムの作成

今回は 5 行 5 列盤に限定せずに n 行 m 列盤のライツアウトを解くプログラムを作りましょう。まず最初に、拡大係数行列を生成する関数 make_matrix() を作ります。

リスト : 拡大係数行列の生成

def make_matrix(ys):
    n, m = ys.shape
    k = n * m
    xs = np.zeros((k, k + 1), dtype=np.uint8)
    for y in range(n):
        for x in range(m):
            z = y * m + x
            for dy, dx in [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]:
                y1 = y + dy
                x1 = x + dx
                if 0 <= x1 < m and 0 <= y1 < n:
                    xs[z, y1 * m + x1] = 1
    xs[:, k] = ys.flatten()
    return xs

make_matrix() の引数 ys はボタンの点灯パターンを表す二次元配列です。ys の shape が (n, m) とすると、拡大係数行列の大きさは (n * m, n * m + 1) になります。プログラムでは n * m の値を変数 k にセットしています。次に for ループの中で、ボタン x, y を押したときに反転するボタンの位置を求め、その位置に対応する係数行列の要素を 1 にセットします。最後に、行列 xs の k 列目 (xs[:, k]) に ys を平坦化したものをセットします。これが連立方程式の左辺式になります。

次はライツアウトを解く関数 lo() を作ります。考え方はガウスの消去法と同じです。

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

def lo(ys):
    zs = make_matrix(ys)
    n = len(zs)
    m = n
    # 前進消去
    for i in range(n): # 解の有無をチェックするため、最後の行まで調べる
        select_pivot(zs, i)
        if zs[i, i] == 0:
            if np.all(zs[i:, i:] == 0):
                # 複数の解がある
                m = i
                break
            return    # 解無し
        for j in range(i + 1, n):
            if zs[j, i] == 1: zs[j, i:] ^= zs[i, i:]

    temp = zs[:, n].copy()
    for j in range(0, 2 ** (n - m)):
        for k in range(m, n): zs[k, n] = (j >> (k - m)) & 1
        # 後退代入
        # zs[m:, n] の値は確定済み
        for k in range(m - 1, -1, -1):
            zs[k, n] ^= (zs[k, k+1:n] @ zs[k+1:, n]) % 2
        yield zs[:, n].reshape(ys.shape)
        zs[:, n] = temp

変数 m は階数を表していて n に初期化します。前進消去で係数 zs[j, i] を 0 にするときは、zs の j 行目と i 行目の排他的論理和を計算するだけです。zs[i, i] が 0 で zs[i:, i:] がすべて 0 のときは複数の解がある場合です。i の値が階数になるので、それを m にセットして break でループを脱出します。zs[i:, i:] の要素に 1 がある場合は解くことができません。return で None を返します。

後退代入をする前に、zs の n 列目の値を temp に保存しておきます。m が n よりも小さい場合、2 ** (n - m) 個の解が存在します。最初の for ループで 0 から 2 ** (n - m) - 1 の値を生成し、それを使って zs[m:n, n] に 0 または 1 をセットします。後退代入も簡単で、内積の結果を mod 2 にして、それと zs[k, n] の排他的論理和をとるだけです。yield で解を出力した後、zs[:, n] に temp を代入して元の値に戻します。

●実行結果

それでは実際に試してみましょう。日月離反 - atobirabon によると、『どんなm×nライツアウトも，すべてのライトがonのときは必ず解ける．』とのことなので、全て点灯した状態を解いてみましょう。ライツアウトの数学的な解析を公開されている作者様に感謝いたします。

リスト : 簡単なテスト

for xs in lo(np.full(9, 1, dtype=np.uint8).reshape((3, 3))):
    print(xs, xs.sum())
for xs in lo(np.full(12, 1, dtype=np.uint8).reshape((3, 4))):
    print(xs, xs.sum())
for xs in lo(np.full(16, 1, dtype=np.uint8).reshape((4, 4))):
    print(xs, xs.sum())
for xs in lo(np.full(20, 1, dtype=np.uint8).reshape((4, 5))):
    print(xs, xs.sum())
for xs in lo(np.full(25, 1, dtype=np.uint8).reshape((5, 5))):
    print(xs, xs.sum())
for xs in lo(np.full(30, 1, dtype=np.uint8).reshape((5, 6))):
    print(xs, xs.sum())
for xs in lo(np.full(36, 1, dtype=np.uint8).reshape((6, 6))):
    print(xs, xs.sum())

実行結果は次のようになりました。ご参考までに、基本変形したあとの拡大係数行列を表示しています。

[[1 1 0 1 0 0 0 0 0 1]
 [0 1 1 0 0 1 0 0 0 1]
 [0 0 1 1 1 0 0 0 0 0]
 [0 0 0 1 0 1 1 0 0 1]
 [0 0 0 0 1 0 1 1 1 1]
 [0 0 0 0 0 1 0 1 0 0]
 [0 0 0 0 0 0 1 0 1 0]
 [0 0 0 0 0 0 0 1 0 0]
 [0 0 0 0 0 0 0 0 1 1]]
[[1 0 1]
 [0 1 0]
 [1 0 1]] 5
[[1 1 0 0 1 0 0 0 0 0 0 0 1]
 [0 1 1 1 0 0 1 0 0 0 0 0 1]
 [0 0 1 0 1 1 0 0 0 0 0 0 0]
 [0 0 0 1 1 1 0 1 0 0 0 0 1]
 [0 0 0 0 1 1 0 1 0 1 0 0 1]
 [0 0 0 0 0 1 1 1 1 0 0 0 0]
 [0 0 0 0 0 0 1 1 0 1 0 1 1]
 [0 0 0 0 0 0 0 1 1 1 1 0 0]
 [0 0 0 0 0 0 0 0 1 0 1 1 1]
 [0 0 0 0 0 0 0 0 0 1 1 1 1]
 [0 0 0 0 0 0 0 0 0 0 1 1 0]
 [0 0 0 0 0 0 0 0 0 0 0 1 1]]
[[1 1 1 1]
 [1 0 0 1]
 [1 1 1 1]] 10
[[1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 1]
 [0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 1]
 [0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 1]
 [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 1 0 1 1 1 0 0 0 1]
 [0 0 0 0 0 0 0 0 0 1 1 1 1 0 1 0 1]
 [0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 1 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]]
[[1 1 1 1]
 [1 0 0 1]
 [1 1 1 1]
 [0 0 0 0]] 10
[[1 0 0 0]
 [0 0 1 1]
 [0 0 1 1]
 [1 0 0 0]] 6
[[0 0 1 0]
 [1 0 0 0]
 [0 0 0 1]
 [0 1 0 0]] 4
[[0 1 0 1]
 [0 0 1 0]
 [1 1 0 1]
 [1 1 0 0]] 8
[[0 1 0 0]
 [0 0 0 1]
 [1 0 0 0]
 [0 0 1 0]] 4
[[0 0 1 1]
 [1 0 1 1]
 [0 1 0 0]
 [1 0 1 0]] 8
[[1 0 0 1]
 [0 0 0 0]
 [0 1 1 0]
 [0 1 1 0]] 6
[[1 1 1 0]
 [1 0 1 0]
 [1 0 1 0]
 [1 1 1 0]] 10
[[0 0 0 1]
 [1 1 0 0]
 [1 1 0 0]
 [0 0 0 1]] 6
[[0 1 1 0]
 [0 1 1 0]
 [0 0 0 0]
 [1 0 0 1]] 6
[[1 1 0 0]
 [1 1 0 1]
 [0 0 1 0]
 [0 1 0 1]] 8
[[1 0 1 1]
 [0 1 1 1]
 [1 1 1 0]
 [1 1 0 1]] 12
[[1 0 1 0]
 [0 1 0 0]
 [1 0 1 1]
 [0 0 1 1]] 8
[[1 1 0 1]
 [1 1 1 0]
 [0 1 1 1]
 [1 0 1 1]] 12
[[0 1 1 1]
 [0 1 0 1]
 [0 1 0 1]
 [0 1 1 1]] 10
[[0 0 0 0]
 [1 1 1 1]
 [1 0 0 1]
 [1 1 1 1]] 10
[[1 1 0 0 0 1 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 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 1 1 1 1 0 1 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 1 1 0 1 1 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 1]
 [0 0 0 0 0 0 0 1 1 1 0 1 0 1 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 1 1 0 1 1 0 1 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 1 1 1 0 0 1 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 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 1]
 [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 1 1 1 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 1 1 1 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 1 0]]
[[0 1 1 1 0]
 [0 1 0 1 0]
 [0 1 0 1 0]
 [0 1 1 1 0]] 10
[[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]]
[[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]] 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]] 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]] 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]] 15
[[1 1 0 0 0 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 1]
 [0 1 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 0 1 0 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 0]
 [0 0 0 1 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 1]
 [0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 1 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 0 0 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 0 1 1 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [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]
 [0 0 0 0 0 0 0 0 0 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 0 0 0 0 0 0 0 1 1 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 1 0 0 0 0 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 0 0 0 0 0 0 0 0 0 1]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 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 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 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 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 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 1 0 1 1 0 0 1 0 0 0 0 0 1]
 [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0 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 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 1 1 1 1 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 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 1 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 0 0 1 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 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 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 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 0 0 0 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 0 0 0 0 0 0 0 0 1 0]]
[[0 0 1 1 0 0]
 [1 0 1 1 0 1]
 [0 1 0 0 1 0]
 [1 0 1 1 0 1]
 [0 0 1 1 0 0]] 14
[[1 1 0 ..., 0 0 1]
 [0 1 1 ..., 0 0 1]
 [0 0 1 ..., 0 0 0]
 ...,
 [0 0 0 ..., 1 1 0]
 [0 0 0 ..., 1 0 0]
 [0 0 0 ..., 0 1 1]]
[[1 0 1 1 0 1]
 [0 1 1 1 1 0]
 [1 1 1 1 1 1]
 [1 1 1 1 1 1]
 [0 1 1 1 1 0]
 [1 0 1 1 0 1]] 28

4 * 4 盤は 16 通り、5 * 5 盤は 4 通りの解があり、3 * 3, 3 * 4, 4 * 5, 6 * 6 の各盤はユニーク解でした。実行速度ですが、試しに 25 * 25 盤の全点灯パターンを解いたところ、解はユニークで手数が 353 手、時間は約 0.3 秒 (Ubuntu 22.04 LTS, WSL2 + WSLg, Windows 10, Intel Core i5-6200U 2.30GHz) でした。Ｃ/Ｃ++ などで書き直すともっと速くなると思いますが、大きな盤面を解くのでなければ、このままで十分なように思います。素の Python ではもっと遅くなるかもしれませんね。

●プログラムリスト

#
# lo.py : ライツアウトの解法
#
#         Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np

# 拡大係数行列の生成
def make_matrix(ys):
    n, m = ys.shape
    k = n * m
    xs = np.zeros((k, k + 1), dtype=np.uint8)
    for y in range(n):
        for x in range(m):
            z = y * m + x
            for dy, dx in [(-1, 0), (0, -1), (0, 0), (0, 1), (1, 0)]:
                y1 = y + dy
                x1 = x + dx
                if 0 <= x1 < m and 0 <= y1 < n:
                    xs[z, y1 * m + x1] = 1
    xs[:, k] = ys.flatten()
    return xs

# ピボット選択
def select_pivot(xs, i):
    k = xs[i:,i].argmax() + i  # 0 と 1 しかないので、絶対値を求める必要はない
    if k != i:
        temp = xs[i].copy()
        xs[i] = xs[k]
        xs[k] = temp

# ライツアウトの解法
def lo(ys):
    zs = make_matrix(ys)
    n = len(zs)
    m = n
    # 前進消去
    for i in range(n):  # 解の有無をチェックするため、最後の行まで調べる
        select_pivot(zs, i)
        if zs[i, i] == 0:
            if np.all(zs[i:, i:] == 0):
                # 複数の解がある
                m = i
                break
            return    # 解無し
        for j in range(i + 1, n):
            if zs[j, i] == 1: zs[j, i:] ^= zs[i, i:]

    temp = zs[:, n].copy()
    for j in range(0, 2 ** (n - m)):
        for k in range(m, n): zs[k, n] = (j >> (k - m)) & 1
        # 後退代入
        # zs[m:, n] の値は確定済み
        for k in range(m - 1, -1, -1):
            zs[k, n] ^= (zs[k, k+1:n] @ zs[k+1:, n]) % 2
        yield zs[:, n].reshape(ys.shape)
        zs[:, n] = temp

●８めくりパズルの解法

「８めくり」はライツアウトに類似のパズルです。ルールは簡単で、あるボタンを押すと周囲のボタンの状態が反転します。つまり、光っているボタンは消灯し、消えていたボタンは点灯します。次の図を見てください。

中央のボタン 4 を押すと、その周囲のボタン 8 個の状態が反転します。押したボタンの状態は反転しません。もう一度同じボタンを押すと、再度ボタンの状態が反転するので、元の状態に戻ります。隅のボタン 0 を押すと 3 個のボタンの状態が反転し、辺にあるボタン 1 を押すと 5 個のボタンの状態が反転します。

８めくりの解法プログラムも簡単に作ることができます。プログラムと実行結果を示します。

#
# turn8.py : ８めくりパズルの解法
#
#            Copyright (C) 2018-2023 Makoto Hiroi
#
import numpy as np

# 拡大係数行列の生成
def make_matrix(ys):
    n, m = ys.shape
    k = n * m
    xs = np.zeros((k, k + 1), dtype=np.uint8)
    for y in range(n):
        for x in range(m):
            z = y * m + x
            for dy, dx in [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]:
                y1 = y + dy
                x1 = x + dx
                if 0 <= x1 < m and 0 <= y1 < n:
                    xs[z, y1 * m + x1] = 1
    xs[:, k] = ys.flatten()
    return xs

# ピボット選択
def select_pivot(xs, i):
    k = xs[i:,i].argmax() + i    # 0 と 1 しかないので、絶対値を求める必要はない
    if k != i:
        temp = xs[i].copy()
        xs[i] = xs[k]
        xs[k] = temp

def turn8(ys):
    zs = make_matrix(ys)
    n = len(zs)
    m = n
    # 前進消去
    for i in range(n):
        select_pivot(zs, i)
        if zs[i, i] == 0:
            if np.all(zs[i:, i:] == 0):
                # 複数の解がある
                m = i
                break
            return    # 解無し
        for j in range(i + 1, n):
            if zs[j, i] == 1: zs[j, i:] ^= zs[i, i:]

    temp = zs[:, n].copy()
    for j in range(0, 2 ** (n - m)):
        for k in range(m, n): zs[k, n] = (j >> (k - m)) & 1
        # 後退代入
        for k in range(m - 1, -1, -1):
            zs[k, n] ^= (zs[k, k+1:n] @ zs[k+1:, n]) % 2
        yield zs[:, n].reshape(ys.shape)
        zs[:, n] = temp

# 簡単なテスト
for x in range(4, 11):
    print("-----", x, x, "-----")
    for xs in turn8(np.full(x * x, 1, dtype=np.uint8).reshape((x, x))):
        print(xs, xs.sum())

----- 4 4 -----
[[1 0 0 1]
 [0 1 1 0]
 [1 0 0 1]
 [0 0 0 0]] 6
[[0 1 1 1]
 [1 1 0 0]
 [1 0 1 0]
 [1 0 0 0]] 8
[[0 0 1 0]
 [0 1 1 1]
 [1 1 1 0]
 [0 1 0 0]] 8
[[1 1 0 0]
 [1 1 0 1]
 [1 1 0 1]
 [1 1 0 0]] 10
[[0 1 0 0]
 [1 1 1 0]
 [0 1 1 1]
 [0 0 1 0]] 8
[[1 0 1 0]
 [0 1 0 0]
 [0 1 0 0]
 [1 0 1 0]] 6
[[1 1 1 1]
 [1 1 1 1]
 [0 0 0 0]
 [0 1 1 0]] 10
[[0 0 0 1]
 [0 1 0 1]
 [0 0 1 1]
 [1 1 1 0]] 8
[[1 1 1 0]
 [0 0 1 1]
 [0 1 0 1]
 [0 0 0 1]] 8
[[0 0 0 0]
 [1 0 0 1]
 [0 1 1 0]
 [1 0 0 1]] 6
[[0 1 0 1]
 [0 0 1 0]
 [0 0 1 0]
 [0 1 0 1]] 6
[[1 0 1 1]
 [1 0 0 0]
 [0 0 0 1]
 [1 1 0 1]] 8
[[0 0 1 1]
 [1 0 1 1]
 [1 0 1 1]
 [0 0 1 1]] 10
[[1 1 0 1]
 [0 0 0 1]
 [1 0 0 0]
 [1 0 1 1]] 8
[[1 0 0 0]
 [1 0 1 0]
 [1 1 0 0]
 [0 1 1 1]] 8
[[0 1 1 0]
 [0 0 0 0]
 [1 1 1 1]
 [1 1 1 1]] 10
----- 5 5 -----
----- 6 6 -----
[[1 0 0 0 0 1]
 [0 1 1 1 1 0]
 [0 1 0 0 1 0]
 [0 1 0 0 1 0]
 [0 1 1 1 1 0]
 [1 0 0 0 0 1]] 16
----- 7 7 -----
----- 8 8 -----
[[0 1 1 0 0 1 1 0]
 [1 1 0 1 1 0 1 1]
 [1 0 1 0 0 1 0 1]
 [0 1 0 0 0 0 1 0]
 [0 1 0 0 0 0 1 0]
 [1 0 1 0 0 1 0 1]
 [1 1 0 1 1 0 1 1]
 [0 1 1 0 0 1 1 0]] 32
----- 9 9 -----
----- 10 10 -----
[[1 0 1 1 1 1 1 1 0 1]
 [0 1 0 1 0 0 1 0 1 0]
 [1 0 0 1 1 1 1 0 0 1]
 [1 1 1 0 0 0 0 1 1 1]
 [1 0 1 0 1 1 0 1 0 1]
 [1 0 1 0 1 1 0 1 0 1]
 [1 1 1 0 0 0 0 1 1 1]
 [1 0 0 1 1 1 1 0 0 1]
 [0 1 0 1 0 0 1 0 1 0]
 [1 0 1 1 1 1 1 1 0 1]] 60

ライツアウトと違って、盤のサイズによっては全点灯パターンが解けないこともあります。実行速度ですが、32 * 32 盤を解くのに 0.79 秒かかりました。なお、deepgreen さんが作成された解法プログラム (turnover3.exe) を使うと、32 * 32 盤でも 94 msec (Windows 10, Intel Core i5-6200U 2.30GHz) で解くことができます。ソースコードも公開されているので、興味のある方は deepgreen さんの Web ページをお読みください。deepgreen さんに感謝いたします。

ところで、拙作のページ Scheme Programming: パズルに挑戦８めくりの解答で最長手数を求めたことがあります。4 * 4 盤の場合、最長手数は 7 手で、局面の総数は全部のボタンが消灯した状態を含めて 4096 通りになりました。全局面の 1 / 16 しかありません。

4 * 6 盤の場合、最長手数は 24 手で、全局面数は 2 ^ 24 = 16777216 通りになります。5 * 5 盤の場合、全消灯に到達できる局面は 2 ^ 25 / 2 = 16777216 通りあり、その中で最長手数は 20 手 (126 通り) になります。興味のある方はいろいろ試してみてください。

Python3 Programming