一道清真好题，只是之前没有接触过，才感觉非常困难。
感谢yn大佬。

我的做法是dp套dp的做法。
首先如果你有一个麻将自动机，
这个自动机有若干节点，以及节点之间的转移边，以及一个起始状态节点，和一个终止状态集合。

终止状态集合中的节点代表胡了。
你可以给这个自动机新添加一种花色，麻将的数量在0~4之间，然后沿着转移边走到下一个状态。

如果你有这样一个自动机，那么考虑一个dp，\(f[i][j][k]\)表示考虑了前\(i\)种花色，选了\(j\)张麻将，在自动机的\(k\)节点时的无序集合数。
然后一个一定长度的无序集合，映射到相同个数的排列数。
就可以愉快的算出答案了。

当然期望按照套路转化成还未结束的概率和。

丑陋的代码

#include <bits/stdc++.h>using namespace std;typedef long long i64;
const int N=105;
const int MM=4000;
const int M=998244353;namespace{void Addt(int &x,int y){(x+=y)>=M?x-=M:0;}int mul(int x,int y){return (i64)x*y%M;}int fp(int x,int y){int ret=1;for (; y; y>>=1,x=mul(x,x))if (y&1) ret=mul(ret,x);return ret;}
}
struct st{int a[3][3];st(){}st(int x){for (int i=0; i<3; ++i)for (int j=0; j<3; ++j)a[i][j]=x;}st operator +(const st &_) const{st ret;for (int i=0; i<3; ++i)for (int j=0; j<3; ++j)ret.a[i][j]=max(a[i][j],_.a[i][j]);return ret;}st operator +(const int &num) const{st ret(-1);for (int i=0; i<3; ++i)for (int j=0; j<3; ++j)for (int k=0; k<3&&i+j+k<=num; ++k){if (a[i][j]==-1) continue;ret.a[j][k]=min(max(ret.a[j][k],a[i][j]+(num-k-i-j)/3+i),4);}return ret;}bool operator <(const st &_) const{for (int i=0; i<3; ++i)for (int j=0; j<3; ++j)if (a[i][j]!=_.a[i][j]) return a[i][j]<_.a[i][j];return 0;}bool operator !=(const st &_) const{for (int i=0; i<3; ++i)for (int j=0; j<3; ++j)if (a[i][j]!=_.a[i][j]) return 1;return 0;}
};
map<st,int> mst;
struct mj{st f,g;int cnt;//number of different dui zibool operator <(const mj &_) const{if (cnt!=_.cnt) return cnt<_.cnt;if (f!=_.f) return f<_.f;return g<_.g;}mj operator +(const int &a) const{mj ret=*this;if (a>=2){ret.cnt=min(ret.cnt+1,7);ret.g=(ret.g+a)+(ret.f+(a-2));}else ret.g=ret.g+a;ret.f=ret.f+a;return ret;}bool ri(){if (cnt>=7) return 1;//carereturn g.a[0][0]>=4;}
}m_j[MM];
map<mj,int> mmj;
int tots,totm;
void Dfscomplicated(const st &&s){if (mst.count(s)) return;mst[s]=++tots;for (int i=0; i<=4; ++i){Dfscomplicated(s+i);}
}
void Dfshdaewr(const mj &&m){//getchar();if (mmj.count(m)) return;//cerr<<"Dfshdaewr"<<" "<<totm<<" "<<m.cnt<<" "<<mst[m.f]<<" "<<mst[m.g]<<" "<<(mmj.find(m)==mmj.end())<<endl;mmj[m]=++totm;//cerr<<mmj[m]<<endl;m_j[totm]=m;for (int i=0; i<=4; ++i){Dfshdaewr(m+i);}
}
st zst(){st ret=st(-1);ret.a[0][0]=0;return ret;
}
mj zmj(){mj ret;ret.f=zst();ret.g=st(-1);ret.cnt=0;return ret;
}
int binomial[5][5],trans[MM][5];
void Prework(){for (int i=0; i<=4; ++i){binomial[i][0]=1;for (int j=1; j<=i; ++j)Addt(binomial[i][j]=binomial[i-1][j-1],binomial[i-1][j]);}//cerr<<"???"<<endl;Dfscomplicated(zst());//cerr<<"????"<<tots<<endl;Dfshdaewr(zmj());//cerr<<"?????"<<totm<<endl;for (int i=1; i<=totm; ++i)for (int j=0; j<=4; ++j){trans[i][j]=mmj[m_j[i]+j];//cerr<<i<<" "<<j<<" "<<trans[i][j]<<endl;}//cerr<<"Preworkend"<<endl;
}
int n,usd[N],f[N][N*4][MM];
int main(){cin>>n;for (int i=1; i<=13; ++i){int x,y;cin>>x>>y;++usd[x];}Prework();//cerr<<"!!!!!"<<binomial[0][0]<<endl;f[0][0][1]=1;for (int i=0; i<n; ++i)for (int j=0; j<=i*4; ++j)for (int d=1; d<=totm; ++d)if (f[i][j][d]){//cerr<<i<<" "<<j<<" "<<d<<" "<<f[i][j][d]<<" "<<m_j[d].ri()<<endl;//getchar();int tmp=f[i][j][d];for (int k=usd[i+1]; k<=4; ++k){//cerr<<"trans:"<<k<<" "<<trans[d][k]<<endl;Addt(f[i+1][j+k][trans[d][k]],mul(tmp,binomial[4-usd[i+1]][k-usd[i+1]]));}}//cerr<<"????"<<endl;int ans=0;for (int j=13; j<=n*4; ++j){int sum=0;int nend=0;for (int i=1; i<=totm; ++i){Addt(sum,f[n][j][i]);//if (f[n][j][i]) cerr<<j<<" "<<i<<" "<<f[n][j][i]<<" "<<ans<<" "<<sum<<" "<<nend<<endl;if (!m_j[i].ri()) Addt(nend,f[n][j][i]);}Addt(ans,mul(nend,fp(sum,M-2)));}cout<<ans;
}