数独を全探索で解くプログラムとその速度向上

# はじめに
数独は楽しいパズルなんですが、難しい問題だと全然解けなくてつらい思いをするので、コンピュータの力を借りて解くことにしました。

最初は解くのに17分以上かかっていたのが、アルゴリズム改善によって最終的に3秒で解くことができた(もちろん環境は同じ)のでその過程をアウトプットします。

使用環境：Google Colaboratory (Python 3.7.10)

数独のルールは以下の通りです。9×9のマトリックスに数字を入れるパズルです。
>
![スクリーンショット 2021-05-24 20.46.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/acd27034-9197-bdcb-e664-7cb5f266f8a7.png)
・まだ数字の入っていないマスに１から９までの数字のどれかを１つずつ入れましょう。０(ゼロ)は使いません。
・タテ列（９列あります)､ヨコ列（９列あります)､太線で囲まれた３×３のブロック（９つあります）のどれにも、１から９までの数字が１つずつ入るようにします。
画像、文章は[ニコリHP](https://www.nikoli.co.jp/ja/puzzles/sudoku/)より引用


# 基本的な解き方(アルゴリズム)
全探索を行います。
具体的にいうと、左上から順番に入れることのできる可能なだけ小さい数字を入れていき、矛盾が生じたら戻ってもう少し大きい数字を入れていく方法です。
図で説明した方がわかりやすいと思います。以下の数独を解くことにします。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/4d752d41-323d-0ced-7082-615d07d15d0a.png" width="40%">
引用元：[本当に解ける人いるの？　フィンランド人数学者が作った “世界一難しい数独” が発表される](https://rocketnews24.com/2012/07/03/22654/)
【注意】この記事の下の方にこの数独の答えが載っています。自力で解きたい方は解いてから読んでください。世界一難しいらしいので大変だと思いますが。

まず、8の右隣のセルに数字を入れます。ここに入れられる最小の数字は1なので1を入れます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/a993a20b-f08a-c618-4689-41cce5074c23.png" width="40%">
同様にして右隣に数字を入れられる最小の数字を入れていきます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/bbf98ce7-4170-e201-230b-2a275292d3b2.png" width="40%">
さて、3の右隣に入れる数字ですが、1~3は同じ行の左の方、4は同じ列の下の方で使用済みなので4の次の数である5を入れます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/c3841486-f69c-4e67-db72-6328a3557878.png" width="40%">
同様にしてどんどん数字を入れていきます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/1a0baddd-6ad4-8694-ef42-c064b8c82876.png" width="40%">
さて、2行目の4の右隣には1~9のどの数字も入れることができません。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/aac52a11-9652-2690-6e4d-81dd0bb8e7d5.png" width="40%">
このような状況になったら、一つ前の数字を変えます。今回のケースでは4の次の数字である5を入れればOKです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/f6259b8e-78d6-87b6-f658-891a69dc3a44.png" width="40%">
どんどん入れていきます。行き詰まったら戻ってやり直します。その繰り返しです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/63432acd-65b5-c71e-eb0e-5e02a93562f3.png" width="40%">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/2b035ea3-37d8-867b-37b5-fd1c06809cdf.png" width="40%">

これをひたすら繰り返せば解くことができます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/f0060730-3925-866e-db52-34eb8521f7fb.png" width="40%">

# 数独を解くプログラム(基本型)
便宜上、各マス目を以下の通り番号つけます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/5b44ff63-81fb-0356-93a9-d4e4b9af3b20.png" width="40%">
基本的な解き方で見たように、以下のような手順でマス目を埋めていきます。
>① 0番目のマスに数字が入っていなければ1を入れる。
② ①の状態で「正しい状態(同じ行、列、3×3ブロックに重複する数字がない)」かチェックする。
③ 「正しい状態」でなければ2,3,4,・・・の順に入れて試してみる。「正しい状態」になったら次のマスに進む。
④ 1番目のマスに1,2,3,・・・の順に「正しい状態」になるまで試す。「正しい状態」になったら次のマスに進む。これを全マスが埋まるまで繰り返す。
⑤ もしk番目のマスに数字が既に入っていたら、そのマスはスルーして次のマスに進む。
⑥ もしk番目のマスに1~9の全ての数字を入れても「正しい状態」にならない場合は、k-1番目のマスの数字を1つ大きくして試す。それで「正しい状態」にならなければさらに1つ大きくして「正しい状態」になるか試す。全ての数字を入れてもダメだったら、さらに前のマスの数字を変える。

これにより数独を解くことができます。
再帰関数を使うとラクチンです。
なお、問題の入力はマス目の順番通りに数字を連ねた(空欄は0にする)文字列としています。
~~最初はリストでやろうとしたけどミュータブルだからなのかごちゃごちゃになった~~


```python:solve_sudoku1.py
sudoku = "800000000003600000070090200050007000000045700000100030001000068008500010090000400" #https://rocketnews24.com/2012/07/03/22654/

import numpy as np

def ck(s): #入力された数独が「正しい」状態か判定します
  for i in range(9):
    ss = [s[i*9+j] for j in range(9)] #横方向にチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
  for i in range(9):
    ss = [s[i%9+9*j] for j in range(9)] #縦方向にチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
  for i in range(9):
    ss = [s[i//3*27+i%3*3+j//3*9+j%3] for j in range(9)] #3×3ブロックでチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
  return True

def dfs(i,t): #全探索を行う部分です
  if i == 81:
    print(np.array(list(t)).reshape(9,9))
  else:
    if t[i] == "0":
      for j in range(1,10):
        t_ = t[:i]+str(j)+t[i+1:]
        if ck(t_):
          dfs(i+1,t_)
    else:
      dfs(i+1,t)

print(np.array(list(sudoku)).reshape(9,9)) #初期状態を表示
print("-----------------------------------------")
dfs(0,sudoku) #解く
```

以下のように出力されます。

```python:出力
[['8' '0' '0' '0' '0' '0' '0' '0' '0']
 ['0' '0' '3' '6' '0' '0' '0' '0' '0']
 ['0' '7' '0' '0' '9' '0' '2' '0' '0']
 ['0' '5' '0' '0' '0' '7' '0' '0' '0']
 ['0' '0' '0' '0' '4' '5' '7' '0' '0']
 ['0' '0' '0' '1' '0' '0' '0' '3' '0']
 ['0' '0' '1' '0' '0' '0' '0' '6' '8']
 ['0' '0' '8' '5' '0' '0' '0' '1' '0']
 ['0' '9' '0' '0' '0' '0' '4' '0' '0']]
-----------------------------------------
[['8' '1' '2' '7' '5' '3' '6' '4' '9']
 ['9' '4' '3' '6' '8' '2' '1' '7' '5']
 ['6' '7' '5' '4' '9' '1' '2' '8' '3']
 ['1' '5' '4' '2' '3' '7' '8' '9' '6']
 ['3' '6' '9' '8' '4' '5' '7' '2' '1']
 ['2' '8' '7' '1' '6' '9' '5' '3' '4']
 ['5' '2' '1' '9' '7' '4' '3' '6' '8']
 ['4' '3' '8' '5' '2' '6' '9' '1' '7']
 ['7' '9' '6' '3' '1' '8' '4' '5' '2']]
```

結果が出力されるまでに要した時間は26秒、探索が全て完了するまでに要した時間は17分43秒でした。
出力されるまでの時間が全探索完了に要する時間より大幅に短いのは、1,2番目の最初に数字を置くマスに入る正解の数字がそれぞれ1,2と小さい数字であることによります。もしこれらのマスに入る数字が9,8だったらおそらく結果出力までに17分くらいかかります。
この記事においては、各プログラムにおける所要時間は探索が全て完了するまでの時間(今回は17分43秒)とします。これは、問題の数字が入れ替わることにより評価が変わることを防ぐためです(この問題でも1,2と9,8が入れ替われば答えが表示されるまでの時間がかなり変わる)。

# 改善その1 「正しい状態」判定の高速化
基本型では　どこかのマス目に数字を入れたときに、「盤面全体が『正しい状態』になっているか」をチェックしています。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/eb32f25a-3834-ea9c-0eb1-faa4709995ce.png" width="50%">
↑盤面全体が「正しい状態」かチェックしている図。全ての行、列、3×3ブロックについても数値の重複がないかチェックしています(3×3ブロックについては描きにくいので省略)。

しかし、k番目にマス目に入れたときにチェックしないといけないのは、そのマス目の所属する行と列とブロックだけのはずです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/f0d8fdf7-913a-e6aa-3f92-8d0afab959e4.png" width="50%">
↑15番目のマスに数字を入れた場合にチェックしなければならない行、列、3×3ブロック

チェック対象を27箇所から3箇所に削減することで、実行時間の短縮が期待されます。

```python:solve_sudoku2.py
sudoku = "800000000003600000070090200050007000000045700000100030001000068008500010090000400" #https://rocketnews24.com/2012/07/03/22654/

import numpy as np

def ck(s,i): #入力された数独が「正しい」状態か、i番目に関してのみ判定します
    ss = [s[i//9*9+j] for j in range(9)] #横方向にチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
    ss = [s[i%9+9*j] for j in range(9)] #縦方向にチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
    ss = [s[i//27*27+i%9//3*3+j//3*9+j%3] for j in range(9)] #3×3ブロックでチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
    return True

def dfs(i,t): #全探索を行う部分です
  if i == 81:
    print(np.array(list(t)).reshape(9,9))
  else:
    if t[i] == "0":
      for j in range(1,10):
        t_ = t[:i]+str(j)+t[i+1:]
        if ck(t_,i):
          dfs(i+1,t_)
    else:
      dfs(i+1,t)

print(np.array(list(sudoku)).reshape(9,9))
print("-----------------------------------------")
dfs(0,sudoku)
```
実行時間は2分40秒(答えが出力されるまで4秒)でした。標準型と比べて実行時間は約85%削減されています。
(出力は標準型と同じなので省略。以下同様)

# 改善その2 入れられる数字の候補が少ない場所から入れていく
こちらのツイートのまんまです。

https://twitter.com/excelspeedup/status/988926717728178176?s=21

① 注目しているマスの所属する行、列、ブロックで既に使われている数字をリストアップする。そのリストに含まれない数字が対象マスに入る数字の候補となる。候補が何種類なのかを調べる。
例えば下の図のオレンジ色のマスでは{3,5,7,8,9}が既に使われているので残りは{1,2,4,6}の4種類となる。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/60b03a79-e8c3-0631-5bdc-cd51d82a6215.png" width="50%">

② ①を全てのマスに対して実行する。ただし既に数字が入っているマスは0とする。(結果は下図)
③ ②の結果のうち、最も候補数の少ないセル(赤字のセル)から数字を入れていく。最も候補数の少ないセルが複数あった場合は番号の若いセルから数字を入れていく。
④ 数字を入れるごとに、①~③を繰り返して次に入れるマスを決定する。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/0a47fa26-322d-2332-d4b0-a965a90ba889.png" width="40%">
数字を入れる順番が変わるだけ(左上から順番に入れる→候補数の少ない順番に入れる)で、全探索の進め方自体は大きく変わりません。

実装は以下の通りとなります。

```python:solve_sudoku3.py
sudoku = "800000000003600000070090200050007000000045700000100030001000068008500010090000400" #https://rocketnews24.com/2012/07/03/22654/

import numpy as np

def ck(s,i): #入力された数独が「正しい」状態か、i番目に関してのみ判定します
    ss = [s[i//9*9+j] for j in range(9)] #横方向にチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
    ss = [s[i%9+9*j] for j in range(9)] #縦方向にチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
    ss = [s[i//27*27+i%9//3*3+j//3*9+j%3] for j in range(9)] #3×3ブロックでチェック
    for j in range(1,10):
        if ss.count(str(j))>=2:
          return False
    return True

def fi(s): #次にどこのマスに数字を入れればいいか(入る候補が最も少ないのはどこか)を判定します
    _ = []
    for k in range(81):
      if s[k] != "0":
        _ += [0]
      else:
        ss_yoko = [s[k//9*9+j] for j in range(9)]
        ss_tate = [s[k%9+9*j] for j in range(9)]
        ss_block = [s[k//27*27+k%9//3*3+j//3*9+j%3] for j in range(9)]
        _ += [len(set(ss_tate+ss_yoko+ss_block))]
    k = _.index(max(_))
    return k

def dfs(i,t): #全探索を行う部分です
  if i == sudoku.count("0"):
    print(np.array(list(t)).reshape(9,9))
  else:
    k = fi(t)
    for j in range(1,10):
      t_ = t[:k]+str(j)+t[k+1:]
      if ck(t_,k):
        dfs(i+1,t_)

print(np.array(list(sudoku)).reshape(9,9))
print("-----------------------------------------")
dfs(0,sudoku)
```

実行時間は7秒(答えが出力されるまで4秒)と、改善その1と比べ96%も実行時間が短くなりました。

#改善その3 入れてみる数字を絞る+「正しい状態」チェックの省略
改善その2では各マスに対して、
・入れることのできる数字(候補)のリストを出して
・その候補数の少ないマスを選んで数字入れを試す
ということを行なっていました。
この手順によれば、候補数の少ないマスがどこかだけではなく「そのマスに入れることのできる数字のリスト」を出すこともできるので、よく考えらたらわざわざ1から9までの数字を全て入れて試す必要はありません。
下図のオレンジ色のマスでいえば、1,2,4,6の4種類のみを試せば事足ります。1から9までの9種類を試す必要はありません。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/60b03a79-e8c3-0631-5bdc-cd51d82a6215.png" width="50%">
さらに、これにより「正しい状態」チェックを省略することが可能です。「入れることのできる数字」しか入れないのですから、その数字を入れた後は「正しい状態」になっていることが保証されているからです。~~改善その1の存在意義~~

実装は以下の通りになります。

```python:solve_sudoku4.py
sudoku = "800000000003600000070090200050007000000045700000100030001000068008500010090000400" #https://rocketnews24.com/2012/07/03/22654/

import numpy as np

def fi(s): #次にどこのマスに数字を入れればいいか(入る候補が最も少ないのはどこか)を判定します
    _ = []
    for k in range(81):
      if s[k] != "0":
        _ += [0]
      else:
        ss_yoko = [s[k//9*9+j] for j in range(9)]
        ss_tate = [s[k%9+9*j] for j in range(9)]
        ss_block = [s[k//27*27+k%9//3*3+j//3*9+j%3] for j in range(9)]
        _ += [len(set(ss_tate+ss_yoko+ss_block))]

    k = _.index(max(_))
    ss_yoko = [s[k//9*9+j] for j in range(9)]
    ss_tate = [s[k%9+9*j] for j in range(9)]
    ss_block = [s[k//27*27+k%9//3*3+j//3*9+j%3] for j in range(9)]
    return k,set(map(str,range(1,10)))-set(ss_tate+ss_yoko+ss_block) #入れることのできる数字の候補も出力しています

def dfs(i,t): #全探索を行う部分です
  if i == sudoku.count("0"):
    print(np.array(list(t)).reshape(9,9))
  else:
    k,next_number = fi(t)
    for j in next_number:
      dfs(i+1,t[:k]+j+t[k+1:])

print(np.array(list(sudoku)).reshape(9,9))
print("-----------------------------------------")
dfs(0,sudoku)
```

実行時間は5秒(答えが出力されるまで2秒)。改善その2と比べ約3割短縮されています(有効数字が怪しい。time関数を使えばちゃんと測れるがめんどいのでcolabで表示される実行時間です)

# 改善その4 入れられる数字を抽出する関数の高速化
改善その3では、入れられる候補の数を全ての空白マスに対して算出していています。
しかしながら、以下の図を見ればわかるとおり行、列またはブロックが同じならばそれらに含まれる数字のセットは共有することが可能です。
(例えば、下図の左と右で行(8)とブロック(3,7,8)の数字セットを共有できる)
数独の状態を文字列で管理していることから各マスの数字を呼び出すのにコストがかかるため、数字セットを共通化できれば計算量の削減が期待されます。
![スクリーンショット 2021-05-25 0.49.37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/c9652075-37aa-2477-56c4-cedbe8de79cb.png)

具体的には、初めに下図のように各行各列各ブロックに対してどの数字が使われているのかのセットをリストアップします(図では各ブロックに対するセットは省略。~~描きづらいので~~)
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/8710963b-f120-d9d3-314f-3419b90445d3.png" width="50%">

次に、全ての空白マスに対して、所属する行、列、ブロックのセットの和集合をとります。
下図のオレンジ色のマスで言えば、{5,7}と{4,9}と{1,4,5,7}の和集合をとり{1,4,5,7,9}となります。このセルに入れられる候補は残りの{2,3,6,8}の4種類です。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/391077/53664397-4257-3dca-5271-d85cddcd24cd.png" width="50%">


実装は以下の通りです。

```python:solve_sudoku5.py
sudoku = "800000000003600000070090200050007000000045700000100030001000068008500010090000400" #https://rocketnews24.com/2012/07/03/22654/

import numpy as np

def fi(s): #次にどこのマスに数字を入れればいいか(入る候補が最も少ないのはどこか)を判定します
    _ = []
    used_yoko = [set([s[k%9+l*9] for l in range(9)]) for k in range(9)] #あらかじめ、各列各行各ブロックで使用済みの数字をリストアップしておきます
    used_tate = [set([s[k*9+l] for l in range(9)]) for k in range(9)]
    used_block = [set([s[k//3*27+k%3*3+l//3*9+l%3] for l in range(9)]) for k in range(9)]
    for k in range(81):
      if s[k] != "0":
        _ += [0]
      else:
        ss = used_yoko[k%9]|used_tate[k//9]|used_block[k//27*3+k%9//3]
        if "0" in ss:
          ss.remove("0")
        _ += [len(set(ss))]
    k_next = _.index(max(_))
    ss_next =used_yoko[k_next%9]|used_tate[k_next//9]|used_block[k_next//27*3+k_next%9//3]
    return k_next,set(map(str,range(1,10)))-set(ss_next)

def dfs(i,t): #全探索を行う部分です
  if i == sudoku.count("0"):
    print(np.array(list(t)).reshape(9,9))
  else:
    k,next_number = fi(t)
    for j in next_number:
      dfs(i+1,t[:k]+j+t[k+1:])

print(np.array(list(sudoku)).reshape(9,9))
print("-----------------------------------------")
dfs(0,sudoku)
```
実行時間は3秒(答えが表示されるまで1秒)にまで短縮できました。

# まとめ
最初は「再帰関数を使った全探索で数独を解こう」と思っていただけでしたが、実行速度があまりにも遅かったので改善を行なったところ、面白い検証ができました。
単純な全探索では17分以上かかっていた数独を解くプログラムを、アルゴリズム改善により3秒にまで実行時間を縮めることが可能となったことから分かるように、アルゴリズムをしっかり考えることの重要性を実感しました。

# 追記(2024/6/13)
今までのプログラムでは、状態遷移の度に毎回数独の状態(81文字の文字列)を生成しており非効率的でした。
そこで数独の状態をlistで表現し、
 - 数字を一つ削除する
 - 数字を一つ挿入する

という二通りの遷移を高速で行えるようにしました。
また、再帰関数ではなくlistによるstackを用いた深さ優先探索としました(pypyだとこっちの方が早いけどpythonだとどうなんだろう)。

```python3solve_sudoku7.py
sudoku_str = "800000000003600000070090200050007000000045700000100030001000068008500010090000400"
sudoku = list(map(int,list(sudoku_str)))


def display(sudoku):
    for i in range(9):
        print(*sudoku[i*9:i*9+9])

def relation_listup():
    r = [[] for _ in range(81)]
    for i in range(81):
        x,y = i//9, i%9
        r[i].extend(list(range(y, y+81, 9)))
        r[i].extend(list(range(9*x, 9*x+9)))
        base_x, base_y = x//3*3, y//3*3
        r[i].extend(list(range(base_x*9+base_y, base_x*9+base_y+3)))
        r[i].extend(list(range(base_x*9+base_y+9, base_x*9+base_y+12)))
        r[i].extend(list(range(base_x*9+base_y+18, base_x*9+base_y+21)))
    for i in range(81):
        r[i] = sorted(list(set(r[i])))
    return r

cell_relationship = relation_listup()

insertable = [list(range(1,10)) for _ in range(81)]
for i in range(81):
    for cell in cell_relationship[i]:
        if sudoku[cell] in insertable[i]:
            insertable[i].remove(sudoku[cell])

tmp_min = 10
for i in range(81):
    if sudoku[i] == 0 and len(insertable[i]) < tmp_min:
        tmp_min = len(insertable[i])
        insert_num = insertable[i]
        insert_pos = i

zeros = sudoku.count(0)
stack = [(True, insert_pos, i,[], zeros) for i in insert_num]

cnt = 0
while stack:
    pre_order, insert_pos, insert_num, restore, zeros = stack.pop()
    if pre_order:
        sudoku[insert_pos] = insert_num
        if zeros == 1:
            cnt += 1
            display(sudoku)
        for cell in cell_relationship[insert_pos]:
            if insert_num in insertable[cell]:
                insertable[cell].remove(insert_num)
                restore.append(cell)
        tmp_min = 10
        for i in range(81):
            if sudoku[i] == 0 and len(insertable[i]) < tmp_min:
                tmp_min = len(insertable[i])
                next_insert_num = insertable[i]
                next_insert_pos = i
        stack.append((False, insert_pos, insert_num, restore, zeros+1))
        stack.extend([(True, next_insert_pos, nin,[], zeros-1) for nin in next_insert_num])
    else:
        sudoku[insert_pos] = 0
        for r in restore:
            insertable[r].append(insert_num)
```

→実行時間は0.3秒になりました(この問題において解が得られるまでには0.1秒)。
数字を挿入する場所をわざわざ全探索しているなど、まだ改善代はありそうです。