这是我第一次写$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 \times 3$的矩阵,其转移本质上是接收一个数字(表示下一种花色的牌的数量),然后变成了另一个$3 \times 3$的矩阵,我们可以把它看成一个自动机的模型。事实上这样的$dp$状态($3 \times 3$的矩阵)并不很多。我们用结构体$State$来形容这个矩阵。为了表示整个当前胡牌状态,我们需要两个$State$分别表示有没有选对子(即$dp_0$和$dp_1$),以及$cnt$。我们把这三个东西放到结构体$Mahjong$中,我们很容易可以得到$Mahjong$的$Trans$关系。经搜索$Mahjong$的状态数总共有$3956$个。
回过来看我们要算的答案:$ans = \sum\limits_{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}}} \binom{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{\sum\limits_{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 + ;
const int M = 4e3 + ;
const int MOD = ; namespace {
int ch[M][], fac[];
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 == )? : Mul(MOD - MOD / x, Inv(MOD % x)); }
void Math_init() {
for (int i = ; i < M; ++i) {
ch[i][] = ;
for (int j = ; j <= min(i, ); ++j) {
ch[i][j] = Add(ch[i - ][j - ], ch[i - ][j]);
}
}
fac[] = ;
for (int i = ; i < ; ++i) {
fac[i] = Mul(fac[i - ], i);
}
}
} bool Chkmax(int &a, int b) {
return (a < b)? (a = b, ) : ();
} struct State {
int dp[][]; // (last last bar, last bar)
State() {
memset(dp, -, sizeof dp);
} friend bool operator < (State a, State b) {
for (int i = ; i < ; ++i)
for (int j = ; j < ; ++j)
if (a.dp[i][j] != b.dp[i][j])
return a.dp[i][j] < b.dp[i][j];
return ;
} friend State Max(State a, State b) {
for (int i = ; i < ; ++i)
for (int j = ; j < ; ++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 = ; i < ; ++i)
for (int j = ; j < ; ++j)
if (~a.dp[i][j])
for (int k = ; k < && i + j + k <= b; ++k)
Chkmax(c.dp[j][k], min(i + a.dp[i][j] + (b - i - j - k) / , ));
return c;
} }; struct Mahjong {
pair<State, State> god;
int cnt;
Mahjong() {
memset(god.first.dp, -, sizeof god.first.dp);
memset(god.second.dp, -, sizeof god.second.dp);
god.first.dp[][] = cnt = ;
} 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 >= ), );
a.god.second = Trans(a.god.second, b);
if (b >= ) {
a.god.second = Max(a.god.second, Trans(a.god.first, b - ));
}
a.god.first = Trans(a.god.first, b);
return a;
} bool right() {
if (cnt == ) return ;
for (int i = ; i < ; ++i)
for (int j = ; j < ; ++j)
if (god.second.dp[i][j] == ) return ;
return ;
} } mahjong[M]; int n, tot;
map<Mahjong, int> idx;
bool win[M];
int org[N], f[N][M][ * N], trans[M][]; 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 = ; i <= ; ++i) {
Dfs_mahjong(Trans(now, i));
}
} int main() {
Math_init();
Dfs_mahjong(Mahjong()); for (int i = ; i <= tot; ++i) {
for (int j = ; j <= ; ++j) {
trans[i][j] = idx[Trans(mahjong[i], j)];
}
} scanf("%d", &n);
for (int i = , x; i < ; ++i) {
scanf("%d%*d", &x);
++org[x];
} f[][][] = ;
for (int i = , cp = ; i < n; ++i) { // consider 1 ... i
cp += org[i + ];
for (int j = ; j <= tot; ++j) { // mahjong j
for (int l = org[i + ]; l <= ; ++l) { // trans
int *nf = f[i + ][trans[j][l]], *ff = f[i][j];
int tmp = Mul(ch[ - org[i + ]][l - org[i + ]], fac[l - org[i + ]]);
for (int k = ; k + l <= * 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 + ]], tmp)));
}
}
}
} int ans = , dw = ;
for (int i = ; i <= * n; ++i) {
int up = ;
for (int j = ; j <= tot; ++j) {
if (!win[j]) Upd(up, f[n][j][i]);
}
Upd(ans, Mul(up, Inv(dw)));
dw = Mul(dw, * n - i);
}
printf("%d\n", ans); return ;
}