[麻雀]天和の確率計算

#はじめに
麻雀における天和の確率は一度は気になることだと思います。天和の確率はらすかるさんがすでに計算されていて約33万分の1です。計算方法はこちらです。

麻雀の数学

本記事ではらすかるさんの考え方を参考に、天和の確率計算アルゴリズムを解説します。

アルゴリズム

麻雀には四面子一雀頭、七対子、国士無双の3つのあがりパターンがあります。これらのパターンごとにあがりの組合せの数を計算して、その和を14枚の牌の組合せの数で割れば天和の確率を求めることができます。四面子一雀頭形を除くあがりパターンと14枚の牌の組合せの総数は簡単に表せます。

• 七対子

{}_{34}C_{7} \times ({}_4C_2)^{7}
\tag{1}

• 国士無双

13 \times {}_4C_2 \times 4^{12}
\tag{2}

• 14枚の牌の組合せ

{}_{136}C_{14}
\tag{3}

四面子一雀頭形の組合せの数は上記のような簡単な式では表せませんのでプログラムで計算します。

#ソースコード

tenhou_prob.cpp

#include <iostream>
#include <iomanip>

const int S = 3;//数牌の種類数
const int H = 7;//字牌の種類数
const int T = 9*S+H;//全ての牌の種類数
const int M = 14;//手牌の枚数

/*
i枚の一色の数牌の組合せに対して、
num[i][0]は全パターン数、
num[i][1]はi/3面子形のパターン数(i mod 3 = 0のとき)と(i-2)/3面子一雀頭形のパターン数(i mod 3 = 2のとき)、
num[i][2]はi/2対子形のパターン数(i mod 2 = 0のとき)
*/
unsigned long long num[M+1][3];

/*
i枚の一色の数牌の組合せに対して、
den[i][0]は全状態数、
den[i][1]はi/3面子形の状態数(i mod 3 = 0のとき)と(i-2)/3面子一雀頭形の状態数(i mod 3 = 2のとき)、
den[i][2]はi/2対子形の状態数(i mod 2 = 0のとき)
*/
unsigned long long den[M+1][3];

unsigned long long ntot  = 0;//(和了していない形を含めた)全パターン数
unsigned long long nlh   = 0;//四面子一雀頭形のパターン数
unsigned long long ntsis = 0;//四面子一雀頭かつ七対子である形(二盃口)のパターン数
unsigned long long nsp   = 0;//(二盃口を含む)七対子のパターン数
unsigned long long nto   = 0;//国士無双のパターン数
unsigned long long dtot  = 0;//(和了していない形を含めた)全状態数
unsigned long long dlh   = 0;//四面子一雀頭形の状態数
unsigned long long dtsis = 0;//四面子一雀頭かつ七対子である形(二盃口)の状態数
unsigned long long dsp   = 0;//(二盃口を含む)七対子のパターン数
unsigned long long dto   = 0;//国士無双のパターン数

//組み合わせの総数nCrの計算
unsigned long long combin(int n, int r)
{
  if(!r) return 1;
  else return (n-r+1)*combin(n,r-1)/r;
}

//nのr乗の計算
unsigned long long pow(int n, int r)
{
  if(!r) return 1;
  else return n*pow(n,r-1);
}

//長さ9のint型配列tに対して面子分解可能か判定する
bool iswh0(const int* t)
{
  int r, a = t[0], b = t[1];

  for(int i=0; i<7; i++){
    if(r=a%3, b>=r && t[i+2]>=r){a=b-r; b=t[i+2]-r;}
    else return false;
  }
  if(a%3==0 && b%3==0) return true;
  else return false;
}

//長さ9のint型配列に対してN面子一雀頭に分解可能か判定する
bool iswh2(int* t)
{
  int p = 0;

  for(int i=0; i<9; i++) p += i*t[i];

  for(int i=p*2%3; i<9; i+=3){
    if(t[i]>=2){
      if(t[i]-=2, iswh0(t)){t[i]+=2; return true;}
      else t[i]+=2;
    }
  }
  return false;
}

//長さ34のint型配列に対して七対子であるか判定する
bool issp(const int* t)
{
  for(int i=0; i<34; i++){
    if(t[i]!=0 && t[i]!=2) return false;
  }
  return true;
}

/*
再帰により一色の数牌からなるm枚以下の牌の組を生成し、num[][]とden[][]に値をセットする
*/
void dealtile(int n, int m, int l, int d)
{
  static int hd[34] = {};
  static const int c[5] = {1,4,6,4,1};
  static int sum[10] = {};

  if(n==9){
    ++num[sum[9]][0];
    den[sum[9]][0] += d;

    if(sum[9]%3==0 && iswh0(hd)){
      ++num[sum[9]][1];
      den[sum[9]][1] += d;

      if(sum[9]%2==0 && issp(hd)){
        ++num[sum[9]][2];
        den[sum[9]][2] += d;
      }
    }
    else if(sum[9]%3==2 && iswh2(hd)){
      ++num[sum[9]][1];
      den[sum[9]][1] += d;

      if(sum[9]%2==0 && issp(hd)){
        ++num[sum[9]][2];
        den[sum[9]][2] += d;
      }
    }
  }
  else{
    for(int i=0; i<=(m>=4 ? 4:m); ++i){
      hd[n]=i;
      sum[n+1] = sum[n]+i;

      dealtile(n+1,m-i,l-4,d*c[i]);
    }
  }
}

