[JZOJ4913]告别

题目大意

给定一个$1$至$n$的排列$A$。有$m$次操作，每次随机选择排列中的一个有序三元组轮换，求$m$次操作之内（包括$m$次）将其变成排列$B$的概率。
结果对$998244353$取模。

$1\le n\le 14,1\le m\le 10^9$

---

题目分析

首先可以发现，我们将$A$和$B$同时乘上同一个置换，从前者转移到后者的概率依然是不变的。因此我们考虑将$A$和$B$乘上其逆置换，使得$B$变为$1$到$n$，那么我们研究的就是新的$A$变为单位置换的问题。
其次，我们可以发现同构的置换可以合并。什么意思呢？我们将每一个排列的置换拆成若干个轮换，这个排列就用所有这些轮换的大小来表示。可以发现，虽然这样一种表示方法可以表示多个置换，但是这些置换对同一个三元组进行轮换都可以转移到同一个表示的状态里面。
那么这样会有多少种状态呢？其实就是$n$的划分数$D_n$。在$n\le14$时不超过$150$。
考虑使用矩阵乘法，那么我们需要预处理两两状态之间的转移：对于一个状态，我们可以构造出一个排列使其满足这个置换，然后在这个排列内枚举三元组，计算其转移到的状态，在转移矩阵累加上相应的概率（$\frac{1}{P_n^3}$）。
由于目标状态是单位置换，每个轮换大小都是$1$，并且这样的表示能唯一表达目标状态，因此答案不会算重。
注意到达了目标状态就不能再转移了，因此我们要特殊处理转移矩阵中从目标状态出发的状态，让其只能以100%的概率转移到自己。
为了方便，我们可以先计算方案总数再乘上$\left(P_n^3\right)^m$的逆元。
至于状态表示我们用哈希或map。
时间复杂度$\mathrm O(D_nn^3\log_2n+D_n^3\log_2m)$。

---

代码实现

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <vector>
#include <map>

using namespace std;

typedef vector<int> V;

const int P=998244353;
const int N=20;
const int S=150;

struct matrix
{
    int num[S][S];
    int r,c;
}TRS,zero,f;

inline matrix operator*(matrix a,matrix b)
{
    matrix ret;
    memset(ret.num,0,sizeof ret.num);
    ret.r=a.r,ret.c=b.c;
    for (int i=0;i<ret.r;i++)
        for (int j=0;j<ret.c;j++)
            for (int k=0;k<a.c;k++)
                (ret.num[i][j]+=1ll*a.num[i][k]*b.num[k][j]%P)%=P;
    return ret;
}

inline matrix operator^(matrix x,int y)
{
    matrix ret=zero;
    for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;
    return ret;
}

int A[N],B[N],trs[N];
map<V,int> id;
V state[S];
int n,m,cnt;

int quick_power(int x,int y)
{
    int ret=1;
    for (;y;y>>=1,x=1ll*x*x%P) if (y&1) ret=1ll*ret*x%P;
    return ret;
}

V tmp,nw;

void dfs(int lst,int sum)
{
    if (sum==n)
    {
        state[++cnt]=tmp,id[tmp]=cnt;
        return;
    }
    for (int i=lst;i<=n-sum;i++)
    {
        tmp.push_back(i);
        dfs(i,sum+i);
        tmp.pop_back();
    }
}

bool vis[N];
int con[N],tot[N];

void makesta()
{
    tot[0]=0;
    for (int i=1;i<=n;i++) vis[i]=0;
    for (int i=1;i<=n;i++)
    {
        if (vis[i]) continue;
        tot[++tot[0]]=0;
        for (int x=i;!vis[x];x=con[x]) vis[x]=1,tot[tot[0]]++;
    }
    sort(tot+1,tot+1+tot[0]);
    nw.clear();
    for (int i=1;i<=tot[0];i++) nw.push_back(tot[i]);
}

void trans()
{
    TRS.r=TRS.c=cnt;
    for (int x=1;x<=cnt;x++)
    {
        tmp=state[x];
        int l=0,r;
        for (vector<int>::iterator it=tmp.begin();it!=tmp.end();it++,l=r)
        {
            r=l+*it;
            for (int i=l+1;i<=r;i++) con[i]=i!=r?i+1:l+1;
        }
        for (int i=1;i<=n;i++)
            for (int j=1;j<=n;j++)
                if (i!=j)
                    for (int k=1;k<=n;k++)
                        if (i!=k&&k!=j)
                        {
                            swap(con[i],con[j]),swap(con[i],con[k]);
                            makesta();
                            TRS.num[x-1][id[nw]-1]++;
                            swap(con[i],con[k]),swap(con[i],con[j]);
                        }
    }
    memset(TRS.num[0],0,sizeof TRS.num[0]);
    TRS.num[0][0]=n*(n-1)*(n-2);
    zero.r=zero.c=cnt;
    for (int i=0;i<cnt;i++) zero.num[i][i]=1;
    memcpy(con+1,A+1,n*sizeof(int));
    makesta();
    f.r=1,f.c=cnt,f.num[0][id[nw]-1]=1;
}

void calc(){f=f*(TRS^m);}

int main()
{
    freopen("goodbye.in","r",stdin),freopen("goodbye.out","w",stdout);
    scanf("%d%d",&n,&m);
    for (int i=1;i<=n;i++) scanf("%d",&A[i]);
    for (int i=1;i<=n;i++) scanf("%d",&B[i]),trs[B[i]]=i;
    for (int i=1;i<=n;i++) A[i]=trs[A[i]];
    dfs(1,0),trans(),calc();
    printf("%d\n",1ll*f.num[0][0]*quick_power(quick_power(n*(n-1)*(n-2),m),P-2)%P);
    fclose(stdin),fclose(stdout);
    return 0;
}