题意

每次随机一个二元组\((i, j)\)，交换\((a_i, a_j)\)，给定一个\(1 \sim n\)排列\(\{a\}\)，求期望交换多少次后，这个排列会变得有序。

\(n \le 20\)

暴力

一个很显然的暴力是，枚举所有的排列，对每个排列枚举所有的交换方案，然后对排列康托展开一下，列方程组求解：（以下设\(x_0\)代表\(\{1, 2, 3, \cdots n\}\)所对应的期望，\(v\)是\(u\)一次交换后能到达的状态）
\[ \begin{cases} x_0 = 0 \\ x_u = \Big(\dfrac{1}{n(n-1)/2}\sum_v x_v\Big) + 1 \end{cases} \]
复杂度\(\mathcal O((n!)^3)\)，可以过40分

优化

打完暴力发现很多状态对应的期望是一样的，考虑找出它们，猜想可能与需要交换的次数有关。

首先发现每个元素有一个要交换到的位置，如果把元素的位置和最终排序后的位置连一条有向边，会形成若干个简单环，因为每个点入、出度都是\(1\)。

现在需要消掉这些环就能使元素到应该在的位置上，即形成\(n\)个自环。

通过手玩规律，发现2条性质：

1. 交换一个环上的\(2\)个数，则这个环会裂成\(2\)个环，设交换了环上的第\(i, j\)个数，则会裂成\(|i-j|\)和\(n-|i-j|\)两个长度的环。
2. 交换不同环上的\(2\)个数，则两个环会合并。

可以手玩两个\(3\)的环的合并，\(5\)的环的分裂。

状态数明显减少，因为期望值只与各个环的大小有关。可以爆搜出所有的集合，用std::map维护就好，计算出每种状态转移到其他状态的方案数就好。

设\(n\)分解成一些整数的和的方案数为\(D(n)\)，则复杂度\(\mathcal O(D^3(n))\)

DP 得\(D(n) = 627\)，可以过。

Code

貌似是跑的最慢的...

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <set>
#include <cassert>
using namespace std;

const int mod = 1000000007;
inline int add(int x, int y){return x+y>=mod ? x+y-mod : x+y;}
inline int sub(int x, int y){return x-y<0 ? x-y+mod : x-y;}
inline int mul(int x, int y){return 1LL * x * y % mod;}
inline int power(int x, int y){
    int res = 1;
    for(; y; y>>=1, x = mul(x, x)) if(y & 1) res = mul(res, x);
    return res;
}
inline int inv(int x){return power(x, mod - 2);}

const int Len = 22, Unk = 635;

int a[Len], n;
int e[Unk][Unk], res[Unk];
int fac[Len];

map<multiset<int>, int> mp;

bool solveEq(int a[Unk][Unk], int x[Unk], int n){
    for(int i=0; i<n; i++){
        int pos = i;
        for(int j=i+1; !a[pos][i] && j < n; j++) pos = j;
        if(!a[i][pos]) return false;
        if(pos != i) swap_ranges(a[i], a[i]+n+1, a[pos]);
        for(int j=i+1; j<n; j++){
            const int rate = mul(a[j][i], inv(a[i][i]));
            for(int k=0; k<=n; k++)
                a[j][k] = sub(a[j][k], mul(a[i][k], rate));
        }
    }
    for(int i=n-1; i>=0; i--){
        for(int j=i+1; j<n; j++)
            a[i][n] = sub(a[i][n], mul(x[j], a[i][j]));
        x[i] = mul(a[i][n], inv(a[i][i]));
    }
    return true;
}

multiset<int> now;
int idx = 0;

void dfs(int nid){
    const multiset<int> back = now;
    for(auto i = back.begin(); i != back.end(); ++i){
        int vi = *i;
        for(auto j = next(i); j != back.end(); ++j){
            int vj = *j;
            now.erase(now.find(vi)); now.erase(now.find(vj));
            now.insert(vi + vj);
            int id;
            if(!mp.count(now)){
                mp[now] = ++idx;
                id = idx;
                dfs(id);
            }
            else id = mp[now];
            e[nid][id] = add(e[nid][id], vi * vj);
            now = back;
        }
        for(int j=1; j*2<=vi; j++){
            now.erase(now.find(vi));
            now.insert(j); now.insert(vi - j);
            int id;
            if(!mp.count(now)){
                mp[now] = ++idx;
                id = idx;
                dfs(id);
            }
            else id = mp[now];
            if(j * 2 != vi) e[nid][id] = add(e[nid][id], vi);
            else e[nid][id] = add(e[nid][id], vi / 2);
            now = back;
        }
    }
}

bool vis[Len];
int cnt = 0;
void getsize(int x){
    ++cnt; vis[x] = true;
    if(vis[a[x]]) return ;
    else getsize(a[x]);
}

int main(){
    scanf("%d", &n);
    for(int i=0; i<n; i++){
        scanf("%d", a + i); --a[i];
        now.insert(1);
    }

    mp[now] = 0;
    dfs(0);
    int nPr = sub(0, inv(n * (n - 1) / 2));
    fill_n(e[0], idx + 1, 0);
    e[0][0] = 1;
    for(int i=1; i<=idx; i++){
        for(int j=0; j<=idx; j++){
            e[i][j] = mul(e[i][j], nPr);
        }
        e[i][i] = 1;
        e[i][idx + 1] = 1;
    }

    assert(solveEq(e, res, idx + 1));
    multiset<int> Ask;
    for(int i=0; i<n; i++)
        if(!vis[i]){
            cnt = 0;
            getsize(i);
            Ask.insert(cnt);
        }

    printf("%d\n", res[mp[Ask]]);
    return 0;
}