【GDOI2017模拟12.3】告别

题目大意

给你两个n的排列A,B（A与B不完全相同），你可以进行最多m次操作，每次对于当前的A，你等概率地在1到n选择3个不同的数按顺序排成一个三元组(i,j,k)，然后A[i]—>A[j],A[j]—>A[k],A[k]—>A[i]地交换。求m次操作内把A变成B的概率对998244353取模。

n≤14 m≤$10^8$

优化状态

如果n的阶乘表示状态明显会很大。
考虑到操作是交换，那么可以对应地去改变表示状态的方法。首先给B重新编成1——n，然后原本的i对应的是它在B中的位置（即输入的B[j]=i）。然后又改变了A数组。
之后一个状态可以看成一堆置换，例如：2,3,1，就是一个大小为3的置换。B数组是1到n的，就是n个大小为1的置换。
之后可以枚举每个置换的大小，然后搞出很多个状态来。（n=14时只有135个状态）。然后每个状态，可以还原出多个排列（当然，n个1只会还原出一个状态，所以不用担心会错）。

转移

首先枚举一个非最终状态的状态S，还原出任意一个合法的排列。
然后枚举三元组，交换过后再转回成状态。
考虑到m比较大，要使用矩阵乘法。
存状态可以用map或者hash

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

using namespace std;

const int N=15,M=145,mo=998244353;

typedef long long LL;

int n,m,s,a[M][M],f[M][M],A[N],B[N],id[N],g[N],New[N],now,Aq[N],cnt,b[M][M];

bool bz[N];

map <int,int> H;

int quick(int x,int y)
{
    if (!y) return 1;
    int t=quick(x,y>>1);
    t=(LL)t*t%mo;
    if (y&1) t=(LL)t*x%mo;
    return t;
}

void init(int x,int y,int la)
{
    if (x==0)
    {
        H.insert(make_pair(y,s++));
        return;
    }
    for (int i=la;i<=x;i++) init(x-i,y*(n+1)+i,i);
}

int get(int *a)
{
    cnt=0;
    memset(bz,0,sizeof(bz));
    for (int i=1;i<=n;i++) if (!bz[i])
    {
        Aq[cnt]=0;
        for (int j=i;!bz[j];j=a[j])
        {
            Aq[cnt]++; bz[j]=1;
        }
        cnt++;
    }
    sort(Aq,Aq+cnt);
    int p=0;
    for (int t=0;t<cnt;t++) p=p*(n+1)+Aq[t];
    return p;
}

void link(int x,int y,int la)
{
    if (x==0)
    {
        if (!now)
        {
            now=1; return;
        }
        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 && j!=k)
                {
                    for (int t=1;t<=n;t++) New[t]=g[t];
                    New[k]=g[j]; New[j]=g[i]; New[i]=g[k];
                    a[now][H.find(get(New))->second]++;
                }
        now++;
        return;
    }
    for (int i=la;i<=x;i++)
    {
        g[x]=x-i+1;
        for (int j=x-i+1;j<x;j++) g[j]=j+1;
        link(x-i,y*(n+1)+i,i);
    }
}

void mulf()
{
    memset(b,0,sizeof(b));
    for (int i=0;i<s;i++)
        for (int j=0;j<s;j++)
            for (int k=0;k<s;k++) b[i][j]=(b[i][j]+(LL)f[i][k]*f[k][j])%mo;
    memcpy(f,b,sizeof(f));
}

void mula()
{
    memset(b,0,sizeof(b));
    for (int i=0;i<s;i++)
        for (int j=0;j<s;j++)
            for (int k=0;k<s;k++) b[i][j]=(b[i][j]+(LL)f[i][k]*a[k][j])%mo;
    memcpy(f,b,sizeof(f));
}

void quick(int x)
{
    if (!x) return;
    quick(x>>1);
    mulf();
    if (x&1) mula();
}

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]);
        id[B[i]]=i;
    }
    for (int i=1;i<=n;i++) A[i]=id[A[i]];
    init(n,0,1); link(n,0,1);
    a[0][0]=n*(n-1)*(n-2);
    for (int i=0;i<s;i++) f[i][i]=1;
    quick(m);
    printf("%d\n",(LL)f[H.find(get(A))->second][0]*quick(quick(n*(n-1)*(n-2),m),mo-2)%mo);
    fclose(stdin); fclose(stdout);
    return 0;
}