# 数独

このページでは、**Hi-QUBO** を使って数独パズルをイジングマシン（アニーリング）によって解く方法を紹介します。  

---

## 数独とは

数独（Sudoku）は、**9×9 グリッドに 1 〜 9 の数字を埋めるパズル**です。  
次のルールを満たすように空欄に数字を入れていきます。

- 各行には 1 〜 9 の数字が一度ずつ現れること  
- 各列には 1 〜 9 の数字が一度ずつ現れること  
- 3×3 のブロック内にも 1 〜 9 の数字が一度ずつ現れること  

ヒントとして初期値がすでに何個か埋まっており、その状態から全てのマスを埋めると完成です。

---

## イジングマシンによる解法の考え方

プログラムでは、数独のルールを **制約条件に変換**し、イジングマシンに解かせます。  
各セルに入る数字は「バイナリ変数（0/1）」として表現し、制約を満たす組み合わせを最適化によって見つけます。

---

## 制約条件の定式化

### 変数定義

- $y_{i,j,k} \in \{0,1\}$ を 3 次元配列で定義  
- $i,j=0,...,8$ が行・列、$k=0,...,8$ が数字（1〜9）を表す  
- $y_{i,j,k}=1$ のとき、セル $(i,j)$ は $k+1$ が入ることを意味する  

これにより変数が **$9 \times 9 \times 9 = 729$ 個** 作られます。

---

### ルールからの制約

数独のルールは次のような one-hot 制約として定義できます：

- **制約 (a)**: 1 マスには 1 個の数字しか入らない  

$$
\sum_{k=0}^8 y_{i,j,k} = 1
$$

- **制約 (b)**: 同じ行に同じ数字は 1 個

$$
\sum_{j=0}^8 y_{i,j,k} = 1
$$

- **制約 (c)**: 同じ列に同じ数字は 1 個 

$$
\sum_{i=0}^8 y_{i,j,k} = 1
$$

- **制約 (d)**: 同じ $3 \times 3$ ブロックに同じ数字は 1 個 

$$
\sum_{(i,j)\in \text{block}} y_{i,j,k} = 1
$$


---

## コード例

以下はサンプルコードです。
:::{container} prog-cpp
```cpp
#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>

int main() {
    const size_t num_size = 3;
    const size_t use_num = num_size*num_size;
    const size_t N = use_num;

    // 初期配置を配列で表記
    std::vector<std::vector<size_t>> initial = 
        {
            {5, 3, 0, 0, 7, 0, 0, 0, 0},
            {6, 0, 0, 1, 9, 5, 0, 0, 0},
            {0, 9, 8, 0, 0, 0, 0, 6, 0},
            {8, 0, 0, 0, 6, 0, 0, 0, 3},
            {4, 0, 0, 8, 0, 3, 0, 0, 1},
            {7, 0, 0, 0, 2, 0, 0, 0, 6},
            {0, 6, 0, 0, 0, 0, 2, 8, 0},
            {0, 0, 0, 4, 1, 9, 0, 0, 5},
            {0, 0, 0, 0, 8, 0, 0, 7, 9}
    };

    // 変数の定義
    auto y = qbpp::var("y", N, N, use_num);

    // 初期配置を変数に反映するための配列
    qbpp::MapList ml;
    for (size_t i=0; i<N; i++){
        for (size_t j=0; j<N; j++){
            size_t k = initial[i][j];
            if ( k != 0) {
                ml.push_back({y[i][j][k-1], 1});
            }
        }  
    }

    // 制約 (a): 1マスには1個の数字しか入らない
    auto c_num_select = qbpp::sum(qbpp::vector_sum(y) == 1);

    // 1~9が、各行と列とブロックの中に何回現れているか数える変数
    auto row = qbpp::vector_sum(y, 1);
    auto column = qbpp::vector_sum(y, 0);
    auto block = qbpp::expr(use_num, N);
    for(size_t k=0; k<use_num; k++){
        for(size_t i=0; i<N; i++){
            size_t i_q = i/num_size;
            size_t i_r = i%num_size;
            for(size_t j=0; j<N; j++){
                size_t j_q = j/num_size;
                size_t j_r = j%num_size;
                block[k][i] += y[num_size*i_q+j_q][num_size*i_r+j_r][k];
            }
        }
    }

    // 制約 (b): 同じ行に同じ数字は 1 個 
    auto c_row = qbpp::sum(row==1);
    // 制約 (c): 同じ列に同じ数字は 1 個 
    auto c_column = qbpp::sum(column==1);
    // 制約 (d): 同じブロックに同じ数字は 1 個 
    auto c_block = qbpp::sum(block==1);

    // 実際の数独の配列を作成する
    auto num = qbpp::expr(use_num);
    for (size_t k=0; k<use_num; k++){
        num[k] = qbpp::toExpr(k+1);
    }

    auto x = qbpp::expr(N, N);
    for(size_t i=0; i<N; i++){
        for(size_t j=0; j<N; j++){
            x[i][j] = qbpp::sum(num*y[i][j]);
        }
    }

    // QUBO式の作成
    auto f = c_row + c_column + c_block + c_num_select;
    // 初期配列をQUBO式に反映(該当する要素を定数化)
    auto g = qbpp::replace(f, ml);
    // QUBO式の展開
    f.simplify_as_binary();
    g.simplify_as_binary();
    // QUBO式を解く
    auto solver = qbpp::easy_solver::EasySolver(g);
    solver.time_limit(30);
    solver.target_energy(0);
    auto sol = solver.search();
    // 初期条件を反映させた解を作成
    qbpp::Sol full_sol(f);
    full_sol.set(sol);
    full_sol.set(ml);

    // 結果の表示
    const std::string hline = "+-------+-------+-------+";

    for (size_t i = 0; i < N; i++) {
        if (i % num_size == 0) {
            std::cout << hline << std::endl;
        }

        for (size_t j = 0; j < N; j++) {
            if (j % num_size == 0) {
                std::cout << "| ";
            }
            std::cout << x[i][j](full_sol) << " ";
        }
        std::cout << "|" << std::endl;
    }
    std::cout << hline << std::endl;
  }
```
:::


