点Nの軌跡

競プロの話または日記

南京錠マインスイーパソルバを書く

ストーリー

「このネジを外せば、この部分がとれるのかな?」

フタを外しました。

1* 2 1
2 3 1*
1* 2 1

「なんだろうね……この数字は?」

フタを外しました。

0 1 1*
1 2 1
1* 1 0

「あっ、また数字の羅列!」

「この数字が穴の位置を示している、ということかな?」

フタを外しました。

2 3 2 2 2 4 3 2 1 1 2 2 3 2 2 2 2 1 2 2 3 3 3 3 1
2 4 3 3 3 5 4 2 2 2 3 2 4 3 3 3 4 4 5 5 6 6 5 4 1
2 4 3 3 3 4 4 2 4 3 5 4 5 4 2 4 4 5 5 6 8 7 6 5 3
1 4 4 5 5 4 4 1 3 2 3 2 4 4 3 3 4 6 7 7 7 6 5 5 3
2 4 3 4 5 4 5 2 4 2 3 2 4 4 4 4 5 6 7 6 6 4 4 4 3
3 4 3 5 6 4 4 2 3 1 1 1 3 3 4 4 6 7 7 5 4 4 4 4 2
3 4 2 3 3 2 3 3 4 2 1 2 3 3 4 5 7 7 6 5 5 6 5 4 2
3 4 3 4 3 2 3 4 3 1 0 3 4 5 4 5 5 5 3 4 4 7 6 7 4
3 5 4 4 2 2 3 5 4 3 2 3 5 5 5 5 6 6 4 4 3 5 5 6 4
3 4 3 3 2 2 4 5 4 2 3 4 6 5 5 5 6 6 5 4 2 2 3 5 4
3 5 3 4 2 3 3 5 5 5 6 5 5 2 4 4 7 6 6 3 2 1 2 2 2
3 5 3 4 2 3 3 5 4 4 5 5 4 1 3 3 6 5 6 3 3 2 2 1 1
3 5 3 4 2 4 4 6 4 5 5 5 4 2 3 2 5 4 5 2 3 3 4 3 2
4 5 3 3 2 3 3 4 2 3 3 4 3 2 2 1 4 4 5 2 2 2 3 4 3
3 3 1 2 3 5 4 3 2 4 5 5 4 3 2 0 2 3 4 3 3 3 3 4 3
3 3 1 2 4 5 3 2 2 5 5 4 2 1 1 0 2 3 3 3 4 4 2 2 2
2 3 1 2 4 6 5 3 3 6 7 5 4 3 3 1 2 4 5 5 5 5 4 4 3
3 4 2 2 3 5 4 3 3 6 6 5 3 3 3 3 5 6 6 5 5 6 6 7 5
2 3 3 3 3 4 4 4 4 5 5 4 3 3 3 3 5 7 8 6 5 5 6 7 5
2 3 5 5 4 3 2 4 4 4 3 4 3 3 3 5 6 6 5 4 5 6 7 7 5
3 4 6 5 4 1 1 3 3 2 1 3 4 5 4 5 5 5 4 3 4 5 6 6 4
4 6 8 6 5 2 1 2 2 3 2 4 4 6 6 8 7 5 2 2 3 5 5 6 4
4 5 6 5 6 3 2 2 2 4 2 4 3 6 5 6 6 6 4 4 4 5 4 5 4
2 3 4 5 7 5 3 2 2 5 3 4 1 3 3 5 5 5 4 5 5 4 2 3 3
1 1 1 3 5 4 2 1 1 3 2 3 1 2 1 2 2 3 3 4 4 3 2 2 2

「これで3つ目だけど……やっぱり内装は同じなのね」

「同じ……?」

「鍵が開いてた錠を元にして、法則を解くしかなさそうよ」

「これじゃ法則がわかっても不可能だろ……」

南京錠マインスイーパソルバを書く。法則はもう知っているので、ざっくり、次のような問題が解きたいということになる。

問題文1

HW列の二次元配列Aが与えられます。

以下の条件を満たすH行の文字列S_{i} (1\le i\le H)を出力してください。

すべてのテストケースについて、解はただひとつ存在することが保証されます。

条件

  •  S_{i}.#からなるW文字の文字列である

  •  S_{ij} (1\le i\le H, 1\le j\le W)およびその8近傍に#はちょうどA_{ij}個存在する

