这题好牛逼啊。。。

虽然说也是 DP 套 DP，但是感觉比 TJOI 那题高明到哪里去了（

我们先考虑如何计算期望。如果设 \(f_i\) 为拿到 \(i\) 张牌后胡的方案数，这个并不是很好做，因为你要考虑前面不胡。

所以我们设 \(f_i\) 为手上有 \(13+i\) 张牌但还没有胡的方案数，那么答案就是：

\[\frac{\sum_{i=1}^{4n-13}f_ii!(4n-13-i)!}{(4n-14)!}+1
\]

发现一个手上牌的顺序不影响是否胡，所以有两个阶乘用来打乱排序（因为牌与牌直接互不相同）。

我们考虑外层的 DP 建立一个“胡牌自动机”，所以需要写一个用来判定是否胡牌的 DP。

加入新的牌就是在自动机上面走一步。我们设 \(dp_{i,x,y,0/1}\) 为加入了前 \(i\) 类牌，现在手上有 \(x\) 对 \((i-1,i)\) 的牌和 \(y\) 张 \(i\)，\(0/1\) 表示当前无/有对子。表示的是在这个状态下面子最多个数。

我们将发现如果 \(x,y\geq 3\)，那么可以自动对出来一个面子，所以有 \(0\leq x,y<3\)。

我们把这两个维度写到矩阵上面。然后交换维度变成 \(dp_{i,0/1,x,y}\)。

考虑加入 \(x\) 张牌时如何从 \(dp_{i-1}\) 转移到 \(dp_i\)。

考虑加入 \(x\) 张牌的时候，先凑到 \((i-2,i-1)\) 里面去，然后自己和自己凑面子，最后加上已经有的面子个数。

如果 \(x\geq 2\)，那么我们可以砍掉一对对子再加进去。

也就是：\((f_{i-1,0},x)\to f_{i,0},(f_{i-1,1},x)\to f_{i,1},(f_{i-1,0},x-2)\to f_{i,1}\)。

然后我们再额外记录对子个数即可。

然后我们用一个 std::map 离散化一下就可以建出“胡牌自动机”了。

这个是主要的，在这个 DFA 上面 DP 反而很简单。

直接考虑当前加 \(x\) 张牌，然后从牌堆里面把这 \(x\) 张牌选出来，丢到下一个 DFA 节点里面去就可以了。

我们可以把所有胡掉的节点压成一个节点，然后标记 \(t=-1\)，可以节省空间。

#include<cstdio>
#include<map>
typedef int ui;
const ui M=105,mod=998244353;
const ui C[5][5]={{1,0,0,0,0},{1,1,0,0,0},{1,2,1,0,0},{1,3,3,1,0},{1,4,6,4,1}};
ui n,m,a[M],fac[M*4],ifac[M*4],dp[2][M*4][2100];
inline ui min(const ui&a,const ui&b){
    return a>b?b:a;
}
inline void max(ui&a,const ui&b){
    if(b>a)a=b;
}
struct HuAutomation{
    ui tot;
    struct Matrix{
        #define For for(ui i=0;i^3;++i)for(ui j=0;j^3;++j)
        ui p[3][3];
        Matrix(){
            For p[i][j]=-1;
        }
        inline ui*operator[](const ui&x){
            return p[x];
        }
        inline bool operator==(Matrix q)const{
            For if(p[i][j]!=q[i][j])return false;return true;
        }
        inline bool operator<(Matrix q)const{
            For if(p[i][j]!=q[i][j])return p[i][j]<q[i][j];return false;
        }
        inline bool CheckHu()const{
            For if(~p[i][j]&&p[i][j]>3)return true;return false;
        }
        inline void Update(Matrix bel,const ui&x){
            For if(~bel[i][j])for(ui k=0;k^3&&i+j+k<=x;++k){
                max(p[j][k],min(i+bel[i][j]+(x-i-j-k)/3,4));
            }
        }
        #undef For
    };
    struct Node{
        ui t;Matrix P[2];
        inline Node(){
            t=0;P[0]=P[1]=Matrix();
        }
        inline bool operator<(Node it)const{
            return t==it.t?P[0]==it.P[0]?P[1]==it.P[1]?0:P[1]<it.P[1]:P[0]<it.P[0]:t<it.t;
        }
        inline Node Hu()const{
            Node x;x.t=-1;return x;
        }
        inline bool CheckHu()const{
            return!~t||t>=7||P[1].CheckHu();
        }
        inline Node operator+(const ui&x){
            if(CheckHu())return Hu();Node ans;
            ans.P[0].Update(P[0],x);ans.P[1].Update(P[1],x);ans.t=t;
            if(x>=2)ans.P[1].Update(P[0],x-2),++ans.t;
            if(ans.CheckHu())ans=Hu();
            return ans;
        }
    }tr[2100];
    ui trans[2100][5];
    inline ui*operator[](const ui&x){
        return trans[x];
    }
    inline Node Hu()const{
        return tr[0].Hu();
    }
    std::map<Node,ui>P;
    #define pos(x) (P.count(x)?P[x]:(tr[P[x]=++tot]=x,tot))
    inline void Build(){
        tr[1].P[0][0][0]=0;P[tr[1]]=++tot;P[tr[2]=Hu()]=++tot;
        for(ui i=1;i<=tot;++i)if(i!=2){
            for(ui t=0;t<=4;++t)trans[i][t]=pos(tr[i]+t);
        }
    }
}Hu;
signed main(){
    ui ans(0);
    Hu.Build();scanf("%u",&n);m=(n<<2)-13;dp[0][0][1]=1;
    fac[0]=fac[1]=ifac[0]=ifac[1]=1;
    for(ui i=2;i<=m;++i)fac[i]=1ull*fac[i-1]*i%mod,ifac[i]=1ull*(mod-mod/i)*ifac[mod%i]%mod;
    for(ui i=2;i<=m;++i)ifac[i]=1ull*ifac[i-1]*ifac[i]%mod;
    for(ui i=1;i<=13;++i){
        ui w,t;scanf("%u%u",&w,&t);++a[w];
    }
    for(ui i=0;i<n;++i){
        const ui&q=i&1,p=q^1;
        for(ui j=0;j<=m;++j)for(ui S=1;S<=Hu.tot;++S)if(dp[q][j][S]){
            for(ui x=0;x<=4-a[i+1];++x){
                dp[p][j+x][Hu[S][a[i+1]+x]]=(dp[p][j+x][Hu[S][a[i+1]+x]]+1ull*C[4-a[i+1]][x]*dp[q][j][S])%mod;
            }
        }
        for(ui i=0;i<=m;++i)for(ui S=1;S<=Hu.tot;++S)dp[q][i][S]=0;
    }
    for(ui i=1;i<=m;++i)for(ui S=1;S<=Hu.tot;++S)if(S!=2){
        ans=(ans+1ull*fac[i]*fac[m-i]%mod*dp[n&1][i][S])%mod;
    }
    printf("%u",1ull*ans*ifac[m]%mod+1);
}