【麻将】计算普通立直麻将规则下天和的概率

在普通立直麻将规则下,庄家天和的概率是多少?

随便上网一查资料就知道,这是 2143179701279708384171528036000\frac{2143179701279}{708384171528036000},约等于 1330530\frac{1}{330530}

那么怎么计算这个东西呢

假如我们枚举庄家的 1414 张配牌,然后逐一检测是否和牌并计算概率。不考虑赤宝牌,那么不同的配牌有方案 326520504500326520504500 种。这是 3434[0,4][0,4] 内的整数之和等于 1414 的方案数,枚举大于 44 的数的个数容斥,可以得到计算式 t=02(1)tC(34,t)C(34+1415t,341)=326520504500\sum_{t=0}^{2}(-1)^tC(34,t)C(34+14-1-5t,34-1)=326520504500。这样枚举量过于巨大,在普通的计算机上运行起码需要数小时的时间。

换一种思路,枚举和牌的牌型。国士无双型和牌型有 1313 种,七对子型有 C(34,7)=5379616C(34,7)=5379616 种,一般型的计算复杂一些,我们可以大致估计一下:3434 种雀头,2121 种顺子,3434 种刻子,我们可以依次枚举雀头和四组面子。面子之间不分顺序,可以规定一个顺序并从小到大枚举,过程中筛去不合法的方案和重复的方案,这样枚举量大约在 34×(34+21)44!1.3×10734\times\frac{(34+21)^4}{4!}\approx 1.3\times 10^7 这个数量级,实现得当的程序运行数秒钟即可得出结果。实际上一般型有 1149865811498658 种和牌型,其中 46684668 种同时也符合七对子型,避免重复计算需要减去。

326520504500326520504500 种牌型里,每种牌型并不是等概率的。概率可以这样计算:假设所有牌都不一样,那配牌有 C(136,14)=4250305029168216000C(136,14)=4250305029168216000 种。如果一种牌型里一种牌出现了 kk 张,那它的方案数就要乘以 C(4,k)C(4,k)。例如,每一种七对子型出现的概率都是 C(4,2)7C(136,14)=2799364250305029168216000\frac{C(4,2)^7}{C(136,14)}=\frac{279936}{4250305029168216000}

剩下的就是写程序了,其实还是没那么容易写的,代码如下:

#include <cstdio>
#include <unordered_set>
#include <algorithm>
typedef long long ll;
const int c4[5] = {1,4,6,4,1};
int p[4]; ll g[4];
// p=牌型数 c=概率乘以C(136,14)
// 0=国士型 1=七对子型 2=一般型 3=一般型内的七对子型(负值)
int a[34]; // 0~26三色数牌,27~33字牌
ll hsh = 0, h[34];
std::unordered_set<ll> T; // 哈希判重
void find(int d, int x) {
	if (!x) {
		if (T.count(hsh)) return;
		T.insert(hsh);
		ll ng = 1;
		bool ct = true;
		for (int i = 0; i < 34; i++) {
			ct &= (a[i] == 0 || a[i] == 2);
			ng *= c4[a[i]];
		}
		p[2]++;
		g[2] += ng;
		if (ct) {
			p[3]--;
			g[3] -= ng;
		}
		return;
	}
	// 34~54三色顺子
	for (int i = d; i < 55; i++) {
		if (i < 34) {
			if (a[i] <= 1) {
				a[i] += 3;
				hsh += 3 * h[i];
				find(i+1, x-1);
				hsh -= 3 * h[i];
				a[i] -= 3;
			}
		} else {
			int y = (i-34) + (i-34)/7 * 2;
			if (a[y] <= 3 && a[y+1] <= 3 && a[y+2] <= 3) {
				a[y]++; a[y+1]++; a[y+2]++;
				hsh += 31 * h[y];
				find(i, x-1);
				hsh -= 31 * h[y];
				a[y]--; a[y+1]--; a[y+2]--;
			}
		}
	}
}
int main()  {
	p[0] = 13;
	g[0] = p[0] * 100663296ll; // C(4,1)^12 * C(4,2)
	p[1] = 5379616; // C(34,7)
	g[1] = p[1] * 279936ll; // C(4,2)^7
	h[0] = 1;
	for (int i = 1; i < 34; i++)
		h[i] = h[i-1] * 5;
	for (int i = 0; i < 34; i++) {
		a[i] += 2;
		hsh += 2 * h[i];
		find(0, 4);
		hsh -= 2 * h[i];
		a[i] -= 2;
	}
	int ap = 0; ll ag = 0;
	for (int i = 0; i < 4; i++) {
		printf("%d %lld\n", p[i], g[i]);
		ap += p[i]; ag += g[i];
	}
	printf("%d %lld\n", ap, ag);
	ll q = 4250305029168216000; // C(136,14);
	ll gcd = std::__gcd(ag, q);
	ag /= gcd; q /= gcd;
	printf("%lld/%lld\n", ag, q);
}

程序输出:

13 1308622848
5379616 1505948184576
11498658 11353128141498
-4668 -1306741248
16873619 12859078207674
2143179701279/708384171528036000

可以看到输出和上面给出的数据都是一致的,那个 2143179701279708384171528036000\frac{2143179701279}{708384171528036000} 就是最终结果了