制約

1\le H,W\le 25

入力例

3 3
2 3 2
2 4 3
1 2 2

出力例

#.#
.#.
..#

解法を考える

bit全探索はO(HW2^{HW})で無理。深さ優先探索も厳しそう。

よく考えると、連立一次方程式が見えてくる。

簡単のために3\times 3の場合を考える。A_{ij}は結局S_{ij}とその8近傍の#の和なので、

a b c
d e f
g h i

と記号をふってやると、

a+b+d+e=A_{11}

a+b+c+d+e+f=A_{12}

...

e+f+h+i=A_{33}

みたいな感じで式が9個できるので、これを解けばいい。一般に式の数と変数の数はHW個。この連立方程式 O((HW)^{3})で解けるので、HWが25ぐらいなら現実的な時間で動く。係数行列を正しく構築してやれば、あとはPythonかなんかにやらせればいい。

問題文2

正整数H,Wが与えられます。上の問題1において、HW列のAのうち矛盾のないもの全てに対して解が一意に定まるならYes、そうでないならNoと出力してください。

制約

1\le H,W \le 10^{18}

解が一意に定まる条件

係数行列A(入力のAではない)とS_{ij}#にするかの\{0,1\}ベクトルx、あと入力を1列に並べ直したベクトルbでもってAx=bとして、x=A^{-1}bで解が求まる(最初はまじめにLU分解でも書こうかと思ったけどめんどくてPythonにやらせることにしたので細かい実装は省略)予定なのだが、気になるのは「どんな盤面サイズに対しても係数行列は逆行列をもつか、解は一意に定まるか」ということである(上の問題文2は要はこれの判定問題である)。ソルバは書いただけではだめでテストしなきゃならないので、こういうことを考えておこうという趣旨(もちろん行列式を求めて確認すればいいが、その辺の理論をちょっとぐらい考えたほうがいいだろう)。

たとえば2\times 2の盤面ではどうだろうか。

2 2
1 1
1 1

1か所だけ#にすればどこでも同じになるので、答えは一意に定まらない。ほかにも

2 5
2 2 1 0 0  
2 2 1 0 0

に対する答えは以下のどっちでも成立してしまう。

#....  .#...
.#...  #....

はっきりとした答えがだせなかったので小さめの係数行列を構築して行列式を求め、それが非零であるかをチェックしてみる。

import numpy as np

table = [[0] * 10 for i in range(10)]

for h in range(1,11):
    for w in range(1,11):
        A = [[0]*(w*h) for i in range(w*h)];
        for a in range(w*h):
            ay = a//w
            ax = a%w
            for b in range(w*h):
                by = b//w
                bx = b%w
                if abs(ay-by) <= 1 and abs(ax-bx) <= 1:
                    A[a][b] += 1
        A = np.array(A)
        d = np.linalg.det(A)
        table[h-1][w-1] = d
        print('%2d x %2d det A = %d' % (h,w,d))

print('h/w ', end='')
for i in range(1,11):
    print('%2d ' % i, end='')
print('')
for h in range(10):
    print(' %2d ' % (h+1), end='')
    for w in range(10):
        print('%2d ' % table[h][w], end='')
    print('')

必要な部分だけ出力を抽出するとこんな感じ。行列式\pm 10で、H\equiv 2(mod\ 3)またはW\equiv 2(mod\ 3)のいずれかが成立していれば0になるっぽい。そうじゃない場合は(ある程度パターンは見えるが)いまいち規則性が分からない。

h/w  1  2  3  4  5  6  7  8  9 10
  1  1  0 -1 -1  0  1  1  0 -1 -1
  2  0  0  0  0  0  0  0  0  0  0
  3 -1  0  1 -1  0  1 -1  0  1 -1
  4 -1  0 -1  1  0  1 -1  0 -1  1
  5  0  0  0  0  0  0  0  0  0  0
  6  1  0  1  1  0  1  1  0  1  1
  7  1  0 -1 -1  0  1  1  0 -1 -1
  8  0  0  0  0  0  0  0  0  0  0
  9 -1  0  1 -1  0  1 -1  0  1 -1
 10 -1  0 -1  1  0  1 -1  0 -1  1

