CodeChef - CODIE

最新推荐文章于 2019-04-19 23:29:22 发布

原创最新推荐文章于 2019-04-19 23:29:22 发布 · 294 阅读

0 ·

CC 4.0 BY-SA版权

DP 同时被 2 个专栏收录

6 篇文章

订阅专栏

矩阵快速幂

1 篇文章

订阅专栏

%%%

参考自：yts1999’s blog

题目链接 : Codie Bird

sol:

考虑朴素的 $d p$ 方程

$d p [i] [j] = d p [i - 1] [j - 1] + d p [i - 1] [j] + d p [i] [j + 1]$

递推矩阵大概长这样。
$T=\begin{bmatrix} 1 & 1 & 0 & 0 & 0 & 0 & \dots & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 0 & \dots & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 & 0 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & 0 & 0 & 0 & 0 & \dots & 1 & 1\\ \end{bmatrix}$

有障碍的时候，将 $[1, a]$ 和 $[b, m]$ 列的矩阵元素全部置为0。将这个转移矩阵设为 $S$ .
设 $v_i=\begin{bmatrix} dp[i][1] & dp[i][2] & \cdots & dp[i][k]\end{bmatrix}$ ，则 $v_n=v_1(\sum\limits_{i=1}^{n-2}T^{i-1}ST^{n-2-i})T$ 。
设 $P_n = \sum\limits_{i=1}^{n} T^{i-1}ST^{n-i}$ ，则 $v_n = v_1P_{n-2}T$ 。
实现上可以用单位矩阵代替 $v_1$ ， $v_n$ 就是最终矩阵的第 $\frac{k}{2}$ 行。

求 $P_n$ 可以用类似倍增的做法。

$P_1 = S$
若n为偶数， $P_n = P_\frac{n}{2} \cdot T^\frac{n}{2} + T^\frac{n}{2} \cdot P_\frac{n}{2}$
若n为奇数， $P_n = P_{n-1} \cdot T + T^{n-1} \cdot S$

实现上可以考虑 $1010_{(2)} = 10_{(10)}$ 这样的求解过程。

code:

#include<iostream>
#include<cstring>
#include<string>
#include<algorithm>
#include<queue>

using namespace std;
typedef long long ll;
const int maxn = 52;
const int mod = 1e9+7;
int n,m;

struct Mat{
    ll a[maxn][maxn];
    void init(int k){
        memset(a,0,sizeof(a));
        for(int i = 1;i<=k;i++) a[i][i] = 1;
    }

    void reset(){
        memset(a,0,sizeof(a));
    }

    Mat operator + (Mat b){
        Mat c;
        c.reset();
        for(int i = 1;i<=m;i++)
        for(int j = 1;j<=m;j++)
            c.a[i][j] = (a[i][j] + b.a[i][j]) % mod;
        return c;
    }

    Mat operator * (Mat b){
        Mat c; 
        for(int i = 1;i<=m;i++)
        for(int j = 1;j<=m;j++){
            c.a[i][j] = 0;
            for(int k = 1;k<=m;k++){
                c.a[i][j] += a[i][k] * b.a[k][j] % mod;
            }
            c.a[i][j] %= mod;
        }
        return c;
    }

    Mat pow(int b){
        Mat ret; ret.init(m);
        Mat a = (*this);
        while(b){
            if(b&1) ret = ret * a;
            a = a * a;
            b>>=1;
        }
        return ret;
    }
};

ll qpow(ll a,int b){
    ll ret = 1;
    while(b){
        if(b&1) ret = ret * a % mod;
        a = a*a % mod;
        b>>=1;
    }
    return ret;
}

int main(){
    int a,b;
    scanf("%d%d%d%d",&n,&m,&a,&b);
    n -= 2;
    Mat T,S;
    T.reset(); S.reset();
    for(int i = 1;i<=m;i++) 
        T.a[i][i] = T.a[i][i-1] = T.a[i][i+1] = 1;
    
    S = T;
    for(int i = 1;i<=m;i++){
        if(i<=a || i>=b){
            for(int j = 1;j<=m;j++){
                S.a[j][i] = 0;
            }
        }
    }
    
    int bit = 0;
    for(int i = n;i>0;i>>=1) bit++;
    Mat P = S; Mat Ts = T;
    for(int i = bit - 2;i>=0;i--){
        P = P * Ts + Ts * P;
        Ts = Ts * Ts;
        if((n>>i)&1){
            P = P * T + Ts * S;
            Ts = Ts * T;
        }
    }
    P = P * T;
    ll ans = P.a[m/2][m/2];
    ans = ans * qpow(n,mod-2) % mod;
    printf("%lld\n",ans);
    return 0;
}