题解：AT_yahoo_procon2019_final_e Espionage

· · 题解

困难题。

考虑这个哈密顿回路长什么样子。

我们考虑按行做一下拆分。

之所以不选列做拆分是因为，你可以在一列里面乱跳，而且 m 可以到 10^9。

你发现，按行做拆分后，底下若干行会被完全填满，并且形成若干条“链”。

考虑这些链对下一行的贡献。

首先，链的两端肯定要继续往上延伸，设有 x 条链，那么就有 2x 个位置必定被下面延伸上来的链占据。剩下还有 n-2x 个位置没啥限制，我们称其为孤立点。你发现往上拓展一行后，这些孤立点必须合并进链中，不然要么这个点你遍历不到，要么你的终点不对。然后链之间也可以互相合并。

那我们不妨以“这一行合并完后有多少条链”来进行 dp。

设 dp_{i,j} 表示前 i 行合并完剩 j 条链的方案数。我们设 f_{i,j} 表示把 i 条链，经过合并，剩下 j 条链，不允许出现孤立点未被合并的方案数。那么有转移：

dp_{i,j} = \sum_{k=0}^n dp_{i-1,k} f_{k,j}

那假如我们把 f 求出来了，这个玩意显然就可以矩阵快速幂了对吧。问题变成怎么求 f。

经典二项式反演，设 g_{i,j} 表示把 i 条链，经过合并，剩下 j 条链的方案数。那么有：

g_{i,j} = \frac{i!}{j!} \binom{i-1}{j-1}

相当于你枚举全排列，然后插板，分成 j 段，每段至少包含一个元素。然后依次合并每段中的链。

那么容斥一下有：

f_{i,j} = \sum_{k=0}^n (-1)^{k} \binom{n-2i}{k} g_{n-i-k,j-k}

相当于你枚举有多少个孤立点未被合并。

好的，那你预处理，然后矩阵快速幂即可。

AC Code

#include<bits/stdc++.h>
#define int long long
using namespace std ;
const int MAXN = 205 ;
const int MOD = 1e9 + 7 ;
int val[MAXN] , inv[MAXN] ;

int ksm(int i , int j) {
    int ans = 1 ;
    while (j) {
        if (j & 1) {
            ans = ans * i % MOD ;
        }
        i = i * i % MOD ;
        j >>= 1 ;
    }

    return ans ;
}

int s[MAXN][MAXN] ;

int C(int i , int j) {
    if (i < j or j < 0) return 0 ;
    return val[i] * inv[j] % MOD * inv[i - j] % MOD ;
}

void add(int& u , int v) {
    u += v ;
    if (u > MOD)    u -= MOD ;
}

struct Matrix {
    int n , m ;
    vector<int> a[MAXN] ;

    void InitZero(int N , int M) {
        n = N , m = M ;
        for (int i = 0 ; i <= N ; i ++) {
            a[i].assign(M + 1 , 0) ;
        }
        return ;
    }

    void Init(int N , int M) {
        n = N , m = M ;
        InitZero(N , M) ;
        for (int i = 0 ; i <= min(N , M) ; i ++) {
            a[i][i] = 1 ;
        }
        return ;
    }

    Matrix operator * (const Matrix& l) const {
        Matrix ans ;
        ans.InitZero(n , l.m) ;
        for (int i = 0 ; i <= n ; i ++) {
            for (int j = 0 ; j <= m ; j ++) {
                for (int k = 0 ; k <= l.m ; k ++) {
                    add(ans.a[i][k] , a[i][j] * l.a[j][k] % MOD) ;
                }
            }
        }

        return ans ;
    }
};

Matrix ksm(Matrix i , int j) {
    Matrix ans ;
    ans.Init(i.n , i.m) ;
    while (j) {
        if (j & 1) {
            ans = ans * i ;
        }
        i = i * i ;
        j >>= 1 ;
    }

    return ans ;
}

signed main() {
    //freopen("building.in" , "r" , stdin) ;
    //freopen("building.out" , "w" , stdout) ;

    ios::sync_with_stdio(0) ;
    cin.tie(0) ;
    cout.tie(0) ;

    int n , m ;
    cin >> n >> m ;

    val[0] = 1 ;
    for (int i = 1 ; i <= n ; i ++) {
        val[i] = val[i - 1] * i % MOD ;
    }

    inv[n] = ksm(val[n] , MOD - 2) ;

    for (int i = n - 1 ; i >= 0 ; i --) {
        inv[i] = inv[i + 1] * (i + 1) % MOD ;
    }

    for (int i = 0 ; i <= n ; i ++) {
        for (int j = 0 ; j <= n ; j ++) {
            s[i][j] = val[i] * inv[j] % MOD * C(i - 1 , j - 1) % MOD ;
        }
    }

    Matrix d ;
    d.InitZero(n / 2 , n / 2) ;

    for (int i = 0 ; i <= n / 2 ; i ++) {
        for (int j = 0 ; j <= n / 2 ; j ++) {
            int cnt = 0 ;
            for (int k = 0 ; k <= n - 2 * i ; k ++) {
                if (n - i - k < 0 or j - k < 0) break ;
                int val = C(n - 2 * i , k) * s[n - i - k][j - k] % MOD ; 
                if (k & 1)  add(cnt , (MOD - val)) ;
                else    add(cnt , val) ;
            }

            d.a[i][j] = cnt ;
        }
    }

    Matrix z ;
    z.InitZero(1 , n / 2) ;
    z.a[1][0] = 1 ;

    z = z * ksm(d , m - 1) ;

    int ans = 0 ;
    for (int i = 0 ; i <= n / 2 ; i ++) {
        add(ans , z.a[1][i] * val[n - 1 - i] % MOD) ;
    }
    cout << ans << "\n" ;

    return 0 ;
}