とりあえずH\equiv 2(mod\ 3)W\equiv 2(mod\ 3)の場合を処理したい。どっちでも一緒なので紙面の都合を考えてWが該当した場合で進める。 A_{1}=\{1,1,...,1\},A_{2}=\{1,1,...,1\}の場合、

  • S_{1}=#.[.#.][.#.]...[.#.]S_{2},S_{3},...=... ... ...

  • S_{1}=[.#.][.#.]...[.#.].#S_{2},S_{3},...=... ... ...

の2パターンが考えられる。[]の中が3文字、外に2文字で、[]の数を変えることで、3で割った余りが2であるようないかなるWに対しても解が一意に定まらないAが存在するとわかる。

ということで、残る問題は「H mod\ 3 \neq 2かつW mod\ 3\neq 2なら係数行列の行列式は必ず非零か」となる。

実際に係数行列がどうなっているのか見ないことには始まらないので見てみる。

import numpy as np

h,w = map(int,input().split())

A = [[0]*(w*h) for i in range(w*h)];
for a in range(w*h):
    ay = a//w
    ax = a%w
    for b in range(w*h):
        by = b//w
        bx = b%w
        if abs(ay-by) <= 1 and abs(ax-bx) <= 1:
            A[a][b] += 1
for y in range(h*w):
    print(A[y])

A = np.array(A)
d = np.linalg.det(A)
print('det A = ',d)

H=3,W=4のときの係数行列はこんな感じ。

[1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0]
[1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0]
[0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0]
[0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0]
[1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0]
[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0]
[0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1]
[0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1]
[0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0]
[0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0]
[0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1]

H=4,W=3はこんな感じ。

[1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0]
[1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
[0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0]
[1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]
[0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0]
[0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0]
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1]

目が疲れそうだが、要は「サイズWの三重対角行列の、サイズHの三重対角行列(特に非零成分はすべて1)」のような感じになっている。主対角線とその上下の成分がすべて1でサイズnの三重対角行列をT_{n}と書くことにすれば、

H=3,W=4の係数行列は

 \displaystyle
\begin{pmatrix}
T_{4}&T_{4}&0 \\
T_{4}&T_{4}&T_{4} \\
0&T_{4}&T_{4} \\
\end{pmatrix}

H=4,W=3の係数行列は

 \displaystyle
\begin{pmatrix}
T_{3}&T_{3}&0&0 \\
T_{3}&T_{3}&T_{3}&0 \\
0&T_{3}&T_{3}&T_{3} \\
0&0&T_{3}&T_{3} \\
\end{pmatrix}

といった感じである。

さて、三重対角行列の行列式は簡単な漸化式で求めることができる。英語版Wikipediaにまとめられているので参考にしつつ、今回は特殊なケースゆえもっと簡単になって、結局|T_{n}|=|T_{n-1}|-|T_{n-2}|, (|T_{0}|=1,|T_{-1}|=0)となる。|T_{n}|は周期6で循環し、以下のようになる。

n 1 2 3 4 5 6
|Tn| 1 0 -1 -1 0 1

で、今問題なのはこれの三重対角行列の行列式が欲しいというところだが、こちらもほぼ同様の性質が現れる。

T_{n}の三重対角行列U_{m}行列式は、実際に計算してみると

|U_{1}|=|T_{n}|

|U_{2}|=0

|U_{3}|=(-1)^{n}|T_{n}|^{3}

|U_{4}|=(-1)^{n}|T_{n}|^{4}

|U_{5}|=0

|U_{6}|=|T_{n}|^{6}

以下想像の通りである(と願いたい)。結局先ほどプログラムを書いて確認した行列式の値はこれに基づいていたということになる。さらに、ここからH mod\ 3 \neq 2かつW mod\ 3 \neq 2なら係数行列に逆行列が存在することがわか(った気にな)る。証明でも何でもない。

ソルバ実装

あとは心ゆくままにソルバと問題ジェネレータを作るだけ。というわけで完成。三重対角行列の性質がいい感じにしてくれるらしく逆行列の成分もどうやら-1か0か1のようで、精度に注意する必要はまずないはず(???)。あまり変な入力をされても壊れちゃうだけなので乱暴はしないでもらえるとありがたい。効率面での改善の余地に関しては、正直その道の人が見ればオーダーがひとつふたつ下がるとか、それぐらいはありそうな気もする。逆行列はどうも高速に計算できそう。

ソルバ。最大25\times 25まで対応できる。それなりな環境で実行すれば25\times 25でも1秒いらないぐらいだと思う。解が複数ある場合にその一例を示す、はさすがに諦めた。

import numpy as np
import sys

print('')
print('input problem size (if HxW -> H W)')
print('')
print('---matrix size---')
H,W = map(int,input().split())

if H%3==2 or W%3==2:
    print('cannot calc')
    sys.exit()

N = H*W

if N>25*25:
    print('too big!')
    sys.exit()

print('')
print('input problem matrix like this')
print('a b c')
print('d e f')
print('g h i')
print('')
print('---matrix---')
B = [list(map(int,input().split())) for i in range(H)]

vb = []
for y in range(H):
    for x in range(W):
        vb.append(B[y][x])
vb = np.array(vb).reshape(N,1)

A = [[0]*N for i in range(N)];
for a in range(N):
    ay = a//W
    ax = a%W
    for b in range(N):
        by = b//W
        bx = b%W
        if abs(ay-by) <= 1 and abs(ax-bx) <= 1:
            A[a][b] += 1

d = np.linalg.det(A)
if -0.5 < d < 0.5:
    print('cannot calc')
    sys.exit()

Ainv = np.linalg.inv(A)

X = np.matmul(Ainv,vb).reshape(H,W)
ans = [''] * H
Bcheck = [[0] * W for i in range(H)]
for y in range(H):
    for x in range(W):
        if X[y][x] > 0.5:
            ans[y] += '#'
            for dy in range(-1,2):
                for dx in range(-1,2):
                    if 0 <= y+dy < H and 0 <= x+dx < W:
                        Bcheck[y+dy][x+dx] += 1
        else:
            ans[y] += '.'

if B!=Bcheck:
    print('solution not found')
    sys.exit()

print('')
print('---answer---')
for y in range(H):
    print(ans[y])

問題ジェネレータ。ネジ存在確率指定、ネジ個数指定、手動ネジ座標設定の3つのモードで簡単問題生成。ほんとは機械的に生成した盤面に人の手を加えられたりするといい気がするけど、それはやりたい人がやって。

import numpy as np
import sys

print('')
print('input problem size (if HxW -> H W)')
print('')
print('---matrix size---')
H,W = map(int,input().split())

if H%3==2 or W%3==2:
    print('cannot calc')
    sys.exit()

N = H*W

if N>25*25:
    print('too big!')
    sys.exit()

print('')
print('select mode')
print('1:probability')
print('2:mine number')
print('3:handmade')

mode = int(input())

A = [[0]*W for i in range(H)]
used  = [[0]*W for i in range(H)]

print('')
if mode==1:
    print('input probability, range:[0.0 ~ 1.0]')
    print('probability:',end='')
    p = float(input())
    arr = np.random.rand(H,W)
    for y in range(H):
        for x in range(W):
            if arr[y][x] < p:
                used[y][x] = 1
                for dy in range(-1,2):
                    for dx in range(-1,2):
                        if 0 <= y+dy < H and 0 <= x+dx < W:
                            A[y+dy][x+dx] += 1
elif mode==2:
    print('input number of mines, range:','[%d ~ %d]' % (0,H*W))
    print('number:',end='')
    n = int(input())
    n = min(H*W,max(0,n))
    while n>0:
        y = np.random.randint(H)
        x = np.random.randint(W)
        if used[y][x]==0:
            used[y][x] = 1
            for dy in range(-1,2):
                for dx in range(-1,2):
                    if 0 <= y+dy < H and 0 <= x+dx < W:
                        A[y+dy][x+dx] += 1
            n -= 1
elif mode==3:
    print('input position, format:y x range:y...[0 ~ %d], x...[0 ~ %d]' % (H-1,W-1))
    print('you can set mines as many as you like by entering multiple lines')
    print('to finish setting, input \'-1 -1\'')
    print('---your input---')
    while True:
        y,x = map(int,input().split())
        if y<0 and x<0:
            break
        y = min(y,H-1)
        x = min(x,W-1)
        if used[y][x]==0:
            used[y][x] = 1
            for dy in range(-1,2):
                for dx in range(-1,2):
                    if 0 <= y+dy < H and 0 <= x+dx < W:
                        A[y+dy][x+dx] += 1
else:
    print('bad input')
    sys.exit()

print('---mines---')
ans = [''] * H
for y in range(H):
    for x in range(W):
        if used[y][x]==1:
            ans[y] += '#'
        else:
            ans[y] += '.'

for y in range(H):
    print(ans[y])

print('---problem---')
print(H,W)
for y in range(H):
    print(*A[y])

問題例(冒頭の問題)

25 25
2 3 2 2 2 4 3 2 1 1 2 2 3 2 2 2 2 1 2 2 3 3 3 3 1
2 4 3 3 3 5 4 2 2 2 3 2 4 3 3 3 4 4 5 5 6 6 5 4 1
2 4 3 3 3 4 4 2 4 3 5 4 5 4 2 4 4 5 5 6 8 7 6 5 3
1 4 4 5 5 4 4 1 3 2 3 2 4 4 3 3 4 6 7 7 7 6 5 5 3
2 4 3 4 5 4 5 2 4 2 3 2 4 4 4 4 5 6 7 6 6 4 4 4 3
3 4 3 5 6 4 4 2 3 1 1 1 3 3 4 4 6 7 7 5 4 4 4 4 2
3 4 2 3 3 2 3 3 4 2 1 2 3 3 4 5 7 7 6 5 5 6 5 4 2
3 4 3 4 3 2 3 4 3 1 0 3 4 5 4 5 5 5 3 4 4 7 6 7 4
3 5 4 4 2 2 3 5 4 3 2 3 5 5 5 5 6 6 4 4 3 5 5 6 4
3 4 3 3 2 2 4 5 4 2 3 4 6 5 5 5 6 6 5 4 2 2 3 5 4
3 5 3 4 2 3 3 5 5 5 6 5 5 2 4 4 7 6 6 3 2 1 2 2 2
3 5 3 4 2 3 3 5 4 4 5 5 4 1 3 3 6 5 6 3 3 2 2 1 1
3 5 3 4 2 4 4 6 4 5 5 5 4 2 3 2 5 4 5 2 3 3 4 3 2
4 5 3 3 2 3 3 4 2 3 3 4 3 2 2 1 4 4 5 2 2 2 3 4 3
3 3 1 2 3 5 4 3 2 4 5 5 4 3 2 0 2 3 4 3 3 3 3 4 3
3 3 1 2 4 5 3 2 2 5 5 4 2 1 1 0 2 3 3 3 4 4 2 2 2
2 3 1 2 4 6 5 3 3 6 7 5 4 3 3 1 2 4 5 5 5 5 4 4 3
3 4 2 2 3 5 4 3 3 6 6 5 3 3 3 3 5 6 6 5 5 6 6 7 5
2 3 3 3 3 4 4 4 4 5 5 4 3 3 3 3 5 7 8 6 5 5 6 7 5
2 3 5 5 4 3 2 4 4 4 3 4 3 3 3 5 6 6 5 4 5 6 7 7 5
3 4 6 5 4 1 1 3 3 2 1 3 4 5 4 5 5 5 4 3 4 5 6 6 4
4 6 8 6 5 2 1 2 2 3 2 4 4 6 6 8 7 5 2 2 3 5 5 6 4
4 5 6 5 6 3 2 2 2 4 2 4 3 6 5 6 6 6 4 4 4 5 4 5 4
2 3 4 5 7 5 3 2 2 5 3 4 1 3 3 5 5 5 4 5 5 4 2 3 3
1 1 1 3 5 4 2 1 1 3 2 3 1 2 1 2 2 3 3 4 4 3 2 2 2

ソルバの解答

..#..##......#.#...#..#..
##..#.#..#.##...#..#.##.#
..#..#...#...#..#######..
..#..#.#.#.##.#.#..###.##
#.####.......#...###..#.#
#....#.#.#...#.#####.#...
#..##..#....#..##.#..##.#
#.#.....#...#..#.##.###.#
#..#..##....###.#...#.###
#.##..#.##.#.#.####....#.
#......#...##..#.###....#
#.#.#.#.####...#.#...#...
#.##..#.#..#...#.##..#...
#.....##..#.##...#.#..###
##..#.....#......#......#
....###..###.#....#.##...
#....#..#.#......#..##..#
#.#..##..###.##..###..###
.#..#.#..##.#..####.#####
...#...##........###.#..#
.####...#..##.###....####
##.#........###.##..#.##.
.###.#..#.#...####..#...#
#..###..#.#.#.#..##.##.##
....##....#.....#..##....