:::{container} prog-python
```python
import pyqbpp as qbpp
num_size = 3
use_num = num_size*num_size
N = use_num

# 初期配置を配列で表記
initial = [
                [5, 3, 0, 0, 7, 0, 0, 0, 0],
                [6, 0, 0, 1, 9, 5, 0, 0, 0],
                [0, 9, 8, 0, 0, 0, 0, 6, 0],
                [8, 0, 0, 0, 6, 0, 0, 0, 3],
                [4, 0, 0, 8, 0, 3, 0, 0, 1],
                [7, 0, 0, 0, 2, 0, 0, 0, 6],
                [0, 6, 0, 0, 0, 0, 2, 8, 0],
                [0, 0, 0, 4, 1, 9, 0, 0, 5],
                [0, 0, 0, 0, 8, 0, 0, 7, 9]
        ]

# 変数の定義
y = qbpp.var("y", N, N, use_num)

# 初期配置を変数に反映するための配列
ml = qbpp.MapList()
for i in range(N):
    for j in range(N):
        k = initial[i][j]
        if k != 0:
            ml.add(y[i][j][k-1], 1)

# 制約 (a): 1マスには1個の数字しか入らない
c_num_select = qbpp.sum(qbpp.vector_sum(y) == 1)

# 1~9が、各行と列とブロックの中に何回現れているか数える変数
row = qbpp.vector_sum(y, 1)
column = qbpp.vector_sum(y, 0)
block = qbpp.expr(use_num, N)
for k in range(use_num):
    for i in range(N):
        i_q = i//num_size
        i_r = i%num_size
        for j in range(N):
            j_q = j//num_size
            j_r = j%num_size
            block[k][i] += y[num_size*i_q+j_q][num_size*i_r+j_r][k]

# 制約 (b): 同じ行に同じ数字は 1 個 
c_row = qbpp.sum(row==1)
# 制約 (c): 同じ列に同じ数字は 1 個 
c_column = qbpp.sum(column==1)
# 制約 (d): 同じブロックに同じ数字は 1 個 
c_block = qbpp.sum(block==1)

# 実際の数独の配列を作成する
num = qbpp.expr(use_num)
for k in range(use_num):
    num[k] = qbpp.toExpr(k+1)

x = qbpp.expr(N, N)
for i in range(N):
    for j in range(N):
        x[i][j] = qbpp.sum(num*y[i][j])

# QUBO式の作成
f = c_row + c_column + c_block + c_num_select
# 初期配列をQUBO式に反映(該当する要素を定数化)
g = qbpp.replace(f, ml)
# QUBO式の展開
f.simplify_as_binary()
g.simplify_as_binary()
# QUBO式を解く
solver = qbpp.EasySolver(g)
solver.time_limit(60)
solver.target_energy(0)
sol = solver.search()
# 初期条件を反映させた解を作成
full_sol = qbpp.Sol(f).set(sol)
full_sol = full_sol.set(ml)

# 結果の表示
hline = "+-------+-------+-------+"

for i in range(N):
    if i % num_size == 0:
        print(hline)
    for j in range(N):
        if j % num_size == 0:
            print("| ", end="")
        print(f"{full_sol(x[i][j])}", end=" ")
    print("|")
print(hline)
```
:::

## 実行結果の例
```
+-------+-------+-------+
| 5 3 4 | 6 7 8 | 9 1 2 |
| 6 7 2 | 1 9 5 | 3 4 8 |
| 1 9 8 | 3 4 2 | 5 6 7 |
+-------+-------+-------+
| 8 5 9 | 7 6 1 | 4 2 3 |
| 4 2 6 | 8 5 3 | 7 9 1 |
| 7 1 3 | 9 2 4 | 8 5 6 |
+-------+-------+-------+
| 9 6 1 | 5 3 7 | 2 8 4 |
| 2 8 7 | 4 1 9 | 6 3 5 |
| 3 4 5 | 2 8 6 | 1 7 9 |
+-------+-------+-------+
```