ZJOI2019神题，间接送我退役的神题233

考场上由于T2写挂去写爆搜的时候已经没多少时间了，所以就写挂了233

这里不多废话直接开始讲正解吧，我们把算法分成两部分

1.建一个“胡牌自动机”

首先我们发现这题不能转化为一般DP问题求解的最大瓶颈就是因为它的状态很诡异

但是我们细细一想，形如\(\{1,1,1,2,3,4,5,6,7,8,9,10,12,12\}\)和\(\{3,3,3,5,6,7,9,10,11,14,15,16,20,20\}\)的本质其实是一样的（都是三个顺子+一个刻子+一个雀头的形式，且相对大小分布相同）

那么考虑到一共只有\(13\)张牌，那么最后的有效集合数必然很少，如果我们可以构造出这些集合那么可以帮助DP了

如果你知道DP套DP的姿势，那么就会套路地明白内层的DP应该用自动机的转移形式，所以我们大致的目标就出现了：

建出一个自动机，其中每一个节点都是一个状态，然后每个节点向子节点连的边就是摸牌之后能转移到的状态

所以我们就先考虑怎么建这个自动机，还是考虑DP

我们设\(f_{0/1,i,j,k}\)表示目前处理完了前\(i\)种牌，还剩下\(j\)组\((i-1,i)\)以及\(k\)张\(i\)，且是否（用\(0/1\)表示）存在雀头时最多的面子数

由于当\(j,k\ge 3\)时可以直接用\(i-1,i\)组成刻子，因此我们状态中的\(j,k\le 2\)

那么我们就开一个\(3\times3\)的矩阵表示状态，单个转移的时候枚举用于拼面子，与\((i-1,i)\)拼顺子，以及用来保留的个数，然后将多余的拿来拼刻子即可

然后考虑每次多出一个数值，我们直接讨论是否要留出一个雀头，分情况转移即可

我们发现这样并没有考虑七对子的情况，这个没关系，我们直接记录雀头的个数，特判了即可

因此剩下终止状态的判断就很简单了，直接把能胡的点作为终止节点结束即可

具体构造的过程可以用一个map来去重，然后用BFS来扩展状态，这个具体看代码

2.期望DP

先说一句，前面由于自动机对于任意数据构造相同，所以总点数是固定的\(2092\)，如前言所述不大

那么考虑DP求解最后的问题，我们用一个经典套路，将求\(\ge\)的化为求\(>\)然后最后加上等于的情况即可

具体到这道题上，其实就是求出\(i\)轮后不胡的方案数然后乘上贡献，最后的总答案加\(1\)即可

那么大致的DP方程就有了，我们设\(f_{i,j,k}\)表示选了前\(i\)种牌，一共用了\(j\)张，此时位于自动机上\(k\)号点的方案数

那么转移其实很简单，枚举第\(i\)种牌选的张数\(t\)，\(f_{i,j,k}\)向\(f_{i+1,j+t,node_k.son_{a_i+t}}\)转移，注意不要忘记乘上组合数\(C_{4-a_i}^t\)

求出\(f\)数组后最后的答案统计就十分简单了，令\(m=4n-13\)：

这个很好理解，因为前后都可以随便选，因此不再赘述

综上，我们得到了一个复杂度上界为\(2092\cdot n^2\)的优秀做法，因为常数很小且跑不满因此足以通过

CODE

#include<cstdio>

#include<cstring>

#include<map>

#define RI register int

#define CI const int&

using namespace std;

const int N=105,MX=2100,mod=998244353;

inline void maxer(int& x,CI y)

{

    if (y>x) x=y;

}

inline int min(CI a,CI b)

{

    return a<b?a:b;

}

struct Matrix

{

    int mat[3][3];

    inline int* operator [] (CI x) { return mat[x]; }

    inline Matrix(void)

    {

        memset(mat,-1,sizeof(mat));

    }

    friend inline bool operator != (Matrix A,Matrix B)

    {

        for (RI i=0;i<3;++i) for (RI j=0;j<3;++j)

        if (A[i][j]!=B[i][j]) return 1; return 0;

    }

    friend inline bool operator < (Matrix A,Matrix B)

    {

        for (RI i=0;i<3;++i) for (RI j=0;j<3;++j)

        if (A[i][j]!=B[i][j]) return A[i][j]<B[i][j];

    }

    inline bool Com_Hu(void)

    {

        for (RI i=0;i<3;++i) for (RI j=0;j<3;++j)

        if (mat[i][j]>=4) return 1; return 0;

    }

