Lights OutをSagemathで解く（その３）

[前回の投稿（その２）](https://qiita.com/masa0599/items/6995363cce5f45888f49)ではLight Outゲーム(5 x 5)をSagemathを使って解きましたが。今回は一般的に（m x n）のLight Outが必ず解けるかを判定する方法を考えます。

### （m x n）のLight Outが必ず解けるか

既に[前々回](https://qiita.com/masa0599/items/0302ef923eaebf96da64)でボタン・ライトの反転規則の行列を作る関数 **makeA(m,x)** を作ってありますので。この行列式が$0$でなければ必ず解けて、解は１つしかないことが分かります。

しかし5 x 5の場合のように行列式が$0$の時には、解けるかを判定したり、解の数を知るために直交基底ベクトルの数= **自由度（行の数ー階数）** を求めた方が役に立ちます。

Sagemathでは階数は **rank()** プロパティを用いてば簡単に求まるので、$1 \le m, n \le 15$で求めてみます。

```python 
m, n = 15,15
print(np.array([[i*j-makeA(i, j).rank() for i in range(1,m+1)] for j in range(1,n+1)]))
# [[0 1 0 0 1 0 0 1 0 0 1 0 0 1 0]
# [1 0 2 0 1 0 2 0 1 0 2 0 1 0 2]
# [0 2 0 0 3 0 0 2 0 0 3 0 0 2 0]
# [0 0 0 4 0 0 0 0 4 0 0 0 0 4 0]
# [1 1 3 0 2 0 4 1 1 0 4 0 1 1 4]
# [0 0 0 0 0 0 0 6 0 0 0 0 0 0 0]
# [0 2 0 0 4 0 0 2 0 0 7 0 0 2 0]
# [1 0 2 0 1 6 2 0 1 0 2 0 7 0 2]
# [0 1 0 4 1 0 0 1 8 0 1 0 0 5 0]
# [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
# [1 2 3 0 4 0 7 2 1 0 6 0 1 2 8]
# [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
# [0 1 0 0 1 0 0 7 0 0 1 0 0 1 0]
# [1 0 2 4 1 0 2 0 5 0 2 0 1 4 2]
# [0 2 0 0 4 0 0 2 0 0 8 0 0 2 0]]
```

このアルゴリズムは$1 \le m, n \le 15$の範囲では問題ないのですが、これがもう少し大きくなると難しくなってきます。効率の良いアルゴリズムが「[群論の味わい －置換群で解き明かすルービックキューブと15パズル－](https://www.amazon.co.jp/dp/4320019415/)」の中にありましたので紹介します。

> **定理 6.5.3** パズルは$GCD(f_{M+1}(x)$,$f_{N+1}(x+1) )=1$の時にどんな初期状態からも解くことが出来る。ここで$f_{n}(x)$はフィボナッチ多項式

まず[フィボナッチ多項式(Wikipedia)](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81%E5%A4%9A%E9%A0%85%E5%BC%8F)ですがこの説明にあるように。フィボナッチ数列を多項式にしたものということです。

```math 
F_n(x)= \left \{
\begin{array}{ll}
0, &if \ n = 0 \\
1, &if \ n = 1 \\
xF_{n-1}(x)+F_{n-2}(x) & if \ n \ge 2
\end{array}
\right.
```

これをSagemathでの実装すると以下のようになります。**剰余環$\mathbb{Z}/2 \mathbb{Z}$** で演算する必要があるので**GF(2)** で定義しています。

プログラムではn = 9までのフィボナッチ多項式を求めていますが、後のGCDの説明のために**factor()** を使って因数分解された形で表示しています。

```python
R = PolynomialRing(GF(2),"x")
x = R.gen()
cache = dict()
def f(t,n):
    if n<2: return R(n)
    else:
        if (t,n) not in cache:
            cache[(t,n)] = f(t,n-2)+t*f(t,n-1)
        return cache[(t,n)]
    
for n in range(1,10):
    print(f"n={n}, f(x,{n})={f(x,n).factor()}")
# n=1, f(x,1)=1
# n=2, f(x,2)=x
# n=3, f(x,3)=(x + 1)^2
# n=4, f(x,4)=x^3
# n=5, f(x,5)=(x^2 + x + 1)^2
# n=6, f(x,6)=x * (x + 1)^4
# n=7, f(x,7)=(x^3 + x^2 + 1)^2
# n=8, f(x,8)=x^7
# n=9, f(x,9)=(x + 1)^2 * (x^3 + x + 1)^2
```

このフィボナッチ多項式を用いて5 x 5のLights Outが必ず解けるかをチェックしてみます。

```math 
f_{5+1}(x) = f_6(x) = x(x+1)^4 \\
f_{5+1}(x+1) = f_6(x+1) = (x+1)x^4 \\  
(\mathbb{Z}/2 \mathbb{Z}なのでx+2=xとなる)\\
したがって　\\
GCD(f_6(x),f_6(x+1))=(x+1)x=x^2+x
```
となり１ではないので必ずは解けないことが分かります。さらにGCDの多項式の**次数2**が行列式の**自由度** に等しくなります。

従ってm x nのLights Outの**自由度** を求める関数**F(m, n)** は以下のように書けます。

```python 
def F(m,n):
    return f(x,m+1).gcd(f(x+1,n+1)).degree()
print(F1(5,5))
# 2
```

この関数**F(m, n)** を使うと　$1 \le m, n \le 500$の範囲の**自由度** を求めることが出来ました。

```python 
m = n = 500
print(np.array([[F(i,j) for i in range(1,m+1)] for j in range(1,n+1)]))
#[[ 0  1  0 ...  0  0  1]
# [ 1  0  2 ...  0  2  0]
# [ 0  2  0 ...  0  0  2]
# ...
# [ 0  0  0 ...  0  0  0]
# [ 0  2  0 ...  0 16  2]
# [ 1  0  2 ...  0  2  0]]
```
（開発環境：Cocalc Sage Worksheet）