[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;
}