    inline void flush(Matrix pre,CI num)

    {

        for (RI i=0;i<3;++i) for (RI j=0;j<3;++j)

        if (~pre[i][j]) for (RI k=0;k<3&&i+j+k<=num;++k)

        maxer(mat[j][k],min(i+pre[i][j]+(num-i-j-k)/3,4));

    }

};

struct Hu_Auto_Node

{

    Matrix p[2]; int cur,ch[5];

    inline Hu_Auto_Node(void)

    {

        memset(ch,0,sizeof(ch)); cur=0; p[0]=p[1]=Matrix();

    }

    inline bool is_Hu(void)

    {

        if (cur>=7) return 1; return p[1].Com_Hu();

    }

    inline void Hu(void)

    {

        memset(ch,0,sizeof(ch)); cur=-1; p[0]=p[1]=Matrix();

    }

    friend inline bool operator < (Hu_Auto_Node A,Hu_Auto_Node B)

    {

        if (A.cur!=B.cur) return A.cur<B.cur;

        if (A.p[0]!=B.p[0]) return A.p[0]<B.p[0];

        if (A.p[1]!=B.p[1]) return A.p[1]<B.p[1]; return 0;

    }

    friend inline Hu_Auto_Node operator + (Hu_Auto_Node A,CI num)

    {

        if (A.is_Hu()) return A.Hu(),A; Hu_Auto_Node s;

        s.p[0].flush(A.p[0],num); s.p[1].flush(A.p[1],num);

        if (num>=2) s.p[1].flush(A.p[0],num-2);

        s.cur=A.cur+(num>=2); if (s.is_Hu()) s.Hu(); return s;

    }

};

class Hu_Automation

{

    private:

        map <Hu_Auto_Node,int> Hash;

        inline void expand(CI id)

        {

            for (RI i=0;i<=4;++i)

            {

                Hu_Auto_Node son=node[id]+i;

                if (!Hash.count(son)) node[Hash[son]=++tot]=son;

                node[id].ch[i]=Hash[son];

            }

        }

    public:

        Hu_Auto_Node node[2100]; int tot;

        inline Hu_Automation(void)

        {

            node[1].p[0][0][0]=0; node[tot=2].cur=-1;

            Hash[node[1]]=1; Hash[node[2]]=2; expand(1);

            for (RI i=3;i<=tot;++i) expand(i);

        }

        /*inline void check(void)

        {

            printf("%d\n",tot); for (RI i=1;i<=tot;++i,putchar('\n'))

            for (RI j=0;j<=4;++j) printf("%d ",node[i].ch[j]);

        }*/

}HA;

int f[2][N<<2][MX],a[N],fact[N<<2],inv[N<<2],n,m,x,y,nw,ans;

inline void inc(int& x,CI y)

{

    if ((x+=y)>=mod) x-=mod;

}

inline int quick_pow(int x,int p=mod-2,int mul=1)

{

    for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;

}

inline void init(CI n)

{

    RI i; for (fact[0]=i=1;i<=n;++i) fact[i]=1LL*fact[i-1]*i%mod;

    for (inv[n]=quick_pow(fact[n]),i=n-1;~i;--i) inv[i]=1LL*inv[i+1]*(i+1)%mod;

}

inline int C(CI n,CI m)

{

    return 1LL*fact[n]*inv[m]%mod*inv[n-m]%mod;

}

inline void calc(CI x,CI y)

{

    inc(ans,1LL*f[nw][x][y]*fact[x]%mod*fact[m-x]%mod);

}

int main()

{

    RI i,j,k,t; for (scanf("%d",&n),i=1;i<=13;++i) scanf("%d%d",&x,&y),++a[x];

    for (init(m=(n<<2)-13),f[0][0][1]=i=1;i<=n;++i)

    for (nw=i&1,memset(f[nw],0,sizeof(f[nw])),j=m;~j;--j)

    for (k=1;k<=HA.tot;++k) if (f[nw^1][j][k]) for (t=0;t<=4-a[i];++t)

    inc(f[nw][j+t][HA.node[k].ch[a[i]+t]],1LL*f[nw^1][j][k]*C(4-a[i],t)%mod);

    for (nw=n&1,i=1;i<=m;++i) for (calc(i,1),j=3;j<=HA.tot;++j) calc(i,j);

    return printf("%d",(1LL*ans*inv[m]%mod+1)%mod),0;

}