Problem 96: 数独

* 本記事は[ProjectEuler](https://projecteuler.net/about)の「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

（2025/05/20 更新「数独をExact cover問題として解く（DLX package 使用）」へのリンクを追加）

### 問題 96. 数独

原文 [Problem 96: Su Doku](https://projecteuler.net/problem=96)

**問題の要約：ファイルから数独の問題50問を読み込んで、その答えの左上の3桁の数字の合計を求めよ**

### PuLPを用いた数理モデルの解法

数独と解くプログラムは世の中にたくさんあると思われますが、その中でも印象的だった「PuLPを用いた数理モデルの解法」を紹介します。詳細は「[パズルでみる組合せ最適化のテクニック](https://qiita.com/SaitoTsutomu/items/f05f4ff05d4ce6099ff3)」等、[Tsutomu Saito](https://qiita.com/SaitoTsutomu)氏の投稿に詳しい説明がありますので参考にして見てください。

### アルゴリズムの簡単な説明

数理モデルの考え方はちょっと独特なので慣れが必要です。ここでは簡単な例として$2 \times 2$の数独を考えてみます。入る数字は１と２のみで各行・列に１と２があるように配置するという超簡単な例です。

| 1| -|
| ----  | --- |
| **-** | **-**|

これに対して以下の8個の変数を作ります。例えば変数V_0_0_1は0行0列の数が1になる確率とみなし、値が1ならば数は1、値が0ならば1以外と考えます。

| 変数 |  1になるのは|
| ----  | --- |
|V_0_0_1 | 0行0列の数が1|
|V_0_0_2 | 0行0列の数が2|
|V_0_1_1 | 0行1列の数が1 |
|V_0_1_2 | 0行1列の数が2 |
|V_1_0_1 | 1行0列の数が1  |
|V_1_0_2 | 1行0列の数が2 |
|V_1_1_1 | 1行1列の数が1 |
|V_1_1_2 | 1行1列の数が2 |

これらの変数に対して条件を連立方程式にしていきます。

```
#--- 行のルール
V_0_0_1+V_0_1_1 == 1    # 0行目には1は1つしかない
V_0_0_2+V_0_1_2 == 1    # 0行目には2は1つしかない
V_1_0_1+V_1_1_1 == 1    # 1行目には1は1つしかない
V_1_0_2+V_1_1_2 == 1    # 1行目には2は1つしかない
#--- 列のルール
V_0_0_1+V_0_0_1 == 1    # 0列目には1は1つしかない
V_1_0_2+V_1_0_2 == 1    # 0列目には2は1つしかない
V_0_1_1+V_0_1_1 == 1    # 1列目には1は1つしかない
V_1_1_2+V_1_1_2 == 1    # 1列目には2は1つしかない
#--- 一つのマスには一つの数
V_0_0_1+V_0_0_2 == 1    # 0行0列には1か2かどちらかが入る
V_0_1_1+V_0_1_2 == 1    # 0行1列には1か2かどちらかが入る
V_1_0_1+V_1_0_2 == 1    # 1行0列には1か2かどちらかが入る
V_1_0_1+V_1_1_2 == 1    # 1行1列には1か2かどちらかが入る
#---
V_0_0_1 == 1            # 0行0列が1になる確率は１
```
これらの連立方程式を**solve**に渡すと答えを返してくれるという仕組みです。

数独はこれを$9 \times 9$にして$3 \times 3$のマスのルールを足したものなので、以下のようなコードになります。数独を解くプログラムがこれだけ簡潔に書けるのは凄いと思います。

* 変数は**addvinvars**を使って自動的に生成している
* 連立方程式の足し算は**lpSum**関数を使う。
* 結果の値は**value**関数をnumpyのarrayに**vectorize**関数を適用して取り出します

Google Colabで**PuLP**を使うためにはインストールが必要なので下の方にコードを添付しました。ファイルの読み込みも添付してあります。

```python 
from ortoolpy import addbinvars
from pulp import LpProblem, LpVariable, LpBinary, lpSum, lpDot, value
import numpy as np
import itertools 

def SuDoku(sd):
  m = LpProblem()
  gd = np.array(addbinvars(9,9,9))  # Variable array
  for i,j in itertools.product(range(9),range(9)):
    m += lpSum(gd[:,i,j]) == 1    # all numbers in a column
    m += lpSum(gd[i,:,j]) == 1    # all numbers in a row
    m += lpSum(gd[i,j,:]) == 1    # 1 number in the grid
    k, l = i//3*3, i%3*3
    m += lpSum(gd[k:k+3,l:l+3,j]) == 1   # all numbers in 3x3 
    c = sd[i*9+j]
    if int(c)>0:
      m += gd[i,j,int(c)-1] == 1   # if a number already assigned
  m.solve()                        # Call solver
  ret = np.vectorize(value)(gd).astype(int).dot(range(1,10))
  return "".join(map(str, list(itertools.chain.from_iterable(ret)))) # flatten the result

print(f"Answer: {sum([int(SuDoku(sd)[:3]) for sd in sdata])}")
```

### PuLpとortoolpyのインストール

```Python 
import pip
!pip install ortoolpy
!pip install pulp
```
### 数独問題ファイルの読み込み

```python 
# show upload dialog
from google.colab import files
uploaded = files.upload()
#---- read numbers from file to name[]
f = open("p096_sudoku.txt")
Flines = f.readlines()
f.close()
sdata = ["".join([Flines[g*10+i].strip() for i in range(1,10)]) for g in range(50)]
print(sdata[0])
```

### 数独をExact cover問題として解く

[「数独をExact cover問題として解く（DLX package 使用）」](https://qiita.com/masa0599/items/78ea14db3beb10966730) に投稿しました。こちらも参考にしてください。

（開発環境：Google Colab）