[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
", Ans.Num[C[(1 << k) - 1]][C[(1 << k) - 1]]);
}
 
signed main() {
    _main();
    return 0;
}