[BZOJ2004][HNOI2010]Bus 公交線路

題目描述

小Z所在的城市有\(N\)個公交車站，排列在一條長(\(N-1\))km的直線上，從左到右依次編號爲\(1\)到\(N\)，相鄰公交車站間的距離均爲\(1\)km。 做爲公交車線路的規劃者，小Z調查了市民的需求，決定按下述規則設計線路：

1.設共\(K\)輛公交車，則\(1\)到\(K\)號站做爲始發站，\(N-K+1\)到N號臺做爲終點站。

2.每一個車站必須被一輛且僅一輛公交車通過（始發站和終點站也算被通過）。

3.公交車只能從編號較小的站臺駛往編號較大的站臺。

4.一輛公交車通過的相鄰兩個站臺間距離不得超過\(P\)km。

在最終設計線路以前，小Z想知道有多少種知足要求的方案。因爲答案可能很大，你只需求出答案對\(30031\)取模的結果。

Input

僅一行包含三個正整數\(N, K ,P\)，分別表示公交車站數，公交車數，相鄰站臺的距離限制。

\(N<=10^9，1<P<=10，K<N，1<K<=P\)

Output

僅包含一個整數，表示知足要求的方案數對\(30031\)取模的結果。

Sample Input

10 3 3
5 2 3
10 2 4

Sample Output

1
3
81

看到\(k,p\)如此之小天然就是狀壓\(dp\)了。

咱們能夠十分想到\(O((n-k)*p*2^p)\)的\(dp\)，即\(dp[i][j]\)表示第\(i\)公里，第\(i-p+1\)到\(i\)公里是否有車的狀態爲\(j\)的方案數。

轉移是很顯然的枚舉一個\(1\)把它放到後面去，因爲直接轉移的話會算重複，因此咱們要限制一位。

咱們只要限制最低位必定要有一個\(1\)在，這要的話就能夠保證狀態的轉移不會重複。

\(dp\)部分見下：

dp[k][(1 << k) - 1] = 1;
ret(i, k, n) drep(K, (1 << P) - 1, 1) if (dp[i][K]) {
    if (!(K & 1))continue;
    ret(p, 0, P) if (K & 1 << p) {
        int Now = (K ^ (1 << p)) % (1 << P - 1);
        Now <<= 1, Now |= 1;
        dp[i + 1][Now] += dp[i][K], Mod(dp[i + 1][Now]);
    }
}

可是！！！這是過不了的。

咱們發現每次dp的轉移都是同樣的

每次都是一個狀態轉移到另外一個固定的狀態，和\(i\)沒有關係。

因而咱們想到了矩陣加速，利用矩陣加速這個遞推。

還沒完！！！

若是把全部狀態都算上去的話，矩陣的大小爲\(1024*1024\)，每次矩陣乘法的複雜度爲\(O(1024*1024*1024)\)沒法承受。

可是咱們有一個限制條件，即狀態\(i\)必須保證\(i\)的最低位爲\(1\)，這樣的話咱們優化到了\(O(512*512*512)\)的複雜度。

好像仍是有一點難過掉，怎麼辦呢？

咱們一開的的狀態\(1\)的個數爲\(k\)，而轉移時並不會改變\(i\)的\(1\)的個數，因此咱們只要把全部二進制位中\(1\)的個數爲\(k\)的狀態利用矩陣轉移便可。

那麼矩陣的大小爲多少呢，爲\(C(p-1,k)\)，最大爲\(C(9,4)=126\)，這樣的話一次矩陣乘法的複雜度爲\(O(126*126*126)\)，咱們就能夠過掉這道題了。

#include <bits/stdc++.h>
 
using namespace std;
 
#define LL long long
#define reg register
#define clr(a,b) memset(a,b,sizeof a)
#define Mod(x) (x>=mod)&&(x-=mod)
#define abs(a) ((a)<0?-(a):(a))
#define debug(x) cerr<<#x<<"="<<x<<endl;
#define debug2(x,y) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<endl;
#define debug3(x,y,z) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl;
#define rep(a,b,c) for(reg int a=(b),a##_end_=(c); a<=a##_end_; ++a)
#define ret(a,b,c) for(reg int a=(b),a##_end_=(c); a<a##_end_; ++a)
#define drep(a,b,c) for(reg int a=(b),a##_end_=(c); a>=a##_end_; --a)
#define erep(i,G,x) for(int i=(G).Head[x]; i; i=(G).Nxt[i])
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize(3,"Ofast","inline")
 
inline int Read(void) {
    int res = 0, f = 1;
    char c;
    while (c = getchar(), c < 48 || c > 57)if (c == '-')f = 0;
    do res = (res << 3) + (res << 1) + (c ^ 48);
    while (c = getchar(), c >= 48 && c <= 57);
    return f ? res : -res;
}
 
template<class T>inline bool Min(T &a, T const&b) {return a > b ? a = b, 1 : 0;}
template<class T>inline bool Max(T &a, T const&b) {return a < b ? a = b, 1 : 0;}
 
const int N = 3e2 + 5, M = 1e5 + 5, K = 10, mod = 30031;
 
bool MOP1;
 
int n, k, P, cnt, Sz[1 << K], B[130], C[1 << K];
 
struct Matrix {
    int Num[N][N];
    inline void clear(void) {clr(Num, 0);}
    inline void Init(void) {rep(i, 1, cnt)Num[i][i] = 1;}
    inline Matrix operator*(Matrix _)const {
        Matrix Ans;
        Ans.clear();
        rep(i, 1, cnt)rep(j, 1, cnt)rep(k, 1, cnt)Ans.Num[i][j] = (Ans.Num[i][j] + 1ll * Num[i][k] * _.Num[k][j]) % mod;
        return Ans;
    }
} us;
 
inline Matrix qpow(Matrix A, int k) {
    Matrix res;
    res.clear(), res.Init();
    while (k) {
        if (k & 1)res = res * A;
        A = A * A, k >>= 1;
    }
    return res;
}
 
bool MOP2;
 
void _main(void) {
    n = Read(), k = Read(), P = Read();
    ret(i, 1, 1 << P) {
        Sz[i] = Sz[i ^ (i & -i)] + 1;
        if (!(i & 1))continue;
        if (Sz[i] == k)B[++cnt] = i, C[i] = cnt;
    }
    Matrix res;
    res.clear();
    rep(i, 1, cnt) {
        int K = B[i];
        ret(p, 0, P) if (K & 1 << p) {
            if (k - P > p)continue;
            int Now = ((K ^ (1 << p)) & ((1 << P - 1) - 1)) << 1 | 1;
            res.Num[i][C[Now]] = 1;
        }
    }
    Matrix Ans = qpow(res, n - k);
    printf("%lld\n", Ans.Num[C[(1 << k) - 1]][C[(1 << k) - 1]]);
}
 
signed main() {
    _main();
    return 0;
}