/*
m枚の牌に対する全パターン数と状態数および、四面子一雀頭形と二盃口のパターン数と状態数を数え上げる
*/
void dealsum(int n, int m, int l)
{
  static const int c[5][3] = {{1,1,1},{1,0,0},{1,1,1},{1,1,0},{1,0,0}};//c[i][j]は字牌がi枚あるときのパターン数および、i/3面子(一雀頭)形とi/2対子形のパターン数
  static const int d[5][3] = {{1,1,1},{4,0,0},{6,6,6},{4,4,0},{1,0,0}};//c[i][j]は字牌がi枚あるときの状態数および、i/3面子(一雀頭)形とi/2対子形の状態数
  static int res[S+H+1][3] = {0};//各色の枚数について3の剰余をカウント
  static unsigned long long p[S+H+1] = {1};//p[i]はi枚の牌の全パターン数
  static unsigned long long q[S+H+1] = {1};//q[i]はi枚の牌のi/3面子(一雀頭)形のパターン数
  static unsigned long long r[S+H+1] = {1};//r[i]はi枚の牌のi/2対子形のパターン数
  static unsigned long long s[S+H+1] = {1};//s[i]はi枚の牌の全状態数
  static unsigned long long t[S+H+1] = {1};//t[i]はi枚の牌のi/3面子(一雀頭)形の状態数
  static unsigned long long u[S+H+1] = {1};//u[i]はi枚の牌のi/2対子形の状態数

  if(n==S+H){
    ntot += p[S+H];
    dtot += s[S+H];

    if(res[S+H][1]==0 && res[S+H][2]==1){
      nlh += q[S+H];
      dlh += t[S+H];
      ntsis += r[S+H];
      dtsis += u[S+H];
    }
  }
  else if(n<H){
    for(int i=(l-4>=m ? 0:m-l+4); i<=(m>=4 ? 4:m); ++i){
      res[n+1][0] = res[n][0];
      res[n+1][1] = res[n][1];
      res[n+1][2] = res[n][2];
      ++res[n+1][i%3];
      p[n+1] = p[n]*c[i][0];
      q[n+1] = q[n]*c[i][1];
      r[n+1] = r[n]*c[i][2];
      s[n+1] = s[n]*d[i][0];
      t[n+1] = t[n]*d[i][1];
      u[n+1] = u[n]*d[i][2];

      dealsum(n+1,m-i,l-4);
    }
  }
  else{
    for(int i=(l-36>=m ? 0:m-l+36); i<=(m>=M ? M:m); ++i){
      res[n+1][0] = res[n][0];
      res[n+1][1] = res[n][1];
      res[n+1][2] = res[n][2];
      ++res[n+1][i%3];
      p[n+1] = p[n]*num[i][0];
      q[n+1] = q[n]*num[i][1];
      r[n+1] = r[n]*num[i][2];
      s[n+1] = s[n]*den[i][0];
      t[n+1] = t[n]*den[i][1];
      u[n+1] = u[n]*den[i][2];

      dealsum(n+1,m-i,l-36);
    }
  }
}

int main()
{
  dealtile(0,M,4*(T-1),1);//num[][]とden[][]に値をセット
  dealsum(0,M,4*T);//パターン数と状態数を数え上げ

  nsp = combin(T,M/2);//七対子のパターン数
  dsp = combin(T,M/2)*pow(6,M/2);//七対子の状態数
  nto = 13;//国士無双のパターン数
  dto = 13*6*pow(4,12);//国士無双の状態数

  std::cout.setf(std::ios::left,std::ios::adjustfield);
  std::cout << std::setw(24) << "Number of Tiles" << std::setw(16) << M << '\n';

  //全てのパターン数と状態数を出力
  std::cout << std::setw(24) << "Total" << std::setw(16) << ntot << std::setw(16) << dtot << '\n';

  //四面子一雀頭形のパターン数と状態数、割合を出力
  std::cout << std::setw(24) << "Legal Hand" << std::setw(16) << nlh << std::setw(16) << dlh << std::setw(16) << 100.0*dlh/dtot << '\n';

  //二盃口のパターン数と状態数、割合を出力
  std::cout << std::setw(24) << "Internal Seven Pairs" << std::setw(16) << ntsis << std::setw(16) << dtsis << std::setw(16) << 100.0*dtsis/dtot << '\n';

  //七対子のパターン数と状態数、割合を出力
  std::cout << std::setw(24) << "Seven Pairs" << std::setw(16) << nsp << std::setw(16) << dsp << std::setw(16) << 100.0*dsp/dtot << '\n';

  //国士無双のパターン数と状態数、割合を出力
  std::cout << std::setw(24) << "Thirteen Orphans" << std::setw(16) << nto << std::setw(16) << dto << std::setw(16) << 100.0*dto/dtot << '\n';
  std::cout.setf(std::ios::scientific);
  std::cout.precision(15);
  //天和の確率を出力
  std::cout << std::setw(24) << "Probability" << std::setw(24) << (dlh+dsp-dtsis+dto)/static_cast<double>(dtot) << '\n';

  return 0;
}

#結果

$ g++ tenhou_prob.cpp
$ ./a.out
Number of Tiles         14
Total                   326520504500    4250305029168216000
Legal Hand              11498658        11353128141498  0.000267113
Internal Seven Pairs    4668            1306741248      3.07446e-08
Seven Pairs             5379616         1505948184576   3.54315e-05
Thirteen Orphans        13              1308622848      3.07889e-08
Probability             3.025448319456385e-06

出力結果の最後の行を見てください。天和の確率が3.025448...e-006と計算されました。これは約33万分の1に等しい確率です。