[麻雀]シャンテン数計算アルゴリズム

# はじめに
本記事では麻雀におけるシャンテン数の計算アルゴリズムについて解説します。このアルゴリズムの特徴として、計算コストが枚数およびシャンテン数によらず高速に計算できることが挙げられます。後半にサンプルプログラムについて説明していますので、早くこのアルゴリズムを試してみたい方は先にご覧になることをおすすめします([サンプルプログラム](#サンプルプログラム))。なお、本記事では厳密さを重視するあまり、もはや麻雀とは思えない議論が展開されていますので注意してください。

# 背景
麻雀では手の進み具合をシャンテン数という言葉で表します。人によって表現の仕方は違いますが、「テンパイするまでに必要な最小のツモ回数」とここでは定義します。

シャンテン数を計算するには七対子と国士無双、四面子一雀頭の3つあがり形に分けて考える必要があります。七対子と国士無双の場合は簡単で、それぞれトイツの数と一九字牌の数を数えればよいです。

```math
\begin{equation}
(七対子のシャンテン数)=6-(トイツ数)
\end{equation}
\tag{1}
```

※ただし牌の種類が7未満の場合はさらに(7-(牌の種類))を引く

```math
\begin{equation}
(国士無双のシャンテン数)=13-(一九字牌の種類数)
\end{equation}
\tag{2}
```

※ただし一九字牌のトイツがある場合はさらに1を引く

少し難しいのが四面子一雀頭の場合です。一般的な方法では、手牌をメンツ、トイツ、ターツに分解してそれらの数から求めます。

```math
\begin{equation}
(四面子一雀頭形のシャンテン数)=8-2(メンツ数)-(トイツ数)-(ターツ数)
\end{equation}
\tag{3}
```

※ただし雀頭がある場合はさらに1を引く

多くの場合手牌の分解パターンは複数存在するので、それらの中から最小のシャンテン数を与えるパターンを採用することになります。ただし、上の式を使うには手牌の分解パターンが確かに和了形を与えることを確認する必要があります。例えば、ターツをシュンツにするために5枚目の牌を要求していないか確認するといった具合です。特に染め手や四枚遣いが含まれる手ではこのような確認を行うことが困難になります。そこでこの問題を回避するために、本記事では独自のアプローチでシャンテン数を求めてみようと思います。

# 手牌の表記
以降では手牌を${\bf h}$のように表し、$i$番目の牌の枚数を$h_i (0 \le i \le 33)$で表します。牌番号の対応は次の表のようになります。

|牌番号|牌|
|---|---|
|0|1マン|
|...|...|
|8|9マン|
|9|1ピン|
|...|...|
|17|9ピン|
|18|1ソウ|
|...|...|
|26|9ソウ|
|27|東|
|28|南|
|29|西|
|30|北|
|31|白|
|32|発|
|33|中|

# 距離
シャンテン数を求める準備として距離という概念を導入します。

手牌${\bf h}$と手牌${\bf g}$について、手牌${\bf g}$から手牌${\bf h}$に変換するためのツモ回数を「距離」として、以下のように定義します。

```math
d({\bf h}, {\bf g}) = \frac{1}{2} \sum_{i=0}^{33} (|h_i - g_i|+h_i - g_i) \tag{4}
```

# シャンテン数を求める方針
冒頭でシャンテン数は「テンパイするまでに必要な最小のツモ回数」と述べました。これを言い換えれば「あがりまでに必要な最小のツモ回数-1」となります。よって、正確にシャンテン数を求めたければ、与えられた手牌と全てのあがり形の距離を求め、その最小値-1をもってシャンテン数とすればよいです。

# 置換数とシャンテン数の定義
全てのあがり形の集合を$W$としたとき、その各要素${\bf w}$と手牌${\bf h}$の距離の最小値を置換数$T({\bf h})$として、以下のように定義します。

```math
T({\bf h}) = \min_{{\bf w} \in W} d({\bf w}, {\bf h})
{\tag 5}
```

シャンテン数$S({\bf h})$は置換数$T({\bf h})$を用いて

```math
S({\bf h}) = T({\bf h})-1
\tag{6}
```

のように定義します。以上でシャンテン数が厳密に定義されました。この定義は四面子一雀頭、七対子、国士無双のあがり形によらずに適用することができます。

# 置換数(シャンテン数)の計算方法
シャンテン数を計算するには(5)式を直接計算できればよいのですが、$W$の要素数は四面子一雀頭形で11498658と膨大なため（[参考](https://qiita.com/tomo_hxx/items/96a4efaacfbb6761cc59)）、多大な時間とメモリが必要になってしまいます。そこで置換数を求める効率的なアルゴリズムを紹介します。ここから四面子一雀頭形に限定して議論を進めます。

議論を簡単にするため${\bf h}$が2色の数牌からなる手牌だとすると、最も近いあがり形は明らかに同じ2色からなる手牌です。このとき、このあがり形の集合$W$は各色の$n$組からなるメンツまたは$n$組のメンツと一雀頭からなる手牌の集合の__直積__となります。ここで$n$組のメンツの集合を$V_n$、$n$組のメンツと一雀頭の集合を$W_n$とすると、

```math
\begin{equation}
W = (V^0_0 \times W^1_4) \cup (V^0_1 \times W^1_3) \cup ... \cup (V^0_4 \times W^1_0) \\\
\cup (W^0_0 \times V^1_4) \cup (W^0_1 \times V^1_3) \cup ... \cup (W^0_4 \times V^1_0)
\end{equation}
\tag{7}
```

となります($\times$は集合の直積を表します)。ここで上付き添え字は何番目の色かを表しています。適宜0をマンズ、1をピンズ、2をソーズ、3を字牌のように読み替えてください。例えば上付き添え字が0のときはマンズ以外の牌が手牌にないことを意味します。(7)式より(5)式は以下のように書き表されます。

```math
\begin{align}
T({\bf h}) &= \min_{{\bf w} \in W} d({\bf w}, {\bf h}) \\\
&= \min \left[
\min_{{\bf v^0} \in V^0_0, {\bf w^1} \in W^1_4} d({\bf v^0}+{\bf w^1}, {\bf h}), \min_{{\bf v^0} \in V^0_1, {\bf w^1} \in W^1_3} d({\bf v^0}+{\bf w^1}, {\bf h}) , ... , \min_{{\bf v^0} \in V^0_4, {\bf w^1} \in W^1_0} d({\bf v^0}+{\bf w^1}, {\bf h}), \\\
\min_{{\bf w^0} \in W^0_0, {\bf v^1} \in V^1_4} d({\bf w^0}+{\bf v^1}, {\bf h}), \min_{{\bf w^0} \in W^0_1, {\bf v^1} \in V^1_3} d({\bf w^0}+{\bf v^1}, {\bf h}) , ... , \min_{{\bf w^0} \in W^0_4, {\bf v^1} \in V^1_0} d({\bf w^0}+{\bf v^1}, {\bf h})
\right]
\end{align}
\tag{8}
```

ここで(8)式の大括弧の中の各項に注目します。例えば第一項は、

```math
\begin{align}
\min_{{\bf v^0} \in V^0_0, {\bf w^1} \in W^1_4} d({\bf v^0}+{\bf w^1}, {\bf h}) &= \min_{{\bf v} \in V^0_0, {\bf w} \in W^1_4} \frac{1}{2} \sum_{i=0}^{34} (|v_i + w_i- h_i|+v_i + w_i - h_i) \\\
&= \min_{{\bf v} \in V^0_0, {\bf w} \in W^1_4} \left[ \frac{1}{2} \sum_{i=0}^{8} (|v_i - h_i|+v_i - h_i) + \frac{1}{2} \sum_{i=9}^{17} (|w_i - h_i|+w_i - h_i) \right] \\\
&= \min_{{\bf v} \in V^0_0} \frac{1}{2} \sum_{i=0}^{8} (|v_i - h_i|+v_i - h_i) + \min_{{\bf w} \in W^1_4} \frac{1}{2} \sum_{i=9}^{17} (|w_i - h_i|+w_i - h_i) \\\
&\equiv u^0_0 + t^1_4
\end{align}
\tag{9}
```

と書くことができます。最後の式変形で新たに$t^0_0$と$u^1_4$という表記を導入しました。ここで__部分置換数__を以下のように定義します。

```math
\begin{align}
u^n_m &= \min_{{\bf v} \in V^n_m} \frac{1}{2} \sum_{i=9n}^{9n+8(6)} (|v_i - h_i|+v_i - h_i) & (0 \le n \le 3, 0 \le m \le 4) \\\
t^n_m &= \min_{{\bf w} \in W^n_m} \frac{1}{2} \sum_{i=9n}^{9n+8(6)} (|w_i - h_i|+w_i - h_i) & (0 \le n \le 3, 0 \le m \le 4)
\end{align}
\tag{10}
```

$9n+8(6)$という表記は、$n=3$すなわち字牌に対しては$9n+6$に読み替えるという意味です。

(8)式の他の項についても同様に書き換えると、

```math
\begin{equation}
T({\bf h}) = \min \left[ t^0_0+u^1_4, t^0_1+u^1_3, ..., t^0_4+u^1_0, \\\
u^0_0+t^1_4, u^0_1+t^1_4, ..., u^0_4+t^1_0 \right]
\end{equation}
\tag{11}
```

となります。つまり、__全ての牌の組合わせに対して__予め部分置換数$t^n_m$と$u^n_m$を計算し(ファイルに保存し置換数の計算前にメモリに格納し)ておけば、__(11)式を計算することで(5)式を直接計算したのと同じ結果を得ることができます。__なお、このときに必要なメモリサイズは

```math
\begin{equation}
(牌の組合わせ) \times (t^n_mとu^n_mの要素数) \times ({\rm int}型のサイズ) = 5^9 \times 10 \times 4 = 78 {\rm MB}
\end{equation}
\tag{12}
```

であり、十分現実的なサイズといえます。

ここで(11)式の意味について考えてみます。2色の手牌のあがり形として、枚数の組合せが$(0,14),(3,11), ...,(12,2),(2,12),(5,9), ...,(14,0)$の10通りの場合が考えられます。その各々の場合に対して、手牌をその枚数のメンツ(と一雀頭)の組に変えるのに必要なツモ数を計算します。求められた10個の値の内、最小値を置換数として求めているのです。

4色の手牌の置換数(シャンテン数)を求めるためには(11)式の計算を繰り返し実行すれば良いです。動的計画法を用います。

```math
\begin{align}
t^{(n+1)}_m &= \min_{0 \le l \le m} \left[ t^{(n)}_l+u^{n+1}_{m-l}\;,\; u^{(n)}_l+t^{n+1}_{m-l} \right] & (0 \le n \le 3\;,\; 0 \le m \le 4)
\tag{13} \\\
u^{(n+1)}_m &= \min_{0 \le l \le m} \left[ u^{(n)}_l+u^{n+1}_{m-l} \right] & (0 \le n \le 3\;,\; 0 \le m \le 4)
\tag{14}
\end{align}
```

```math
\begin{align}
t^{(0)}_m &= t^0_m & (0 \le m \le 4) \\\
u^{(0)}_m &= u^0_m & (0 \le m \le 4)
\end{align}
\tag{15}
```

このとき置換数は

```math
\begin{equation}
T({\bf h}) = t^{(3)}_4
\end{equation}
\tag{16}
```

となります。$t^{(3)}_3, ..., t^{(3)}_0, u^{(3)}_4, ..., u^{(3)}_0$を計算する必要はありません。

(13)式と(14)式は次のように書き換えることができます。

```math
\begin{align}
t^{(n+1,\:l+1)}_m &= \min \left[ \min \left[ t^{(n+1,\:l)}_m\;,\;t^{(n)}_{l+1}+u^{n+1}_{m-l-1} \right]\;,\; u^{(n)}_{l+1}+t^{n+1}_{m-l-1} \right] \\\
t^{(n+1,\:0)}_m &= \min \left[ t^{(n)}_0+u^{n+1}_m \;,\; u^{(n)}_0+t^{n+1}_m \right] \\\
t^{(n+1)}_m &= t^{(n+1,\:m)}_m
\end{align}
(0 \le n \le 3\;,\; 0 \le m \le 4\;,\; 0 \le l < m)
\tag{17}
```

```math
\begin{align}
u^{(n+1,\:l+1)}_m &= \min \left[ u^{(n+1,\:l)}_m\;,\;u^{(n)}_{l+1}+u^{n+1}_{m-l-1} \right] \\\
u^{(n+1,\:0)}_m &= u^{(n)}_0+u^{n+1}_m \\\
u^{(n+1)}_m &= u^{(n+1,\:m)}_m
\end{align}
(0 \le n \le 3\;,\; 0 \le m \le 4\;,\; 0 \le l < m)
\tag{18}
```

次のソースコードは(17)式と(18)式に基づいて実装しました。

# ソースコード
(13)式の計算を行うC++コードは以下のようになります。事前に部分置換数を(10)式により計算しておきます。

```cpp
/*
lhs[5]からlhs[9]にt^{(n)}_0からt^{(n)}_4を, lhs[0]からlhs[4]にu^{(n)}_0からu^{(n)}_4をセットする。
rhs[5]からrhs[9]にt^{n+1}_0からt^{n+1}_4を, lhs[0]からlhs[4]にu^{n+1}_0からu^{n+1}_4をセットする。
lhs[5]からlhs[9]をt^{(n+1)}_0からt^{(n+1)}_4で, lhs[0]からlhs[4]をu^{(n+1)}_0からu^{(n+1)}_4で置き換える。
*/

using Vec  = std::vector<int>;

void add(Vec& lhs, const Vec& rhs)
{
  //t_4からt_0について計算
  for(int j=9; j>=5; --j){
    int sht = std::min(lhs[j]+rhs[0], lhs[0]+rhs[j]);
    
    for(int k=5; k<j; ++k){
      sht = std::min({sht, lhs[k]+rhs[j-k], lhs[j-k]+rhs[k]});
    }
    lhs[j] = sht;
  }

  //u_4からu_0について計算
  for(int j=4; j>=0; --j){
    int sht = lhs[j]+rhs[0];
    
    for(int k=0; k<j; ++k){
      sht = std::min(sht, lhs[k]+rhs[j-k]);
    }
    lhs[j] = sht;
  }
}
```


# サンプルプログラム
以上のアルゴリズムを用いて、シャンテン数を計算するクラスと、配牌時のシャンテン数をモンテカルロ法で計算するサンプルプログラムをC++で作成しました。詳しい使い方はリンク先を参照してください。
[shanten-number-calculator](https://github.com/tomohxx/shanten-number-calculator) (GitHub)
サンプルプログラムの実行例として、親の配牌時シャンテン数と子の配牌時シャンテン数を計算してみます。

- ビルド

```shell
$ make
```

- 親の配牌時シャンテン数(1億回シミュレーション)

```shell
$ ./sample.out 14 100000000
=========================RESULT=========================
-1  315         0.000315
0   69900       0.0699
1   2333875     2.33387
2   19496567    19.4966
3   43932086    43.9321
4   28515676    28.5157
5   5496180     5.49618
6   155401      0.155401
Number of Tiles         14
Total                   100000000
Time (msec.)            99710
Expected Value          3.15599
```

- 子の配牌時シャンテン数(1億回シミュレーション)

```shell
./sample.out 13 100000000
=========================RESULT=========================
-1  0           0
0   8196        0.008196
1   621271      0.621271
2   9357183     9.35718
3   36197590    36.1976
4   39876941    39.8769
5   13102859    13.1029
6   835960      0.83596
Number of Tiles         13
Total                   100000000
Time (msec.)            90403
Expected Value          3.57966
```

親の配牌時シャンテン数の期待値は3.15599、子の配牌時シャンテン数の期待値は3.57966と計算されました。ちなみに厳密な値は親で3.15593、子で3.57967です(小数点以下5桁目まで)。

# おわりに
~~次は有効牌を計算するアルゴリズムについて書きたいと思います。~~
有効牌を計算するアルゴリズムについては[こちら](https://github.com/tomohxx/necessary-and-unnecessary-tiles-calculator)をご覧ください。ソースコードのみで解説はありません。