这是我第一次写$dp ; of ; dp$，大致思路参考了xyx的做法，可能有些地方不太一样，但也许会详细一点。
考虑给你一副牌，如何判断这副牌是否是胡的。
容易发现相同的顺子不会选三个以上，于是考虑花色从小到大进行$dp$。令$dp_{0/1, i, j, k}$表示是否有对子，考虑了前$i$种花色的牌，以花色$i - 1$的牌为开头的顺子准备选$j$个，以花色$i$的牌为开头的顺子准备选$k$个，此时能选的最大的面子数。转移只要枚举以花色$i + 1$的牌开头的顺子准备选几个，剩下的牌组成刻子就行了(具体细节可以看代码中结构体$State$中的$Trans$函数)。
判定胡牌的条件有两种：$exists (j,k), dp_{1, n, j, k} ge 4$；记$cnt$为牌数$ge 2$的花色数，$cnt ge 7$。
可以发现这个$dp$和状态$i$这一维关系不是很大，$dp_{0}$和$dp_{1}$的转移方式也都是一样的。只考虑有$dp_{*,j,k}$这个$3 imes 3$的矩阵，其转移本质上是接收一个数字（表示下一种花色的牌的数量），然后变成了另一个$3 imes 3$的矩阵，我们可以把它看成一个自动机的模型。事实上这样的$dp$状态（$3 imes 3$的矩阵）并不很多。我们用结构体$State$来形容这个矩阵。为了表示整个当前胡牌状态，我们需要两个$State$分别表示有没有选对子（即$dp_0$和$dp_1$），以及$cnt$。我们把这三个东西放到结构体$Mahjong$中，我们很容易可以得到$Mahjong$的$Trans$关系。经搜索$Mahjong$的状态数总共有$3956$个。
回过来看我们要算的答案：$ans = sumlimits_{a = 13}^{4n} p(a)$，其中$p(a)$表示你总共摸了$a$张牌之后仍然不能胡的概率。于是我们只要算摸了$a$张牌仍然不能胡的排列数即可。令$f_{i, j, k}$表示我们考虑了前$i$种花色的牌，当前胡牌状态为$j$（$j$是一个$Mahjong$类），已经摸了$k$张牌的排列数有多少。转移枚举摸了花色$i + 1$的牌数$z$，则由$f_{i, j, k}$转移到$f_{i + 1, trans(j,z), k + z}$，乘上的系数就是$(4 - org_{i + 1})^{underline{z - org{i + 1}}} inom{k + z - sum_{i + 1}}{z - org_{i + 1}}$。$org_i$表示原有的$13$张牌中花色为$i$的有几张，$sum$则是$org$的前缀和，这个式子的意思就是我们需要在没被选过的$4 - org_{i + 1}$张牌中选$z - org_{i + 1}$张的排列，并且插入到之前的排列中，但前$13$张牌的顺序是固定的。
最后$p(a)$是很好算的，$p(a) = frac{sumlimits_{j ; can ; not ; win } f_{n, j, a}}{ (4n - 13)^{underline{a - 13}} }$。
复杂度为$O(3956 * n^2)$。
#include <bits/stdc++.h> using namespace std; const int N = 1e2 + 5; const int M = 4e3 + 5; const int MOD = 998244353; namespace { int ch[M][5], fac[5]; int Add(int a, int b) { return (a += b) >= MOD? a - MOD : a; } void Upd(int &a, int b) { a = Add(a, b); } int Mul(int a, int b) { return (long long)a * b % MOD; } int Inv(int x) { return (x == 1)? 1 : Mul(MOD - MOD / x, Inv(MOD % x)); } void Math_init() { for (int i = 0; i < M; ++i) { ch[i][0] = 1; for (int j = 1; j <= min(i, 4); ++j) { ch[i][j] = Add(ch[i - 1][j - 1], ch[i - 1][j]); } } fac[0] = 1; for (int i = 1; i < 5; ++i) { fac[i] = Mul(fac[i - 1], i); } } } bool Chkmax(int &a, int b) { return (a < b)? (a = b, 1) : (0); } struct State { int dp[3][3]; // (last last bar, last bar) State() { memset(dp, -1, sizeof dp); } friend bool operator < (State a, State b) { for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) if (a.dp[i][j] != b.dp[i][j]) return a.dp[i][j] < b.dp[i][j]; return 0; } friend State Max(State a, State b) { for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) a.dp[i][j] = max(a.dp[i][j], b.dp[i][j]); return a; } friend State Trans(State a, int b) { State c; for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) if (~a.dp[i][j]) for (int k = 0; k < 3 && i + j + k <= b; ++k) Chkmax(c.dp[j][k], min(i + a.dp[i][j] + (b - i - j - k) / 3, 4)); return c; } }; struct Mahjong { pair<State, State> god; int cnt; Mahjong() { memset(god.first.dp, -1, sizeof god.first.dp); memset(god.second.dp, -1, sizeof god.second.dp); god.first.dp[0][0] = cnt = 0; } friend bool operator < (Mahjong a, Mahjong b) { return a.cnt != b.cnt? a.cnt < b.cnt : a.god < b.god; } friend Mahjong Trans(Mahjong a, int b) { a.cnt = min(a.cnt + (b >= 2), 7); a.god.second = Trans(a.god.second, b); if (b >= 2) { a.god.second = Max(a.god.second, Trans(a.god.first, b - 2)); } a.god.first = Trans(a.god.first, b); return a; } bool right() { if (cnt == 7) return 1; for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) if (god.second.dp[i][j] == 4) return 1; return 0; } } mahjong[M]; int n, tot; map<Mahjong, int> idx; bool win[M]; int org[N], f[N][M][4 * N], trans[M][5]; void Dfs_mahjong(Mahjong now) { if (idx.find(now) != idx.end()) return; mahjong[++tot] = now; win[tot] = now.right(); idx[now] = tot; for (int i = 0; i <= 4; ++i) { Dfs_mahjong(Trans(now, i)); } } int main() { Math_init(); Dfs_mahjong(Mahjong()); for (int i = 1; i <= tot; ++i) { for (int j = 0; j <= 4; ++j) { trans[i][j] = idx[Trans(mahjong[i], j)]; } } scanf("%d", &n); for (int i = 0, x; i < 13; ++i) { scanf("%d%*d", &x); ++org[x]; } f[0][1][0] = 1; for (int i = 0, cp = 0; i < n; ++i) { // consider 1 ... i cp += org[i + 1]; for (int j = 1; j <= tot; ++j) { // mahjong j for (int l = org[i + 1]; l <= 4; ++l) { // trans int *nf = f[i + 1][trans[j][l]], *ff = f[i][j]; int tmp = Mul(ch[4 - org[i + 1]][l - org[i + 1]], fac[l - org[i + 1]]); for (int k = 0; k + l <= 4 * n; ++k) { // have chosen k cards if (!ff[k]) continue; Upd(nf[k + l], Mul(ff[k], Mul(ch[k + l - cp][l - org[i + 1]], tmp))); } } } } int ans = 0, dw = 1; for (int i = 13; i <= 4 * n; ++i) { int up = 0; for (int j = 1; j <= tot; ++j) { if (!win[j]) Upd(up, f[n][j][i]); } Upd(ans, Mul(up, Inv(dw))); dw = Mul(dw, 4 * n - i); } printf("%d ", ans); return 0; }